912算法模板。 适用于机考参考之用。
本文集成于:https://github.com/stellarkey/912_project
因为模板整理时比较仓促,本文中大量代码为搬运,并非原创,特此说明。
References ACM在线模版-f-zyj
ACM Standard Code Library
https://github.com/LittileNiQ/ACM
https://github.com/eecrazy/ACM
ACM卡常数(各种玄学优化)
……etc.
“考试期间严禁使用任何个人电子设备,可携带并查阅任何纸质材料”
912 SCL
通用 I/O模板 Fast_I/O 当大量IO时,一般用scanf
/printf
进行输入输出即可解决。
输入封装 避免重复书写scanf()
。
1 2 3 4 5 int readint () { int x; scanf ("%d" , &x); return x; } vector <int > vc;vc.pushback(readint());
手写I/O优化 手写readint()
的高速输入法。
1 2 3 4 5 6 inline int freadint () { int x=0 ,f=1 ;char ch=getchar(); while (!isdigit (ch)){if (ch=='-' ) f=-1 ;ch=getchar();} while (isdigit (ch)){x=x*10 +ch-48 ;ch=getchar();} return x*f; }
或模板型:
1 2 3 4 5 6 7 8 9 template <class T >T read (T &x ){ x=0 ;char ch=0 ; while (ch < '0' || ch > '9' )ch=getchar(); while (ch >= '0' && ch <= '9' ){x=x*10 +ch-'0' ;ch=getchar();} return x; } int n=read (n),p=read (p),k=read (k);
I/O同时优化:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 int read () { int x = 0 , w = 1 ; char ch = 0 ; while (ch < '0' || ch > '9' ) { if (ch == '-' ) w = -1 ; ch = getchar(); } while (ch >= '0' && ch <= '9' ) { x = x * 10 + (ch - '0' ); ch = getchar(); } return x * w; } inline void write (int x) { static int sta[35 ]; int top = 0 ; do { sta[top++] = x % 10 , x /= 10 ; } while (x); while (top) putchar (sta[--top] + 48 ); } int num = read (); write (num);
光速I/O 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 namespace IO {const int MAXSIZE = 1 << 20 ;char buf[MAXSIZE], *p1, *p2;#define gc() \ (p1 == p2 && (p2 = (p1 = buf) + fread(buf, 1 , MAXSIZE, stdin ), p1 == p2) \ ? EOF \ : *p1++) inline int rd () { int x = 0 , f = 1 ; char c = gc(); while (!isdigit (c)) { if (c == '-' ) f = -1 ; c = gc(); } while (isdigit (c)) x = x * 10 + (c ^ 48 ), c = gc(); return x * f; } char pbuf[1 << 20 ], *pp = pbuf;inline void push (const char &c) { if (pp - pbuf == 1 << 20 ) fwrite(pbuf, 1 , 1 << 20 , stdout ), pp = pbuf; *pp++ = c; } inline void write (int x) { static int sta[35 ]; int top = 0 ; do { sta[top++] = x % 10 , x /= 10 ; } while (x); while (top) push(sta[--top] + '0' ); } }
I/O重定向 1 2 3 4 5 6 #define _OJ_ #ifndef _OJ_ freopen("OJ.in" , "r" , stdin ); freopen("OJ.out" , "w" , stdout ); #endif
多组测试数据 1 2 3 4 5 6 7 8 9 10 11 int n; while (scanf ("%d" , &n) != EOF){ } while (cin >> n != NULL ){ } while (gets(STRING)){ }
循环宏 1 2 #define _for(i,a,b) for(int i =(a); i < (b); ++i) #define _rep(i,a,b) for(int i =(a); i <= (b); ++i)
C++ I/O模板 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 #include <bits/stdc++.h> using namespace std ;#define _OJ_ #define MAXN 2000005 #define INF 1000000007 #define _for(i,a,b) for(int i =(a); i < (b); ++i) #define _rep(i,a,b) for(int i =(a); i <= (b); ++i) typedef long long ll;typedef pair<string , string > pss;typedef pair<int , int > pii;typedef pair<string , int > psi;const double PI = acos (-1.0 );const double eps = 1e-6 ;const int MOD = 1e7 + 7 ;int main () { ios_base::sync_with_stdio(false ); cin .tie(NULL ); #ifndef _OJ_ freopen("OJ.in" , "r" , stdin ); freopen("OJ.out" , "w" , stdout ); #endif int n, m; cin >> n >> m; vector <int > A; int tmp; _for(i, 0 , n) {cin >> tmp; A.push_back(tmp); } int a[MAXN]; _for(i, 0 , n) {cin >>a[i];} #ifndef _OJ_ cout <<"A: " << A <<endl ; _for(i, 0 , n) cout << a[i] << " \n" [i == n - 1 ]; _rep(i, 1 , n) cout << a[i] << " \n" [i == n]; #endif return 0 ; }
Java I/O模板 1 2 3 4 5 6 7 8 9 10 11 import java.io.*;import java.util.*;import java.math.BigInteger; import java.math.BigDecimal; public class Main { public static void main (String args[]) throws Exception { Scanner cin=new Scanner(System.in); int a = cin.nextInt(), b = cin.nextInt(); System.out.println(a+b); } }
Python I/O模板 1 2 3 4 5 6 7 8 9 10 a, b = map(int, input().split(' ' )) print(a + b) import syssys.stdin = open("OJ.in" , mode='r' ) sys.stdout = open("OJ.out" , mode='w' ) A = list(map(int, input().split())) print(A)
1 2 3 4 import sysif __name__ == '__main__' : a, b = sys.stdin.readline().split(' ' ) print(int(a) + int(b))
多组测试数据:
1 2 3 4 5 6 while True : try : a, b = map(int, raw_input().strip().split()) print(a + b) except EOFError: break
复杂度适定表 时间限制(千万)
参考时间一般为1000ms
。
一般而言,可行的算法复杂度在数字意义上不能超过一千万 ,即$10^7$。
数据范围
允许算法复杂度
常用思路
$\infty$
$O(logn)$
二分、数论、打表
$10^7\longrightarrow10^8$
$O(n)$
在线、减治
$10^5\longrightarrow10^6$
$O(nlogn)$
排序、分治
$10^3\longrightarrow10^4$
$O(n^2)$
动态规划
$200\longrightarrow500$
$O(n^3)$
模拟
$20\longrightarrow50$
$O(2^n)$
BFS/DFS
$12$
$O(n!)$
组合
空间限制(千万) 参考内存一般为256MB
。大约相当于64M个int
(可开千万 级数组 ),其余类型可换算。
机试中通常不对空间做特别要求,因此常常适用“空间换时间”。
栈溢出 递归过深/局部数组达百万级常引发爆栈,此时需改写方法,或手动加栈:
1 2 #pragma comment(linker, "/STACK:1024000000,1024000000" )
C++头文件 万能头 bits/stdc++.h
常规头
头文件
说明
iostream
数据流输入/输出,如cin,cout等
string
字符串类
vector
STL 动态数组容器
queue
STL 队列容器
set
STL 集合容器
stack
STL 堆栈容器
algorithm
STL 通用算法,如find等
deque
STL 双端队列容器
list
STL 线性列表容器
map
STL 映射容器
iterator
STL 迭代器
utility
STL 通用模板类,如pair,make_pair等
bitset
STL 位集容器
stdexcept
标准异常类
new
动态内存分配
memory
STL通过分配器进行的内存分配,如auto_ptr等
numeric
STL常用的数字操作,如accumulate等
random
随机数
regex
正则表达式
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 #include <cstdio> #include <cstring> #include <algorithm> #include <iostream> #include <string> #include <vector> #include <stack> #include <bitset> #include <cstdlib> #include <cmath> #include <set> #include <list> #include <deque> #include <map> #include <queue> #include <functional> // less/greater等
模拟
模拟题一般不考察效率,因此在策略上应倾向于采用Python解决。
枚举 子集枚举 1 2 3 4 5 6 7 8 9 10 11 12 void print_subset (int n, bool * B, int cur) { if (cur == n) { for (int i = 0 ; i < cur; ++i) if (B[i]) printf ("%d " , i); printf ("\n" ); return ; } B[cur] = 1 ; print_subset(n ,B, cur+1 ); B[cur] = 0 ; print_subset(n ,B, cur+1 ); }
亦可采用STL中的 bitset
位图。
二进制法:
1 2 3 4 5 6 7 8 void print_subset (int n, int s) { for (int i = 0 ; i < n; ++i) if (s & (1 <<i)) printf ("%d " , i); printf ("\n" ); } for (int i = 0 ; i < (1 << n); ++i) print_subset(n, i);
回溯法 即剪枝版的深度优先搜索。
1 2 3 4 5 6 7 8 9 10 11 12 13 void DFS (int cur) { if (){ } else { for (...){ EXECUTE(i); vis[i] = true ; dfs(cur+1 ); vis[i] = false ; } } }
扫描法 扫描法往往在枚举时维护一些重要的量,从而简化计算。例如,从左到右考虑数组的各个元素,也可以说从左到右“扫描”。
尺取法/滑动窗口 $O(n)$ 端点单调的双指针。
例:输入一个序列,找到最长连续子序列,保证其内部没有相同元素。
利用滑动窗口+set
。set
的单次插/查/删都是$O(logn)$。总复杂度为$O(nlogn)$。
也可以利用map
预处理,利用序列A处理序列B,B中元素是“A中当前元素的上一个相同元素的下标”。总复杂度为$O(nlogn+n)=O(nlogn)$。
众数 $O(n)$ 1 2 3 4 5 6 7 8 9 int majEleCandidate (vector <int > A) { int maj; for (int c = 0 , i = 0 ; i < A.size (); ++i) if ( 0 == c){ maj = A[i]; c = 1 ; } else maj == A[i] ? c++ : c--; return maj; }
二分法 $O(logn)$ 见到“最大最小 ”同时出现时,考虑二分法。二分法的基本思路:“猜 ”,如果能够事先确定答案的范围,那么就可以二分来不断猜测正确答案,利用轴点是否能满足某个条件来缩减区间。
例:把一个长度为N的序列A划分成k个连续子序列,使得每个子序列的和S[i]
中的最大值尽量小。
二分范围:$[0, \;\sum A[i]\;]$。总复杂度为$O(Nlog\{\sum A[i]\})$。
每轮,判断对于轴点x
,是否存在划分,使得max{s[i]} <= x
?这不难判断。只需扫描一遍,每次尽量往右划分即可。
二分剪枝 有时,可以提前确定并缩小二分法的范围。(如19年真题-第2题-高精度开根)
通过对二分的上下界的准确定界预处理,使得效率增加。(当M越大,效果越好 )
二分剪枝$O(\cfrac{1}{M}logN)$ $\times$ 高精度乘法$O((logN)^2)$ $\times$ 快速幂$O(logM)$。(Python自带快速幂)
a**(1/b)
表示a的b次根,等价于pow(a, 1/b)
,存在精度误差。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 m = int(input()) n = int(input()) le, ri = 0 , 1 while pow(ri, m) < n: le, ri= ri, 2 *ri print(le, ri) while le < ri - 1 : mid = (le + ri) >> 1 t = pow(mid, m) if t <= n: le = mid else : ri = mid print(le)
折半枚举(中途相遇法) $O(\sqrt M)$ 对n
个自由变量的组合问题,可以化为两个n/2
的子问题分步求解。
如,给定4个集合,从中依次抽取4个元素a,b,c,d
,使得a+b+c+d=0
,问:有多少种选法?
先计算a+b
耗费$O(n^2)$并保存到新集M中;同理计算-(c+d)
的组合耗费$O(n^2)$保存到新集N。此时只需再按照N中元素查找M集即可,耗费$O(n^2logn)$。
模拟退火 如果能够给出状态的能量函数(能量越低状态越优秀),则可以采用模拟退火。
有n个重物,每个重物系在一条足够长的绳子上。每条绳子自上而下穿过桌面上的洞,然后系在一起。图中X处就是公共的绳结。假设绳子是完全弹性的(不会造成能量损失),桌子足够高(因而重物不会垂到地上),且忽略所有的摩擦。问绳结X最终平衡于何处。
文件的第一行为一个正整数n(1≤n≤1000),表示重物和洞的数目。接下来的n行,每行是3个整数:Xi.Yi.Wi,分别表示第i个洞的坐标以及第 i个重物的重量。(-10000≤x,y≤10000, 0<w≤1000 )
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 #include <bits/stdc++.h> #define down 0.996 using namespace std ;int n;struct node {int x;int y;int w;}object[2005 ]; double ansx,ansy,answ;double energy (double x,double y) { double r=0 ,dx,dy; for (int a=1 ;a<=n;a++) { dx=x-object[a].x; dy=y-object[a].y; r+=sqrt (dx*dx+dy*dy)*object[a].w; } return r; } void sa () { double t=3000 ; while (t>1e-15 ) { double ex=ansx+(rand()*2 -RAND_MAX)*t; double ey=ansy+(rand()*2 -RAND_MAX)*t; double ew=energy(ex,ey); double de=ew-answ; if (de<0 ) { ansx=ex; ansy=ey; answ=ew; } else if (exp (-de/t)*RAND_MAX>rand()) { ansx=ex; ansy=ey; } t*=down; } } void solve () { sa(); sa(); sa(); sa(); } int main () { cin >>n; for (int a=1 ;a<=n;a++) { scanf ("%d%d%d" ,&object[a].x,&object[a].y,&object[a].w); ansx+=object[a].x; ansy+=object[a].y; } ansx/=n; ansy/=n; answ=energy(ansx,ansy); solve(); printf ("%.3lf %.3lf\n" ,ansx,ansy); return 0 ; }
语法解析 表达式求值
Python支持自动表达式求值:print(eval(input()) % 10000)
。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 #include <bits/stdc++.h> #define For(i,l,r) for(int i=(l);i<=(r);i++) using namespace std ; const int prior[7 ][7 ]={{ 1 , 1 ,-1 ,-1 ,-1 ,1 ,-1 }, { 1 , 1 ,-1 ,-1 ,-1 ,1 ,-1 }, { 1 , 1 , 1 , 1 ,-1 ,1 ,-1 }, { 1 , 1 , 1 , 1 ,-1 ,1 ,-1 }, {-1 ,-1 ,-1 ,-1 ,-1 ,0 ,-1 }, { 0 , 0 , 0 , 0 , 0 ,0 , 0 }, { 1 , 1 , 1 , 1 ,-1 ,1 , 1 }}; char strh[2000 ]={0 },strz[2000 ]={0 }; char opst[10000 ],top=0 ; int numst[10000 ],t=0 ;int f (char op) { switch (op){ case '+' : return 0 ; case '-' : return 1 ; case '*' : return 2 ; case '/' : return 3 ; case '(' : return 4 ; case ')' : return 5 ; case '^' : return 6 ; } } void convert (char strz[],char strh[]) { int len=strlen (strz); int j=0 ; opst[++top]='(' ; For(i,0 ,len-1 ){ if (isdigit (strz[i])) strh[j++]=strz[i]; else { while (prior[f(opst[top])][f(strz[i])]>0 ) strh[j++]=opst[top--]; if (prior[f(opst[top])][f(strz[i])]==0 ) top--; else opst[++top]=strz[i]; } } } void print (char s[]) { int len=strlen (s); For(i,0 ,len-1 ){ printf ("%c " ,s[i]); } printf ("\n" ); } void cal (char strh[]) { int len=strlen (strh); print (strh); For(i,0 ,len-1 ) if (isdigit (strh[i])) numst[++t]=strh[i]-'0' ; else { t--; switch (strh[i]){ case '+' :numst[t]=numst[t]+numst[t+1 ];break ; case '-' :numst[t]=numst[t]-numst[t+1 ];break ; case '*' :numst[t]=numst[t]*numst[t+1 ];break ; case '/' :numst[t]=numst[t]/numst[t+1 ];break ; case '^' :numst[t]=pow (numst[t],numst[t+1 ]);break ; } For(j,1 ,t) printf ("%d " ,numst[j]); print (strh+i+1 ); } } int main () { gets(strz); int len=strlen (strz); strz[len]=')' ; strz[len+1 ]='\0' ; convert(strz,strh); cal(strh); return 0 ; }
字符串 KMP算法 $O(n)$ 模式串匹配。多个模式串时采用AC自动机。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 int next[MAXN];int KMP (string text , string pattern) { int i = next[0 ] = -1 , j = 0 ; int n = text .size (), m = pattern.size (); while (j < m){ if (i == -1 || pattern[j] == pattern[i]){ i++; j++; next[j] = i; } else { i = next[i]; } } i = j = 0 ; while (i < n && j < m) { if (j == -1 || text [i] == pattern[j]){ i++; j++; } else { j = next[j]; } } if (j == m) return i - j + 1 ; else return -1 ; }
Z函数/扩展KMP $O(n)$ 对于一般的 KMP 只需要求所有 extend[i] == m
的位置,那么 ExKMP 就是需要求出这个 extend[i]
数组。
求解T
与S
对的每一个后缀的最长公共前缀 LCP。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 #include <bits/stdc++.h> #define For(i, l, r) for (register int i = (l), i##end = (int)(r); i <= i##end; ++i) #define Fordown(i, r, l) for (register int i = (r), i##end = (int)(l); i >= i##end; --i) #define Rep(i, r) for (register int i = (0), i##end = (int)(r); i < i##end; ++i) #define Set(a, v) memset(a, v, sizeof(a)) #define Cpy(a, b) memcpy(a, b, sizeof(a)) #define debug(x) cout << #x << ": " << (x) << endl #define next Next using namespace std ;template <typename T> inline bool chkmin (T &a, T b) { return b < a ? a = b, 1 : 0 ; }template <typename T> inline bool chkmax (T &a, T b) { return b > a ? a = b, 1 : 0 ; }inline int read () { int x(0), sgn(1); char ch(getchar()); for (; !isdigit (ch); ch = getchar()) if (ch == '-' ) sgn = -1 ; for (; isdigit (ch); ch = getchar()) x = (x * 10 ) + (ch ^ 48 ); return x * sgn; } void File () {#ifdef zjp_shadow freopen ("1461.in" , "r" , stdin ); freopen ("1461.out" , "w" , stdout ); #endif } const int N = 1e6 + 1e3 ;void Get_Next (char *S, int *next) { int lenS = strlen (S + 1 ), p = 1 , pos; next[1 ] = lenS; while (p + 1 <= lenS && S[p] == S[p + 1 ]) ++ p; next[pos = 2 ] = p - 1 ; For (i, 3 , lenS) { int len = next[i - pos + 1 ]; if (len + i < p + 1 ) next[i] = len; else { int j = max (p - i + 1 , 0 ); while (i + j <= lenS && S[j + 1 ] == S[i + j]) ++ j; p = i + (next[pos = i] = j) - 1 ; } } } void ExKMP (char *S, char *T, int *next, int *extend) { int lenS = strlen (S + 1 ), lenT = strlen (T + 1 ), p = 1 , pos; while (p <= lenT && S[p] == T[p]) ++ p; p = extend[pos = 1 ] = p - 1 ; For (i, 2 , lenS) { int len = next[i - pos + 1 ]; if (len + i < p + 1 ) extend[i] = len; else { int j = max (p - i + 1 , 0 ); while (i + j <= lenS && j <= lenT && T[j + 1 ] == S[i + j]) ++ j; p = i + (extend[pos = i] = j) - 1 ; } } } char S[N], T[N]; int next[N], extend[N]; int main () { File (); scanf ("%s" , S + 1 ); scanf ("%s" , T + 1 ); Get_Next(T, next); ExKMP(S, T, next, extend); For (i, 1 , strlen (S + 1 )) printf ("%d%c" , extend[i], i == iend ? '\n' : ' ' ); return 0 ; }
Manacher算法 $O(n)$
求最长回文子串。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 char Ma[MAXN*2 ];int Mp[MAXN*2 ];void Manacher (char s[],int len) { int l=0 ; Ma[l++]='$' ; Ma[l++]='#' ; for (int i=0 ;i<len;i++){ Ma[l++]=s[i]; Ma[l++]='#' ; } Ma[l]=0 ; int mx=0 ,id=0 ; for (int i=0 ;i<l;i++){ Mp[i]=mx>i?min (Mp[2 *id − i],mx − i):1 ; while (Ma[i+Mp[i]]==Ma[i − Mp[i]])Mp[i]++; if (i+Mp[i]>mx){ mx=i+Mp[i]; id=i; } } } char s[MAXN];int main () { while (scanf ("%s" ,s)==1 ){ int len=strlen (s); Manacher(s,len); int ans=0 ; for (int i=0 ;i<2 *len+2 ;i++) ans=max (ans,Mp[i] − 1 ); printf ("%d\n" ,ans); } return 0 ; }
最小表示法 $O(n)$ 求解与字符串S
循环同构的最小字典序的字符串。
如对于S = deabc
,最小表示位置为2
(a
),最小字符串为abcde
。
1 2 3 4 5 6 7 8 9 10 11 12 13 int min_cyclic_string (string sec) { int k = 0 , i = 0 , j = 1 , n = sec.size (); while (k < n && i < n && j < n) { if (sec[(i + k) % n] == sec[(j + k) % n]) { k++; } else { sec[(i + k) % n] > sec[(j + k) % n] ? i = i + k + 1 : j = j + k + 1 ; if (i == j) i++; k = 0 ; } } return min (i, j); }
Lyhdon分解法 $O(n)$ Lyndon 串:对于字符串 S,如果 S 的字典序严格小于 S 的所有后缀的字典序,我们称 S 是简单串,或者 Lyndon 串 。举一些例子,a
, b
, ab
, aab
, abb
, ababb
, abcd
都是 Lyndon 串。如果 S 是 Lyndon 串,当且仅当 S 的字典序严格小于它的所有非平凡的循环同构串。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 string min_cyclic_string (string s) { s += s; int n = s.size (); int i = 0 , ans = 0 ; while (i < n / 2 ) { ans = i; int j = i + 1 , k = i; while (j < n && s[k] <= s[j]) { if (s[k] < s[j]) k = i; else k++; j++; } while (i <= k) i += j - k; } return s.substr(ans, n / 2 ); }
Lyhdon分解 Lyndon 分解:串 S 的 Lyndon 分解(划分)记为W1 W2 ... Wk
,其中所有 Wi
为简单串,并且他们的字典序按照非严格单减排序,即 W1 >= W2 >= ... >= Wk
。可以发现,这样的分解存在且唯一。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 vector <string > duval (string const & s) { int n = s.size (), i = 0 ; vector <string > factorization; while (i < n) { int j = i + 1 , k = i; while (j < n && s[k] <= s[j]) { if (s[k] < s[j]) k = i; else k++; j++; } while (i <= k) { factorization.push_back(s.substr(i, j - k)); i += j - k; } } return factorization; } int main () { string s = "fghacdabcde" ; vector <string > factorization = duval(s); for (int i = 0 ; i < factorization.size (); ++i) cout << factorization[i] <<" \n" [i == factorization.size () - 1 ]; return 0 ; }
日期 给定日期,
week()
:计算星期几
kthDay()
:计算是该年的第几天
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 struct Date { int year, month, day; Date(int y, int m, int d): year(y), month(m), day(d) {} int week () { int y = year, m = month, d = day; if (m==1 ||m==2 ) m+=12 ,y--; int w = (d+2 *m+3 *(m+1 )/5 + y + y/4 - y/100 + y/400 ) % 7 ; return ++w; } int MONTH[13 ] = {0 , 31 , 28 , 31 , 30 , 31 , 30 , 31 , 31 , 30 , 31 , 30 , 31 }; int kthDay () { if (year % 4 == 0 && (year % 100 || year % 400 == 0 )) MONTH[2 ] = 29 ; int days = 0 ; _for(i, 1 , month) days += MONTH[i]; return days += day; } };
排序
尽量借助STL
实现排序。
快速排序 $\Theta(nlogn)+O(1)$ LUG分类快排 1 2 3 4 5 6 7 8 9 10 11 12 void quicksort (vector <int > &a, int le, int ri) { if (le >= ri - 1 ) return ; int i = le - 1 , j = le - 1 , pivot = ri - 1 , r = (rand()%(ri-le))+le; swap(a[pivot], a[r]); while (++j < pivot){ if (a[j] < a[pivot]) swap(a[++i], a[j]); } swap(a[i + 1 ], a[pivot]); quicksort(a, le, i+1 ); quicksort(a, i+2 , ri); }
LUG双向快排 1 2 3 4 5 6 7 8 9 10 void quicksort (vector <int > &a, int le, int ri) { int i = le, j = ri - 1 , Pivot = a[(rand()%(ri-le))+le]; while (i <= j){ while (a[i]<Pivot) i++; while (a[j]>Pivot) j--; if (i<=j) swap(a[i++],a[j--]); } if (le < j ) quicksort(a, le, j+1 ); if (i < ri-1 ) quicksort(a, i, ri); }
二分查找 $O(logn)$
亦可采用STL算法,如:binary_search(beg, end, e, cmp)
/ lower_bound(a, a+a.size(), x, cmp)
等。
1 2 3 4 5 6 7 8 9 10 11 int binarySearch (vector <int > &a, int e) { int le = 0 , ri = a.size (); while (ri - le > 1 ){ int mid = (le + ri) >> 1 ; if (a[mid] <= e) le = mid; else ri = mid; } return le; }
归并排序 $O(nlogn)+O(n)$ 1 2 3 4 5 6 7 8 9 10 11 12 13 void mergesort (vector <int >& a, int le, int ri) { if (le >= ri - 1 ) return ; int mid = (le + ri) >> 1 ; mergesort(a, le, mid); mergesort(a, mid, ri); vector <int > b; int p = le, q = mid; while (p < mid && q < ri) if (a[p] <= a[q]) b.push_back(a[p++]); else b.push_back(a[q++]); while (p < mid) b.push_back(a[p++]); while (q < ri) b.push_back(a[q++]); for (int i = le; i < ri; ++i) a[i] = b[i-le]; }
STL写法:(或者直接调用stable_sort()
)
1 2 3 4 5 6 void mergesort (vector <int >& a, int le, int ri) { if (le >= ri - 1 ) return ; int mid = (le + ri) >> 1 ; mergesort(a, le, mid); mergesort(a, mid, ri); inplace_merge(a.begin () + le, a.begin () + mid, a.begin ()+ ri); }
求逆序对 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 long long mergesort (vector <int >& a, long long le, long long ri) { long long ans = 0 ; if (le >= ri - 1 ) return 0 ; long long mid = (le + ri) >> 1 ; ans += mergesort(a, le, mid); ans += mergesort(a, mid, ri); vector <long long > b; long long p = le, q = mid; long long count = 0 ; while (p < mid && q < ri){ if (a[p] <= a[q]) b.push_back(a[p++]); else { b.push_back(a[q++]); count += (mid - p); } } while (p < mid) b.push_back(a[p++]); while (q < ri) b.push_back(a[q++]); for (long long i = le; i < ri; ++i) a[i] = b[i-le]; ans += count; return ans; }
动态规划
DP问题的复杂度由两部分组成:状态数$\times$状态转移复杂度。
动态规划的难点在于选取合适的状态 。
做题时明确三个核心:状态方程设计 、边界条件 、递推顺序和约束 。
动态规划的特点:最优子结构 + 子问题重复。
DP有两种状态转移策略:
【刷表法】:更新dp[i]
影响的所有后继(必须保证每个状态依赖的所有影响相互独立)
【填表法】:利用dp[i]
的所有前驱,一次性状态转移计算出dp[i]
贪心 区间问题 不相交区间
区间选点
区间覆盖
线性基 给定n个整数(数字可能重复),求在这些数中选取任意个【不一定连续】,使得他们的异或和最大。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 #include <cstdio> #define maxn 105 typedef long long Lovelive;Lovelive base[maxn]; int main () { Lovelive n, tmp; scanf ("%lld" , &n); for (int i = 0 ; i < n; ++i) { scanf ("%lld" , &tmp); for (int j = 50 ; j >= 0 ; --j) { if (tmp & (1L L << j)) { if (!base[j]) base[j] = tmp; tmp ^= base[j]; } } } Lovelive ans = 0 ; for (int i = 50 ; i >= 0 ; --i) { if (ans < (ans ^ base[i])) ans ^= base[i]; } printf ("%lld\n" , ans); }
DAG模型
很多问题可以转化为DAG上的最长路、最短路或路径计数问题。
对于具有二元关系 的特征,可以用图来建模,即用一条有向边来对应表示一个关系。(这样我们便可以先预处理出一个DAG,将问题转化为DAG模型的求解问题)
DAG最短路/最长路
先将图建立好,假设用邻接矩阵保存在矩阵G中。
记忆化搜索求解DAG $O(e)$ 1 2 3 4 5 6 7 8 9 10 11 int d[MAXN];bool vis[MAXN]; memset (vis, 0 , sizeof (vis)); int DP (int i) { int & ans = d[i]; if (vis[i]) return ans; ans = 1 ; vis[i] = true ; for (int j = 1 ; j <= n; ++j) if (G[i][j]) ans = max (ans, DP(j) + 1 ); return ans; }
多段图 多段图是特殊的DAG,其结点可以划分成若干个阶段 ,每个阶段只由上一个阶段所决定。
多阶段决策问题,阶段定义了天然的计算顺序,因此更适宜采用递推法。如,背包问题。
背包问题(链式结构DP) 01背包 $O(nV) + O(V)$ dp[i][j]
表示总体积不超过j
时前i
个物品所能达到的最大价值。第i
件物品体积为w[i]
,价值为v[i]
。
转移时注意保证j-w[i]
非负,若为负值不能进行转移。
1 2 3 4 5 6 7 int dp[2 ][MAXN]; int DP (vector <int > w, vector <int > v, int n, int V) { for (int i = 0 ; i < n; ++i) for (int j = V; j >= 0 ; --j) dp[(i+1 )&1 ][j] = (j < w[i]) ? (dp[i&1 ][j]) : max (dp[i&1 ][j], dp[i&1 ][j-w[i]] + v[i]); return dp[n&1 ][V]; }
降维版本:(DAG模型 )
1 2 3 4 5 6 7 int dp[MAXN]; int DP (vector <int > w, vector <int > v, int n, int V) { for (int i = 0 ; i < n; ++i) for (int j = V; j >= w[i]; --j) dp[j] = max (dp[j], dp[j-w[i]] + v[i]); return dp[V]; }
恰好装满 背包版本:
dp[i][j]
表示总体积恰好为j
时前i
个物品所能达到的最大价值。
初始值改变:dp[0][0] = 0
,dp[0][j] = -INF/不存在
。(其他步骤不变)
完全背包 $O(nV) + O(V)$ 每个物品均有无限个。
1 2 3 4 5 6 7 int dp[MAXN]; int DP (vector <int > w, vector <int > v, int n, int V) { for (int i = 0 ; i < n; ++i) for (int j = w[i]; j <= V; ++j) dp[j] = max (dp[j], dp[j-w[i]] + v[i]); return dp[V]; }
01背包每件物品最多选择一次,因此倒序保证更新时前缀总是没有放入过i
物品;完全背包每件物品可选无数次,因此正序更新反而恰好满足题意。
多重背包 $O(V\times \sum \text{log}_2k_i)$ 每个物品有有限个。第i
种物品体积为w[i]
,价值为v[i]
,个数为k[i]
。
(1)直接当作01背包处理:$O(V*\sum k_i)$
(2)二进制拆分+01背包:$O(V\times \sum \text{log}_2k_i)$
每种物品拆分为多组,每组包含物品个数为:1,2,4,…,k-2^c+1;其中c为使得k-2^c+1>0的最大整数。对这些若干组物品的不同组合,可以得到0-k之间任意件物品的价值重量和。对分组后的物品进行01背包,即可得到多重背包的解。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 int __v[MAXN], __w[MAXN], cnt = 0 ; for (int i = 0 ; i < n; ++i){ int v, w, k; scanf ("%d%d%d" , &w, &v, &k); int c = 1 ; while (k - c > 0 ){ k -= c; __w[cnt] = c*w; __v[cnt++] = c*v; c >>= 1 ; } __w[cnt] = k*w; __v[cnt++] = k*v; }
线性结构DP 最大子段和 $O(n)$ d[i]
表示以A[i]
结尾的连续序列的最大和。
1 2 3 4 5 6 7 8 9 int dp[MAXN];int MCS (vector <int >& A) { int ret = dp[0 ] = A[0 ]; for (int i = 1 ; i < A.size (); ++i){ dp[i] = max (A[i], dp[i - 1 ] + A[i]); ret = max (ret, dp[i]); } return ret; }
计算最大子段和所在区间 $O(n)$ 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 int dp[MAXN];pair<int , int > seq[MAXN]; pair<int, int> MCS(vector<int>& A){ dp[0 ] = A[0 ]; seq[0 ] = make_pair(0 , 0 ); for (int i = 1 ; i < A.size (); ++i){ dp[i] = max (A[i], dp[i - 1 ] + A[i]); if (A[i] < dp[i - 1 ] + A[i]) seq[i].first = seq[i - 1 ].first, seq[i].second = i; else seq[i].first = seq[i].second = i; } int ret = 0 ; for (int i = 0 ; i < A.size (); ++i) if (dp[ret] < dp[i]) ret = i; return seq[ret]; }
最大双子段和 打擂法 $O(n)$ 给定一个长度为n的整数序列,要求从中选出两个连续子序列,使得这两个连续子序列的序列和之和最大,最终只需输出最大和。一个连续子序列的和为该子序列中所有数之和。每个连续子序列的最小长度为1,并且两个连续子序列之间至少间隔一个数 。
预处理前i - 1
个数的最大和,后i + 1
个数的最大和。那么只需枚举i
即可求出最大双序列和。
1 2 3 4 5 6 7 8 9 10 11 12 ll pre[MAXN], succ[MAXN]; ll MDCS (vector <ll>& A) { int n = A.size (); pre[0 ] = A[0 ]; succ[n - 1 ] = A[n - 1 ]; _for(i, 1 , n) pre[i] = max (A[i], pre[i - 1 ] + A[i]); _for(i, 1 , n) pre[i] = max (pre[i-1 ], pre[i]); for (int i = n - 2 ; i >= 0 ; --i) succ[i] = max (A[i], succ[i+1 ] + A[i]); for (int i = n - 2 ; i >= 0 ; --i) succ[i] = max (succ[i+1 ], succ[i]); ll ret = pre[0 ] + succ[2 ]; _rep(i, 2 , n - 2 ) ret = max (ret, pre[i-1 ] + succ[i+1 ]); return ret; }
最大m子段和 $O(mn)$ 给定由n个整数(可能为负整数)组成的序列,以及一个正整数m,要求确定序列m个不相交子段,使得这m个子段的总和达到最大,求出最大和。
1 2 3 4 5 6 7 8 9 10 11 12 ll f[MAXN][5 ][2 ]; ll MCMS (vector <ll>& A, ll m) { ll n = A.size (); memset (f,-0x7f ,sizeof (f)); f[1 ][0 ][0 ] = 0 ; f[1 ][1 ][1 ] = A[0 ]; _rep(i, 2 , n) _rep(j, 0 , m){ f[i][j][0 ] = max (f[i-1 ][j][1 ], f[i-1 ][j][0 ]); if (j > 0 ) f[i][j][1 ] = max (f[i-1 ][j][1 ], f[i-1 ][j-1 ][0 ]) + A[i - 1 ]; } return max (f[n][m][0 ], f[n][m][1 ]); }
受限最大子段和 $O(n)$ 计算:子段长度不大于m 的最大子段和;子段长度不小于m 的最大子段和。
利用前缀和+单调队列+滑动窗口。动态维护窗口内的sum[i]
的最值。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 int sum[MAXN];int RMCS (vector <int >& A, int m) { int ret= A[0 ], n = A.size (); sum[0 ] = 0 ; sum[1 ] = A[0 ]; _rep(i, 2 , n) sum[i] = sum[i - 1 ] + A[i - 1 ]; deque <int > MQ; MQ.push_back(0 ); _rep(i, 1 , n){ if (i - MQ.front() >= m + 1 ) MQ.pop_front(); while (!MQ.empty() && sum[MQ.back()] > sum[i]) MQ.pop_back(); MQ.push_back(i); ret = max (ret, sum[i] - sum[MQ.front()]); } return ret; }
受限最小子段长度 $O(n)$ 给出了N个正整数(10 <N <100 000)的序列, 每个正整数均小于或等于10000,并给出了一个正整数S(S <100 000 000)。 编写程序以查找序列中连续元素的子序列的最小长度,其总和大于或等于S。
尺取法。
1 2 3 4 5 6 7 8 9 10 11 12 13 int RMCL (vector <int >& A, int S) { int ret = 1 , sum = A[0 ], i = 0 , j = 1 ; while (j < A.size ()){ while (sum < S && j < A.size ()){ sum += A[j++]; } while (sum >= S){ ret = max (ret, j - i); sum -= A[i++]; } } return ret; }
环状最大子段和 $O(2n)$
第一种方法:将序列倍长,就可以做一遍长度不大于n的受限最大子段和。
第二种方法:最大子段和 + 打擂法计算最大跨区子段和(预处理前缀、后缀和)。
第三种方法:最大子段和 + 最小子段和。【此处实现】
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 int dp[MAXN];int MCS (vector <int >& A) { int ret = dp[0 ] = A[0 ]; for (int i = 1 ; i < A.size (); ++i){ dp[i] = max (A[i], dp[i - 1 ] + A[i]); ret = max (ret, dp[i]); } return ret; } int CMCS (vector <int >& A) { int sum = accumulate(A.begin (), A.end (), 0 ); int ret = MCS(A); _for(i, 0 , A.size ()) A[i] = - A[i]; return max (ret, sum + MCS(A)); }
环状最大双子段和
最大双子段和 + 最小双子段和。(注意子段之间不相邻)
缺点:无法适用于单一段(++++----++++
、------
),此时会退化成环状最大子段和。需要特判。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 int pre[MAXN], succ[MAXN];int MDCS (vector <int >& A) { int n = A.size (); pre[0 ] = A[0 ]; succ[n - 1 ] = A[n - 1 ]; _for(i, 1 , n) pre[i] = max (A[i], pre[i - 1 ] + A[i]); _for(i, 1 , n) pre[i] = max (pre[i-1 ], pre[i]); for (int i = n - 2 ; i >= 0 ; --i) succ[i] = max (A[i], succ[i+1 ] + A[i]); for (int i = n - 2 ; i >= 0 ; --i) succ[i] = max (succ[i+1 ], succ[i]); int ret = pre[0 ] + succ[2 ]; _rep(i, 2 , n - 2 ) ret = max (ret, pre[i-1 ] + succ[i+1 ]); return ret; } int CMDCS (vector <int >& A) { int sum = accumulate(A.begin (), A.end (), 0 ); int ret = MDCS(A); _for(i, 0 , A.size ()) A[i] = - A[i]; return max (ret, sum + MDCS(A)); }
环状最大双子段和:【全能版】适用于单一段(++++----++++
、------
)
给出一段长度为 n 的环状序列 a,即认为 $a_1$ 和 $a_n$ 是相邻的,选出其中连续不重叠且非空 的两段使得这两段和最大。(如果形如+++++++
,取全局SUM值)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 #include <cstdio> #include <algorithm> #include <cstring> using namespace std ;const int INF=0x3f3f3f3f ;const int N=2e5 +5 ;int n,f[N],g[N],a[N],sum=0 ,tot=0 ;int query () { int res=-INF; for (int i=1 ;i<=n;i++)f[i]=max(f[i-1 ],0 )+a[i]; for (int i=n;i>0 ;i--)g[i]=max(g[i+1 ],0 )+a[i]; for (int i=1 ;i<=n;i++)f[i]=max(f[i-1 ],f[i]); for (int i=n;i>0 ;i--)g[i]=max(g[i+1 ],g[i]); for (int i=1 ;i<n;i++)res=max(res,f[i]+g[i+1 ]); return res; } int main () { scanf ("%d" ,&n); memset (f,~0x3f ,sizeof (f));memset (g,~0x3f ,sizeof (g)); for (int i=1 ;i<=n;i++)scanf ("%d" ,&a[i]),sum+=a[i],tot+=a[i]>0 ; int t1=query(); if (tot==1 ){ printf ("%d" ,t1); }else { for (int i=1 ;i<=n;i++)a[i]=-a[i]; int t2=sum+query(); if (!t2)t2=-INF; printf ("%d" ,max(t1,t2)); } return 0 ; }
环状最大合并 $O(n^3)$ 题目大意:在一个环上有 $n$ 个数 $a_1,a_2,…,a_n$,进行 $n-1$ 次合并操作,每次操作将相邻的两堆合并成一堆,能获得新的一堆中的石子数量的和的得分。你需要最大化你的得分。
记 f[i][j]
为区间 [i, j]
合并为一堆时的最大可能得分。
环的处理:倍增原序列。
1 2 3 4 5 6 for (len = 1 ; len <= n; len++) for (i = 1 ; i <= 2 * n - 1 ; i++) { int j = len + i - 1 ; for (k = i; k < j && k <= 2 * n - 1 ; k++) f[i][j] = max (f[i][j], f[i][k] + f[k + 1 ][j] + sum[j] - sum[i - 1 ]); }
最大子矩阵和 $O(n^3)$ 枚举行i
和j
,
i == j
:等价于求第i
行元素的最大子段和。
i != j
:将i
到j
行压缩为1行元素,从而转化为求最大子段和。(可用前缀和预处理优化压缩效率)
最长上升子序列LIS $O(n^2) + O(n)$
F[i]
表示[0,i]
子序列的LIS长度或其他标量。
1 2 3 4 5 6 7 8 int F[MAXN];int LIS (vector <int > a) { for (int i = 0 ; i < a.size (); ++i){ F[i] = 1 ; for (int j = 0 ; j < i; ++j) if (a[j] < a[i]) F[i] = max (F[i], F[j] + 1 ); } }
LIS 二分DP $O(nlogn) + O(n)$
dp
数组表示长度为i+1
的LIS结尾元素的最小值。
1 2 3 4 5 6 7 8 int dp[MAXN];void LIS (vector <int > a) { fill (dp, dp+n, INF); for (int i = 0 ; i < n; i++) { *lower_bound(dp, dp + n, a[i]) = a[i]; } printf ("%d\n" , lower_bound(dp, dp + n, INF) - dp); }
最长公共子序列LCS $O(n^2) + O(n)$ 1 2 3 4 5 6 7 int dp[2 ][MAXN]; int LCS (string A, string B) { for (int i = 0 ; i < A.size (); ++i) for (int j = 0 ; j < B.size (); ++j) dp[(i+1 )&1 ][j+1 ] = (A[i]==B[j]) ? (dp[i&1 ][j]+1 ) : max (dp[(i+1 )&1 ][j], dp[i&1 ][j+1 ]); return dp[A.size ()&1 ][B.size ()]; }
LCS方案计数
如果需要打印路径,则在计算时要记录局部最优时的前驱。
LCS最优路径只须(用指针)记录贡献者,进行回溯即可。
LCS计数取三种状态中被继承的若干状态,若左、上共同贡献则进行容斥 。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 int dp[2 ][MAXN], sum[2 ][MAXN];int LCS (string A, string B) { fill (sum[0 ], sum[0 ] + MAXN + 1 , 1 ); for (int i = 0 ; i < A.size (); ++i){ for (int j = 0 ; j < B.size (); ++j){ dp[(i+1 )&1 ][j+1 ] = (A[i] == B[j]) ? (dp[i&1 ][j] + 1 ) : max (dp[(i+1 )&1 ][j], dp[i&1 ][j+1 ]); int a = dp[i&1 ][j], b = dp[(i+1 )&1 ][j], c = dp[i&1 ][j+1 ] , Max = dp[(i+1 )&1 ][j+1 ], & Now = sum[(i+1 )&1 ][j+1 ] = 0 ; if (A[i] == B[j] && a + 1 == Max) Now += sum[i&1 ][j]; if (b == Max) Now += sum[(i+1 )&1 ][j]; if (c == Max) Now += sum[i&1 ][j+1 ]; if (a == Max) Now -= sum[i&1 ][j]; } } return sum[(A.size ()) &1 ][B.size ()]; }
最优三角剖分 $O(n^3)$ d(i,j)
表示子多边形i,I+1, ..., j-1, j
的最优解。
边界$d(i, i+1) = 0$,原问题的解为$d(0,n-1)$。
树形DP 树形DP的主要实现方式是DFS,在DFS中DP。
最大独立集 树上节点的一个最大子集,其中任何两个节点均不相邻,称为最大独立集。
状态转移方程:(gs(i)
是i
的孙节点集合,s(i)
是i
的儿子节点集合)
节点i
的决策:选和不选。如果选,求出 1 + 孙节点的d(j)
即可;若不选,只需求出儿子节点的d(j)
。
动态最大权独立集
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 #include <cstdio> #include <algorithm> using namespace std ;const int N=1e5 +10 ;int n;int m;int v[2 *N];int x[2 *N];int ct;int al[N];int siz[N];int h[N];int we[N];inline void add (int u,int V) {v[++ct]=V;x[ct]=al[u];al[u]=ct;}inline int dfs1 (int u) { siz[u]=1 ;int mx=0 ; for (int i=al[u];i;i=x[i]) if (siz[v[i]]==0 ){siz[u]+=dfs1(v[i]);if (mx<siz[v[i]])mx=siz[v[i]],h[u]=v[i];} return siz[u]; } struct mar //矩阵类 { int mp[2 ][2 ]; mar(){mp[0 ][0 ]=mp[0 ][1 ]=mp[1 ][0 ]=mp[1 ][1 ]=-0x3f3f3f3f ;} mar(int x){mp[0 ][0 ]=mp[1 ][1 ]=0 ;mp[1 ][0 ]=mp[0 ][1 ]=-0x3f3f3f3f ;} friend mar operator *(mar a,mar b) { mar c;for (int i=0 ;i<2 ;++i) for (int k=0 ;k<2 ;++k) for (int j=0 ;j<2 ;++j)c.mp[i][j]=max (c.mp[i][j],a.mp[i][k]+b.mp[k][j]); return c; } inline int gmx () {return max (max (mp[0 ][0 ],mp[0 ][1 ]),max (mp[1 ][0 ],mp[1 ][1 ]));} inline int * operator [](const int & x){return mp[x];} }; struct bst { int s[N][2 ];int fa[N];int st[N];int tp;int lsiz[N];bool book[N];int root; mar mul[N];mar w[N];bst(){w[0 ]=mul[0 ]=mar(1 );} inline void ud (const int & x) {mul[x]=mul[s[x][0 ]]*w[x]*mul[s[x][1 ]];} inline void gtw (const int & x,const int & v) {w[x][1 ][0 ]+=mul[v].gmx();w[x][0 ][0 ]=w[x][1 ][0 ];w[x][0 ][1 ]+=max (mul[v][0 ][0 ],mul[v][1 ][0 ]);fa[v]=x;} inline void ih () {for (int i=1 ;i<=n;i++)w[i][0 ][1 ]=we[i],w[i][0 ][0 ]=w[i][1 ][0 ]=0 ;} inline bool isr (const int & p) {return (s[fa[p]][1 ]!=p)&&(s[fa[p]][0 ]!=p);} inline int sbuild (const int & l,const int & r) { if (l>r)return 0 ;int tot=0 ;for (int i=l;i<=r;i++)tot+=lsiz[st[i]]; for (int i=l,ns=lsiz[st[i]];i<=r;i++,ns+=lsiz[st[i]]) if (2 *ns>=tot) { int rs=sbuild(l,i-1 );int ls=sbuild(i+1 ,r);s[st[i]][0 ]=ls;s[st[i]][1 ]=rs; fa[ls]=st[i];fa[rs]=st[i];ud(st[i]);return st[i]; } } inline int build (int p) { for (int t=p;t;t=h[t])book[t]=true ; for (int t=p;t;t=h[t]) for (int i=al[t];i;i=x[i])if (!book[v[i]])gtw(t,build(v[i])); tp=0 ;for (int t=p;t;t=h[t])st[++tp]=t; for (int t=p;t;t=h[t])lsiz[t]=siz[t]-siz[h[t]];return sbuild(1 ,tp); } inline void modify (int p,int W) { w[p][0 ][1 ]+=W-we[p];we[p]=W; for (int t=p;t;t=fa[t]) if (isr(t)&&fa[t]) { w[fa[t]][0 ][0 ]-=mul[t].gmx();w[fa[t]][1 ][0 ]=w[fa[t]][0 ][0 ]; w[fa[t]][0 ][1 ]-=max (mul[t][0 ][0 ],mul[t][1 ][0 ]);ud(t); w[fa[t]][0 ][0 ]+=mul[t].gmx();w[fa[t]][1 ][0 ]=w[fa[t]][0 ][0 ]; w[fa[t]][0 ][1 ]+=max (mul[t][0 ][0 ],mul[t][1 ][0 ]); }else ud(t); } }bst; int main () { scanf ("%d%d" ,&n,&m); for (int i=1 ;i<=n;i++)scanf ("%d" ,&we[i]); for (int i=1 ,u,v;i<n;i++){scanf ("%d%d" ,&u,&v);add(u,v);add(v,u);} dfs1(1 );bst.ih();bst.root=bst.build(1 ); for (int i=1 ,p,w;i<=m;i++) {scanf ("%d%d" ,&p,&w);bst.modify(p,w);printf ("%d\n" ,bst.mul[bst.root].gmx());} return 0 ; }
缩点 Tarjan算法 给定一个 n 个点 m 条边有向图,每个点有一个权值,求一条路径,使路径经过的点权值之和最大 。你只需要求出这个权值和。允许多次经过一条边或者一个点,但是,重复经过的点,权值只计算一次。
若将环缩成点,则原图变为DAG。就可以拓扑排序然后DP辣。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 #include <iostream> #include <cstdio> #include <cstring> #include <queue> #include <algorithm> using namespace std ;#define N 100010 struct edge { int u, v, nxt; }e[N]; queue <int > q;int head[N], cnt, n, m, x, y, vis[N], f[N], c[N], s[N], r[N], ans = 1 , W[N];inline void add (int u, int v) { e[++cnt].u = u; e[cnt].v = v; e[cnt].nxt = head[u]; head[u] = cnt; } int anc (int x) { if (x == f[x]) return x; return f[x] = anc(f[x]); } void dfs (int u) { for (int i = head[u]; i; i = e[i].nxt){ int v = e[i].v; if (!vis[v]){ vis[v] = vis[u] + 1 ; dfs(v); } int fu = anc(u), fv = anc(v); if (vis[fv] > 0 ){ if (vis[fu] < vis[fv]) f[fv] = fu; else f[fu] = fv; } } vis[u] = -1 ; return ; } int main () { scanf ("%d%d" , &n, &m); for (int i = 1 ; i <= n; ++i) f[i] = i, scanf ("%d" , W + i); for (int i = 1 ; i <= m; ++i){ scanf ("%d%d" , &x, &y); add(x, y); } for (int i = 1 ; i <= n; ++i){ if (!vis[i]){ vis[i] = 1 ; dfs(i); } } for (int i = 1 ; i <= n; ++i) c[anc(i)] += W[i]; memset (head, 0 , sizeof (head)); cnt = 0 ; for (int i = 1 ; i <= m; ++i){ x = f[e[i].u], y = f[e[i].v]; if (x != y) add(x, y), r[y]++; } for (int i = 1 ; i <= n; ++i){ if (i == f[i] && !r[i]) q.push(i), s[i] = c[i]; } while (!q.empty()){ int u = q.front(); q.pop(); for (int i = head[u]; i; i = e[i].nxt){ int v = e[i].v; s[v] = max (s[v], s[u] + c[v]); if (--r[v] == 0 ) q.push(v); } ans = max (ans, s[u]); } printf ("%d" , ans); return 0 ; }
最小树形图 给定包含 n 个结点, m 条有向边的一个图。试求一棵以结点 r 为根的最小树形图,并输出最小树形图每条边的权值之和,如果没有以 r 为根的最小树形图,输出 −1。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 #include <iostream> #include <cstdio> #include <cstring> #include <cstdlib> #define il inline #define ri register #define Size 100050 using namespace std ;int fa[Size],cnt,is[Size];il int find (int ) ;il void read(int&),Union(int,int); struct leftist { struct point { int l,r,d,v,t,to; }a[Size]={{0 ,0 ,-1 ,0 ,0 ,0 }}; int r[Size]; il void merge (int &x,int &y) { if (!x||!y){x^=y;return ;} if (a[x].v>a[y].v)x^=y^=x^=y; a[y].t-=a[x].t,a[y].v-=a[x].t; merge(a[x].r,y); if (a[a[x].l].d<a[a[x].r].d) a[x].l^=a[x].r^=a[x].l^=a[x].r; a[x].d=a[a[x].r].d+1 ; } il void spread (int &p) { a[a[p].l].t+=a[p].t,a[a[p].r].t+=a[p].t; a[a[p].l].v+=a[p].t,a[a[p].r].v+=a[p].t; a[p].t=0 ; } il void pop (int &x) { spread(x),merge(a[x].l,a[x].r),x=a[x].l; } il point *top (int &x) { while (r[x]&&!(find (a[r[x]].to)^x))pop(r[x]); if (!r[x])puts ("-1" ),exit (0 ); a[r[x]].to=find (a[r[x]].to); return &a[r[x]]; } }L; int pre[Size];int main () { int n,m,r,ans(0 );leftist::point *temp; read (n),read (m),read (r),cnt=n,is[r]=r; for (int i(1 ),u,v,w;i<=m;++i) read (u),read (v),read (w), L.a[i]={0 ,0 ,0 ,w,0 ,u}, L.merge(L.r[v],u=i); for (int i(1 );i<=n<<1 ;++i)fa[i]=i; for (int i(1 ),j(i);i<=n;j=++i) while (!is[j]){ while (!is[j]) is[j]=i,j=(temp=L.top(j))->to, ans+=temp->v;if (is[j]^i)break ; while (~is[j]) is[j]=-1 ,j=pre[j]=(temp=L.top(j))->to, temp->t-=temp->v,temp->v=0 ;++cnt; while (is[j]^i)is[j]=i,Union(j,cnt),j=pre[j]; j=cnt; }return printf ("%d" ,ans),0 ; } il void Union (int u,int v) { if ((u=find (u))^(v=find (v))) L.merge(L.r[v],L.r[u]),fa[u]=v; } il int find (int x) { return x^fa[x]?fa[x]=find (fa[x]):x; } il void read (int &x) { x^=x;ri char c;while (c=getchar(),c<'0' ||c>'9' ); while (c>='0' &&c<='9' )x=(x<<1 )+(x<<3 )+(c^48 ),c=getchar(); }
树形DP例题 上司舞会
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 #include <iostream> #include <cstdio> using namespace std ;int const maxn=10011 ;int f[maxn][5 ],fa[maxn],num[maxn],son[maxn][maxn];bool kid[maxn];int n,v,x,y;void dfs (int u) { if (!son[u][0 ]){ f[u][1 ]=num[u]; f[u][0 ]=0 ; return ; } else { for (int i=1 ;i<=son[u][0 ];i++){ dfs(son[u][i]); f[u][0 ]+=max (f[son[u][i]][0 ],f[son[u][i]][1 ]); f[u][1 ]+=f[son[u][i]][0 ]; } f[u][1 ]+=num[u]; } } int main () { cin >>n; for (int i=1 ;i<=n;i++){ cin >>v; num[i]=v; } for (int i=1 ;i<n;i++){ cin >>x>>y; son[y][0 ]++; fa[x]++; son[y][son[y][0 ]]=x; } cin >>x>>y; int root; for (int i=1 ;i<=n;i++) if (!fa[i]){ root=i; break ; } dfs(root); cout <<max (f[root][0 ],f[root][1 ]); return 0 ; }
最大子树和 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 #include <bits/stdc++.h> using namespace std ;struct node { int next,to; }tree[100000 ]; int n,a[100000 ],head[100000 ],cnt,f[100000 ],ans; void add (int x,int y) { tree[++cnt].next=head[x]; tree[cnt].to=y; head[x]=cnt; } void dfs (int u,int fa) { f[u]=a[u]; for (int i=head[u];i;i=tree[i].next) { int v=tree[i].to; if (v!=fa) { dfs(v,u); f[u]+=max(0 ,f[v]); } } ans=max(ans,f[u]); } int main () { scanf ("%d" ,&n); for (int i=1 ;i<=n;i++) scanf ("%d" ,&a[i]); for (int i=1 ,x,y;i<n;i++) { scanf ("%d%d" ,&x,&y); add(x,y);add(y,x); } dfs(1 ,0 ); printf ("%d" ,ans); }
状压DP 常见于将一个集合的状态压缩为一个数。
数位DP 求出在给定区间 $[A,B]$ 内,符合条件 $f(i)$ 的数 $i$ 的个数。条件 $f(i)$ 一般与数的大小无关,而与数的组成有关。
数位DP = DFS + 记忆化。由于数是按位dp,数的大小对复杂度的影响很小。
数位DP的简单状态方程:
定义$f(i,lim)$表示当前将要考虑的是从高到低的第 $i$ 位。$\text{max}\;x$ 表示就是当前能取到的最高位。因为如果 $lim=1$,那么你在这一位上取的值一定不能大于求解的数字上该位的值,否则则没有限制。
例:(给定两个正整数 $a$ 和 $b$,求在 $[a,b]$ 中的所有整数中,每个数码(digit)各出现了多少次。)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 #include <bits/stdc++.h> using namespace std ;long long a,b;long long ten[20 ],f[20 ];long long cnta[20 ],cntb[20 ];void solve (long long x,long long *cnt) { long long num[20 ]={0 }; int len=0 ; while (x){ num[++len]=x%10 ; x=x/10 ; } for (int i=len; i>=1 ; i--){ for (int j=0 ; j<=9 ; j++) cnt[j] += f[i-1 ] * num[i]; for (int j=0 ; j<num[i]; j++) cnt[j] += ten[i-1 ]; long long num2=0 ; for (int j= i-1 ; j>=1 ; j--){ num2 = num2 * 10 + num[j]; } cnt[num[i]]+=num2+1 ; cnt[0 ]-= ten[i-1 ]; } } int main () { scanf ("%lld %lld" ,&a,&b); ten[0 ]=1 ; for (int i=1 ;i<=15 ;i++) { f[i]=f[i-1 ]*10 +ten[i-1 ]; ten[i]=10 *ten[i-1 ]; } solve(a-1 ,cnta); solve(b,cntb); for (int i=0 ;i<=9 ;i++) printf ("%lld " ,cntb[i]-cnta[i]); }
数学 快速幂 $O(logn)$
减治法。对于矩阵幂只需重载乘法 即可。见矩阵部分。
Python自带快速幂pow(a, n)
。
1 2 3 4 5 6 int qPower (int a, int n) { if (n == 0 ) return 1 ; if (n == 1 ) return a; int base = qPower(a, n/2 ); base *= base; return (n % 2 ) ? a*base : base; }
迭代快速幂 $O(logn)+O(1)$ 1 2 3 4 5 6 7 8 9 int qPower (int a, int n) { int ans = 1 , base = a; while (n){ if (n & 1 ) ans *= base; base *= base; n >>= 1 ; } return ans; }
取模快速幂 1 2 3 4 5 6 7 8 9 10 int qPower (int a, int n, int mod) { int ans = 1 , base = a; while (n){ if (n & 1 ) ans *= base; base *= base; ans %= mod; base %= mod; n >>= 1 ; } return ans; }
费马小定理求模幂 $O(logP)$ 计算模意义下的大规模的指数。p必须为质数 。
因为 $a^{p-1}\equiv 1(\text{mod}\,\,p)$,所以 $a^{b}\equiv a^{b\;\text{mod}\;(p-1)}\;(\text{mod}\,\,p)$。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 int qPower (int a, int n, int mod) { int ans = 1 , base = a; while (n){ if (n & 1 ) ans *= base; base *= base; ans %= mod; base %= mod; n >>= 1 ; } return ans; } int qBigPower (int a, int n, int p) { return qPower(a, n % (p - 1 ), p); }
扩展欧拉定理求模幂 $O(log\;\varphi(P))$ 计算模意义下的超大规模的指数($b\leq 10^{10000000}$)。p可不为质数。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 #include <bits/stdc++.h> #define ll long long using namespace std ;ll a,m,b; inline ll read (ll m) { register ll x=0 ,f=0 ;char ch=getchar(); while (!isdigit (ch)) ch=getchar(); while (isdigit (ch)){ x=x*10 +ch-'0' ; if (x>=m) f=1 ; x%=m;ch=getchar(); } return x+(f==1 ?m:0 ); } ll phi (ll n) { ll ans=n,m=sqrt (n); for (ll i=2 ;i<=m;i++){ if (n%i==0 ){ ans=ans/i*(i-1 ); while (n%i==0 ) n/=i; } } if (n>1 ) ans=ans/n*(n-1 ); return ans; } ll qBigPower (ll a,ll b,ll p) { ll ret=1 ; for (;b;b>>=1 ,a=a*a%p) if (b&1 ) ret=ret*a%p; return ret; } int main () { scanf ("%lld%lld" ,&a,&m); b=read (phi(m)); printf ("%lld\n" ,qBigPower(a,b,m)); return 0 ; }
矩阵快速幂/矩阵加速 $O(logn)$
需要使用矩阵类。
当n的规模很大时,矩阵快速幂可以用来加快计算固定递推公式的数列 。
所有的有限常递推公式,都可以使用矩阵快速幂。写成矩阵形式即可。
1 2 3 4 5 6 int fib (int n) { Matrix F (2 , 1 ) ; F.a[0 ][0 ] = 1 ; Matrix A (2 , 2 ) ; A.a[0 ][0 ] = A.a[1 ][0 ] = A.a[0 ][1 ] = 1 ; F = qPower(A, n) * F; return F.a[1 ][0 ]; }
光速幂 $O(T)$ 对$a$不变的$T$次不同$t$的查询,快速幂需要$O(Tlogn)$。而光速幂只需$O(T)$。取$n=max\{t_i\}$,
首先分别对$a^\sqrt{n}$和$a$的$0$到$\sqrt{n}$次幂进行预处理。这样之后每次计算$a^t$都可以在$O(1)$时间内得到。
快速乘 有时long long
取模乘法中途会溢出,因此需要快速乘处理,
1 2 3 4 5 6 7 8 9 ll qMul (ll n, ll k, ll mod) { ll ans = 0 ; while (k){ if (k & 1 ) ans = (ans + n) % mod; k >>= 1 ; n = (n + n) % mod; } return ans; }
公约数/公倍数 $O(logn)$ 最大公约数 gcd
(欧几里德算法) 1 2 3 inline int gcd (int a, int b) { return b == 0 ? a : gcd(b, a%b); }
STL:__gcd(a,b)
。(部分OJ不支持)
1 2 3 4 5 6 7 8 int kgcd (int a, int b) { if (a == 0 ) return b; if (b == 0 ) return a; if (!(a & 1 ) && !(b & 1 )) return kgcd(a>>1 , b>>1 ) << 1 ; else if (!(b & 1 )) return kgcd(a, b>>1 ); else if (!(a & 1 )) return kgcd(a>>1 , b); else return kgcd(abs (a - b), min (a, b)); }
最小公倍数 lcm
1 2 3 int lcm (int a, int b) { return a*b / gcd(a, b); }
扩展欧几里得算法 exgcd $ax+by=gcd(a,b)$ 扩展欧几里得算法用于求解形如$ax+by=gcd(a,b)$的不定方程特解。
裴蜀定理/贝祖定理:
设$a,b$是不全为零的整数,则存在整数$x,y$,使得$ax+by=\text{gcd}(a,b)$。
1 2 3 4 5 void exgcd (int a, int b, int & x, int & y) { if (b == 0 ) { x = 1 , y = 0 ; return ;} exgcd(b, a % b, y, x); y -= a / b * x; }
$b=0$时,$\left\{\begin{array}{l} x=1 \\ y=0 \end{array}\right.$是方程的一组特解。否则递归求解$x’$和$y’$,并且有$\left\{\begin{array}{l} x=y^{\prime} \\ y=x^{\prime}-a / b * y^{\prime} \end{array}\right.$。
求解 $ax+by=c$
线性同余方程 $ax\equiv b(\text{mod}\,\,n)$ 上述方程等价于:$ax-b = n\times c$。即等价于:$ax+nc = b$。即等价于$ax+by=c$。可用exgcd
方法求解。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 int ex_gcd (int a, int b, int & x, int & y) { if (b == 0 ) { x = 1 ; y = 0 ; return a; } int d = ex_gcd(b, a % b, x, y); int temp = x; x = y; y = temp - a / b * y; return d; } bool liEu (int a, int b, int n, int & x, int & y) { int d = ex_gcd(a, b, x, y); if (n % d != 0 ) return 0 ; int k = n / d; x *= k; y *= k; return 1 ; }
扩展裴蜀定理
任何可行解必然是 $gcd(X_1,X_2,…,X_n)$ 的倍数,故 $S = gcd(X_1,X_2,…,X_n)$ 。
1 2 3 4 5 int exPeiShu (int a[], int n) { int ret = 0 ; _for(i, 0 , n) ret = gcd(ret, abs (a[i])); return ret; }
二次剩余/模意义下开根 求解 $x^{2} \equiv N(\bmod p)$
保证p是奇素数 。
输入格式
第一行一个整数T表示数据组数。
接下来T行,每行一个N一个p。
输出格式
输出共T行。
对于每一行输出:
若有解,则按 $\bmod ~p$ 后递增的顺序输出在 $\bmod~ p$ 意义下的全部解.
若两解相同,只输出其中一个;
若无解,则输出Hola!
;
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 #include <bits/stdc++.h> using namespace std ;typedef long long ll;ll w; struct num { ll x,y; }; num mul (num a,num b,ll p) { num ans={0 ,0 }; ans.x=((a.x*b.x%p+a.y*b.y%p*w%p)%p+p)%p; ans.y=((a.x*b.y%p+a.y*b.x%p)%p+p)%p; return ans; } ll powwR (ll a,ll b,ll p) { ll ans=1 ; while (b){ if (b&1 )ans=1l l*ans%p*a%p; a=a%p*a%p; b>>=1 ; } return ans%p; } ll powwi (num a,ll b,ll p) { num ans={1 ,0 }; while (b){ if (b&1 )ans=mul(ans,a,p); a=mul(a,a,p); b>>=1 ; } return ans.x%p; } ll solve (ll n,ll p) { n%=p; if (p==2 )return n; if (powwR(n,(p-1 )/2 ,p)==p-1 )return -1 ; ll a; while (1 ) { a=rand()%p; w=((a*a%p-n)%p+p)%p; if (powwR(w,(p-1 )/2 ,p)==p-1 )break ; } num x={a,1 }; return powwi(x,(p+1 )/2 ,p); } int main () { srand(time(0 )); int t; scanf ("%d" ,&t); while (t--) { ll n,p; scanf ("%lld%lld" ,&n,&p); if (!n){ printf ("0\n" );continue ; } ll ans1=solve(n,p),ans2; if (ans1==-1 )printf ("Hola!\n" ); else { ans2=p-ans1; if (ans1>ans2)swap(ans1,ans2); if (ans1==ans2)printf ("%lld\n" ,ans1); else printf ("%lld %lld\n" ,ans1,ans2); } } }
BSGS离散对数 $a^x \equiv b \pmod p$ $O(\sqrt n)$ 北上广深大步小步算法(Bady Step Giant Step)。
给定一个质数 $p$,以及一个整数 $b$,一个整数 $n$,现在要求你计算一个最小的 $l$,满足$b^l \equiv n \pmod p$。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 #include <bits/stdc++.h> #define ll long long using namespace std ;ll p,b,n,m; map <ll,int >vis; ll quick_pow (ll a,ll b,ll p) { ll res=1 ; for (;b;b>>=1 ){ if (b&1 )res=res*a%p; a=a*a%p; } return res; } int main () { scanf ("%lld%lld%lld" ,&p,&b,&n); m=ceil (sqrt (p)); for (ll i=0 ,t=n;i<=m;i++,t=t*b%p)vis[t]=i; for (ll i=1 ,tt=quick_pow(b,m,p),t=tt;i<=m;i++,t=t*tt%p) if (vis.count(t))printf ("%lld\n" ,i*m-vis[t]),exit (0 ); puts ("no solution" ); return 0 ; }
扩展BSGS算法 $a^x \equiv b \pmod p$
$p$ 可以不为素数。计算 $x$ 的最小自然数解。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 #include <bits/stdc++.h> #define ll long long using namespace std ; unordered_map <ll, int > H ;int N, M, P, ans ; inline ll gcd (ll a, ll b) { if (!b) return a ; return gcd(b, a % b) ; } inline ll expow (ll a, ll b, ll mod) { ll res = 1 ; while (b) res = ((b & 1 )?res * a % mod : res), a = a * a % mod, b >>= 1 ; return res ; } inline ll exgcd (ll &x, ll &y, ll a, ll b) { if (!b){ x = 1 , y = 0 ; return a ; } ll t = exgcd(y, x, b, a % b) ; y -= x * (a / b) ; return t ; } inline ll BSGS (ll a, ll b, ll mod, ll qaq) { H.clear () ; ll Q, p = ceil (sqrt (mod)), x, y ; exgcd(x, y, qaq, mod), b = (b * x % mod + mod) % mod, Q = expow(a, p, mod), exgcd(x, y, Q, mod), Q = (x % mod + mod) % mod ; for (ll i = 1 , j = 0 ; j <= p ; ++ j, i = i * a % mod) if (!H.count(i)) H[i] = j ; for (ll i = b, j = 0 ; j <= p ; ++ j, i = i * Q % mod) if (H[i]) return j * p + H[i] ; return -1 ; } inline ll exBSGS () { ll qaq = 1 ; ll k = 0 , qwq = 1 ; if (M == 1 ) return 0 ; while ((qwq = gcd(N, P)) > 1 ){ if (M % qwq) return -1 ; ++ k, M /= qwq, P /= qwq, qaq = qaq * (N / qwq) % P ; if (qaq == M) return k ; } return (qwq = BSGS(N, M, P, qaq)) == -1 ? -1 : qwq + k ; } int main () { while (cin >> N){ scanf ("%d%d" , &P, &M); if (!N && !M && !P) return 0 ; N %= P, M %= P, ans = exBSGS() ; if (ans < 0 ) puts ("No Solution" ) ; else cout << ans << '\n' ; } }
类欧几里得算法 $O(logn)$ 这个算法用于求一条直线下整点个数,我们定义:
1 2 3 4 5 6 7 ll f (ll a, ll b, ll c, ll n) { if (!a) return b / c * (n + 1 ); if (a >= c || b >= c) return f(a % c, b % c, c, n) + (a / c) * n * (n + 1 ) / 2 + (b / c) * (n + 1 ); ll m = (a * n + b) / c; return n * m - f(c, c - b - 1 , a, m - 1 ); }
中国剩余定理 CRT 中国剩余定理 (Chinese Remainder Theorem, CRT) 可求解如下形式的一元线性同余方程组:
其中 $m_1,m_2,…,m_k$ 两两互质。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 void exgcd (ll a, ll b, ll& x, ll& y) { if (b == 0 ) { x = 1 , y = 0 ; return ;} exgcd(b, a % b, y, x); y -= a / b * x; } int CRT (int a[], int m[], int n) { int M = 1 ; int ans = 0 ; for (int i = 1 ; i <= n; i++) { M *= m[i]; } for (int i = 1 ; i <= n; i++) { int x, y; int Mi = M / m[i]; exgcd(Mi, m[i], x, y); ans = (ans + Mi * x * a[i]) % M; } if (ans < 0 ) ans += M; return ans; }
某些计数问题或数论问题出于加长代码、增加难度、或者是一些其他不可告人的原因,给出的模数: 不是质数 !
但是对其质因数分解会发现它没有平方因子,也就是该模数是由一些不重复的质数相乘得到。
那么我们可以分别对这些模数进行计算,最后用 CRT 合并答案。
扩展中国剩余定理 当 $m_1,m_2,…,m_k$ 不满足 两两互质时使用。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 typedef long long ll;const int maxn=1e5 +5 ;int n; ll exgcd (ll a,ll b,ll &x,ll &y) { if (!b){x=1 ,y=0 ;return a;} ll re=exgcd(b,a%b,x,y),tmp=x; x=y,y=tmp-(a/b)*y; return re; } ll m[maxn],a[maxn]; ll exCRT () { ll M=m[1 ],A=a[1 ],t,d,x,y;int i; for (i=2 ;i<=n;i++){ d=exgcd(M,m[i],x,y); if ((a[i]-A)%d)return -1 ; x*=(a[i]-A)/d,t=m[i]/d,x=(x%t+t)%t; A=M*x+A,M=M/d*m[i],A%=M; } A=(A%M+M)%M; return A; } int main () { int i,j; while (scanf ("%d" ,&n)!=EOF){ for (i=1 ;i<=n;i++)scanf ("%lld%lld" ,&m[i],&a[i]); printf ("%lld\n" ,exCRT()); } return 0 ; }
快速乘-防long long
溢出版:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 typedef long long ll;const int maxn=1e5 +5 ;ll qMul (ll n, ll k, ll mod) { ll ans = 0 ; while (k){ if (k & 1 ) ans = (ans + n) % mod; k >>= 1 ; n = (n + n) % mod; } return ans; } int n; ll exgcd (ll a,ll b,ll &x,ll &y) { if (!b){x=1 ,y=0 ;return a;} ll re=exgcd(b,a%b,x,y),tmp=x; x=y,y=tmp-(a/b)*y; return re; } ll m[maxn],a[maxn]; ll exCRT () { ll M=m[1 ], ans=a[1 ], x, y; for (ll i=2 ;i<=n;++i){ ll c=((a[i]-ans)%m[i]+m[i])%m[i]; ll g=exgcd(M,m[i],x,y); ll p=m[i]/g; x=qMul(x,c/g,p); ans+=x*M; M*=p; ans=(ans%M+M)%M; } return ans; } int main () { int i,j; while (scanf ("%d" ,&n)!=EOF){ for (i=1 ;i<=n;i++)scanf ("%lld%lld" ,&m[i],&a[i]); printf ("%lld\n" ,exCRT()); } return 0 ; }
乘法逆元 扩展欧几里得求逆元 $ax\equiv 1(\text{mod}\,\,b)$ 如果一个线性同余方程$ax\equiv 1(\text{mod}\,\,b)$,则$x$称为$a\;\text{mod}\;b$的逆元,记为$a^{-1}$。
当$gcd(a, b) = 1$时,该方程有唯一解;否则无解。
1 2 3 4 5 6 7 8 9 10 void exgcd (ll a, ll b, ll& x, ll& y) { if (b == 0 ) { x = 1 , y = 0 ; return ;} exgcd(b, a % b, y, x); y -= a / b * x; } ll inverse (ll a, ll b) { ll x, y; exgcd(a, b, x, y); return (x % b + b) % b; }
费马小定理 $a^{m-1}\equiv 1(\text{mod}\,\,m)$ $O(logP)$ 设$m$是素数 ,$a$是任意整数且$a\not\equiv 0(\text{mod}\,\,m)$,则$a^{m-1}=a\times a^{m-2}\equiv 1(\text{mod}\,\,m)$。
故,$a^{m-2}\;\%\;m$ 就是 $a$ 模 $m$ 的逆元。
1 2 3 4 5 6 7 8 9 10 ll fpm (ll x, ll power, ll mod) { x %= mod; ll ans = 1 ; for (; power; power >>= 1 , (x *= x) %= mod) if (power & 1 ) (ans *= x) %= mod; return ans; } ll inverse (ll a, ll p) { return fpm(a, p - 2 , p); }
线性求逆元 $O(n)$ 用于求一连串数字对于一个mod p的逆元。
比如:给定 n,p, 求 1∼n 中所有整数在模 p 意义下的乘法逆元。
1 2 3 4 5 6 int inv[MAXN];void getinv (int p) { inv[1 ] = 1 ; for (int i = 2 ; i < p; ++ i) inv[i] = (p - p / i) * inv[p % i] % p; }
阶乘逆元 $O(n)$ 根据:$\cfrac{1}{(n+1)!}\times (n+1)=\cfrac{1}{n!}$。可以倒推求出 $1!,2!,…,n!$ 阶乘逆元。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 ll fpm (ll x, ll power, ll mod) { x %= mod; ll ans = 1 ; for (; power; power >>= 1 , (x *= x) %= mod) if (power & 1 ) (ans *= x) %= mod; return ans; } ll inv[MAXN]; void getFracInv (ll n, ll p) { ll frac = 1 ; for (ll i = 1 ; i <= n; ++i) frac = frac * i % p; inv[n] = fpm(frac, p - 2 , p); for (ll i = n - 1 ; i >= 1 ; --i) inv[i] = ( inv[i+1 ] * (i+1 ) ) % p; }
素数筛法 素数验证 $O(\sqrt N)$ 1 2 3 4 5 bool is_prime (int number) { for (int i = 2 ; i*i <= number; i+=1 ) if (!(number % i)) return 0 ; return 1 ; }
快速素数验证 $O(\cfrac{1}{3}\sqrt N)$ 1 2 3 4 5 6 7 8 bool quick_is_prime (int x) { if (x==2 ||x==3 ) return true ; if (x%6 !=1 &&x%6 !=5 ) return false ; for (int i=5 ;i<=sqrt (x);i+=6 ) if (x%i==0 ||x%(i+2 )==0 ) return false ; return true ; }
Miller-Rabbin素性测试 $O(klog^3n)$ 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 bool millerRabbin (int n, int test_time = 10 ) { if (n < 3 ) return n == 2 ; int a = n - 1 , b = 0 ; while (a % 2 == 0 ) a /= 2 , ++b; for (int i = 1 , j; i <= test_time; ++i) { int x = rand() % (n - 2 ) + 2 , v = quickPow(x, a, n); if (v == 1 || v == n - 1 ) continue ; for (j = 0 ; j < b; ++j) { v = (long long )v * v % n; if (v == n - 1 ) break ; } if (j >= b) return 0 ; } return 1 ; }
Pollard Rho素性测试 $O(n^{1/4})$ 对于每个数字检验是否是质数,是质数就输出Prime
;如果不是质数,输出它最大的质因子是哪个。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 #include <bits/stdc++.h> #define Tp template<typename Ty> #define Ts template<typename Ty,typename... Ar> #define Reg register #define RI Reg int #define RL Reg LL #define Con const #define CI Con int& #define CL Con LL& #define I inline #define W while #define LL long long #define Gmax(x,y) (x<(y)&&(x=(y))) #define abs(x) ((x)<0?-(x):(x)) #define hl_AK_NOI true using namespace std ;I LL Qmul (CL x,CL y,CL X) { RL k=(LL)((1.0L *x*y)/(1.0L *X)),t=x*y-k*X; t-=X;W(t<0 ) t+=X;return t; } class MillerRabin //MillerRabin 判素数板子{ private : #define Pcnt 12 Con int P[Pcnt]={2 ,3 ,5 ,7 ,11 ,13 ,17 ,19 ,61 ,2333 ,4567 ,24251 }; I LL Qpow (RL x,RL y,CL X) { RL t=1 ;W(y) y&1 &&(t=Qmul(t,x,X)),x=Qmul(x,x,X),y>>=1 ; return t; } I bool Check (CL x,CI p) { if (!(x%p)||Qpow(p%x,x-1 ,x)^1 ) return false ; RL k=x-1 ,t;W(!(k&1 )) { if ((t=Qpow(p%x,k>>=1 ,x))^1 &&t^(x-1 )) return false ; if (!(t^(x-1 ))) return true ; }return true ; } public : bool IsPrime (CL x) { if (x<2 ) return false ; for (RI i=0 ;i^Pcnt;++i) {if (!(x^P[i])) return true ;if (!Check(x,P[i])) return false ;} return true ; } }; class PollardRho //PollardRho 分解质因数{ private : #define Rand(x) (1LL*rand()*rand()*rand()*rand()%(x)+1) LL ans;MillerRabin MR; I LL gcd (CL x,CL y) {return y?gcd(y,x%y):x;} I LL Work (CL x,CI y) { RI t=0 ,k=1 ;RL v0=Rand(x-1 ),v=v0,d,s=1 ;W(hl_AK_NOI) { if (v=(Qmul(v,v,x)+y)%x,s=Qmul(s,abs (v-v0),x),!(v^v0)||!s) return x; if (++t==k) {if ((d=gcd(s,x))^1 ) return d;v0=v,k<<=1 ;} } } I void Resolve (RL x,RI t) { if (!(x^1 )||x<=ans) return ;if (MR.IsPrime(x)) return (void )Gmax(ans,x); RL y=x;W((y=Work(x,t--))==x);W(!(x%y)) x/=y;Resolve(x,t),Resolve(y,t); } public : I PollardRho () {srand(20050521 );} I LL GetMax (CL x) {return ans=0 ,Resolve(x,302627441 ),ans;} }P; int main () { RI Ttot;RL n,res;scanf ("%d" ,&Ttot); W(Ttot--) scanf ("%lld" ,&n),(res=P.GetMax(n))^n?printf ("%lld\n" ,res):puts ("Prime" ); return 0 ; }
素数统计 洲阁筛 $O(n^{2/3} / logn)$
用低于线性的复杂度统计 $[1,n]$ 内的素数个数。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 #include <bits/stdc++.h> #define rep(i,x,y) for (int i=(x); i<=(y); i++) #define per(i,x,y) for (int i=(x); i>=(y); i--) #define ll long long #define ld long double #define inf 1000000000 using namespace std ;#define N 350005 int p[N/10 ],tot,res[N]; bitset <N> vis; ll n;void pre (int n) { rep (i,2 ,n){ res[i]=res[i-1 ]; if (!vis[i]) p[++tot]=i,res[i]++; for (int j=1 ; j<=tot && (ll)i*p[j]<=n; j++){ vis[i*p[j]]=1 ; if (i%p[j]==0 ) break ; } } } #define M 350005 int sn,pos,cnt,last[M<<1 ]; ll g[M<<1 ],value[M<<1 ];ll cal (ll n) { cnt=0 ; sn=(ll)sqrt ((ld)n); pos=upper_bound(p+1 ,p+1 +tot,sn)-p-1 ; for (ll i=n; i>=1 ; i=n/(n/i+1 )) value[++cnt]=n/i; ll k; rep (i,1 ,cnt) g[i]=value[i],last[i]=0 ; rep (i,1 ,pos) per (j,cnt,1 ){ k=value[j]/p[i]; if (k<p[i]) break ; k=k<sn?k:cnt-n/k+1 ; g[j]-=g[k]-(i-last[k]-1 ); last[j]=i; } return res[sn]+g[cnt]-1 ; } int main () {#ifdef local freopen("test.in" ,"r" ,stdin ); freopen("test.out" ,"w" ,stdout ); #endif pre(350000 ); cin >>n; cout <<cal(n)<<endl ; return 0 ; }
Meissel Lehmer筛法 $O(n^(3/4) / log n)$
用低于线性的复杂度统计 $[L,R]$ 内的素数个数。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 #include <bits/stdc++.h> #define INF 2000000005 #define lowbit(a) ((a)&-(a)) #define FAIL -INF const long long MAXN = 6893911 ;long long p[MAXN], cnt;bool mark[MAXN];int pi[MAXN];void init () { long long i, j; for (i = 2 ; i < MAXN; i++){ if (!mark[i]) p[cnt++] = i; pi[i] = pi[i - 1 ] + !mark[i]; for (j = 0 ; p[j] * i < MAXN&&j < cnt; j++){ mark[p[j] * i] = true ; if (i%p[j] == 0 ) break ; } } } int f (long long n, int m) { if (n == 0 ) return 0 ; if (m == 0 ) return n - n / 2 ; return f(n, m - 1 ) - f(n / p[m], m - 1 ); } int Pi (long long N) ;int p2 (long long n, int m) { int ans = 0 ; for (int i = m + 1 ; (long long )p[i] * p[i] <= n; i++) ans += Pi(n / p[i]) - Pi(p[i]) + 1 ; return ans; } int p3 (long long n, int m) { int ans = 0 ; for (int i = m + 1 ; (long long )p[i] * p[i] * p[i] <= n; i++) ans += p2(n / p[i], i - 1 ); return ans; } int Pi (long long N) { if (N < MAXN) return pi[N]; int lim = pow (N, 0.25 ) + 1 ; int i; for (i = 0 ; p[i] <= lim; i++); int ans = i + f(N, i - 1 ) - 1 - p2(N, i - 1 ) - p3(N, i - 1 ); return ans; } int main () { long long L, R; scanf ("%lld %lld" , &L, &R); init(); printf ("%d" , Pi(R) - Pi(L-1 )); return 0 ; }
埃式筛法 $O(nloglogn)$
初始化:默认所有数为素数(也可以附加设定0
,1
不是素数)
从i=2
开始,遍历最近的素数(到根号n
为止)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 bool prime[MAXN];void Eratosthenes (int b) { prime[0 ] = prime[1 ] = 0 ; for (int i = 2 ; i <= b; ++i) prime[i] = 1 ; for (int i = 2 ; i*i <= b; ++i) if (prime[i]) for (int j = i*i; j <= b; j+=i){ prime[j] = 0 ; } } bool is_prime (int number) { return prime[number]; }
欧拉线筛 $O(n)$ 注意到埃式筛法质因子较多的合数会被反复标0
(累计$O(loglogn)$?),这显然是重复的工作。
而每一个合数的最小/最大质因子都是唯一的。
枚举每一个数 i
若i
是素数,先存入 素数表p_numbers
(显然素数表单增)
用J
遍历素数表(生成关于i*J
的乘法元集,并筛掉)
先筛掉 i*J
(默认筛 )
若i%J==0
(即i*J=A*J*J
)
记J’>J
是素数表中任意一个比J
更大的素数
接下来要筛的任何一个对象i*J'=(A*J)*J'=(A*J')*J
,都会在之后被一个比i
更大的i'=(A*J')
默认筛掉
而且只要被筛对象合法,i'=(A*J')<i*J'
肯定会被遍历到
而且能保证不会重复筛
假设i'*J'
之前被i
筛过,那么必然有i'*J'=(i*J)*J'=i*(J*J')
,而(J*J')
是合数,矛盾
显然无须继续遍历,提前退出 即可
欧拉筛能保证合数总是被它最大的因子
/最小的质因子
筛掉 。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 bool prime[MAXN];int p_numbers[MAXN], size = 0 ; void Euler (int Length) { prime[0 ] = prime[1 ] = 0 ; for (int i = 2 ; i <= Length; ++i) prime[i] = 1 ; for (int i = 2 ; i <= Length; ++i){ if (prime[i]){ p_numbers[size ++] = i; } for (int j = 0 ; j <size && i*p_numbers[j] < Length; j++){ prime[ i*p_numbers[j] ] = 0 ; if (i % p_numbers[j] == 0 ) break ; } } }
因数分解(唯一分解定理) 利用素数筛法筛出素数表,再遍历表寻找素因数。复杂度取决于采用的筛法。
整数的因数分解形式表示如下:
1 2 3 4 5 bool prime[MAXN]; int p_numbers[MAXN], size = 0 ; int e[MAXN]; e = {1 ,0 ,2 ,0 ,0 ,0 ,...}
因数分解:(须先进行欧拉线筛)
1 2 3 4 5 6 7 8 9 10 11 12 13 int e[MAXN];int factor (int n) { memset (e, 0 , sizeof (e)); for (int i = 0 ; i < size && n != 1 ; ++i) while (n % p_numbers[i] == 0 ){ n /= p_numbers[i]; e[i] += 1 ; } int ret = 1 ; if (n > 1 ) ret *= 2 ; _for(i, 0 , size ) ret *= e[i] + 1 ; return ret; }
积性函数 如果已知一个函数为数论函数,且 $f(1)=1$,并且满足以下条件,若对于任意的两个互质 的正整数p ,q 都满足$f(p\cdot q)=f(p)\cdot f(q)$,那么则称这个函数为积性函数 。
常见积性函数:
$\mu(n)$:莫比乌斯函数。
$\varphi(n)$:欧拉函数。表示不大于且 $n$ 与 $n$ 互质的正整数个数。
$d(n)$:约数个数。表示 $n$ 的约数的个数。
$\sigma(n)$:约数和函数。表示 $n$ 的各个约数之和。
$f(x)=x^k$:幂函数。
任意积性函数都可以线性筛。
约数个数 $d(n)$ $O(n)$ 已知唯一分解式:$n=p_1^{e_1}\times p_2^{e_2}\times p_3^{e_3}\times …\times p_n^{e_n}$。求$n$的正约数的个数。
根据乘法原理,求$n$的正约数的个数为:(利用唯一分解定理预处理)
1 2 3 4 5 6 7 8 9 10 11 12 int d[N],a[N],pri[N],tot,zhi[N];void sieve () { zhi[1 ]=d[1 ]=1 ; for (int i=2 ;i<=n;i++){ if (!zhi[i]) pri[++tot]=i,d[i]=2 ,a[i]=1 ; for (int j=1 ;j<=tot&&i*pri[j]<=n;j++){ zhi[i*pri[j]]=1 ; if (i%pri[j]) d[i*pri[j]]=d[i]*d[pri[j]],a[i*pri[j]]=1 ; else {d[i*pri[j]]=d[i]/(a[i]+1 )*(a[i]+2 );a[i*pri[j]]=a[i]+1 ;break ;} } } }
约数和 $\sigma(n)$ $O(n)$
注意可能溢出。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 int low[N],sum[N],sigma[N],pri[N],tot,zhi[N];void sieve () { zhi[1 ]=low[1 ]=sum[1 ]=sigma[1 ]=1 ; for (int i=2 ;i<=n;i++){ if (!zhi[i]) low[i]=pri[++tot]=i,sum[i]=sigma[i]=i+1 ; for (int j=1 ;j<=tot&&i*pri[j]<=n;j++){ zhi[i*pri[j]]=1 ; if (i%pri[j]==0 ) { low[i*pri[j]]=low[i]*pri[j]; sum[i*pri[j]]=sum[i]+low[i*pri[j]]; sigma[i*pri[j]]=sigma[i]/sum[i]*sum[i*pri[j]]; break ; } low[i*pri[j]]=pri[j]; sum[i*pri[j]]=pri[j]+1 ; sigma[i*pri[j]]=sigma[i]*sigma[pri[j]]; } } }
欧拉函数 $\varphi(n)$ $O(nlogn)$ 定义:欧拉函数$\varphi(n)$,表示$\leq n$的数中与$n$互质的数的数目。
欧拉函数的求解方法:
$\varphi(1) = 1$
若n是素数p的k次幂,$\varphi(n)=p^k-p^{k-1}=(p-1)p^{k-1}$
若m、n互质,$\varphi(mn)=\varphi(m)\varphi(n)$
欧拉函数的递推式:
令p为N的最小质因数,若$p^2|N$,则$\varphi(N)=\varphi(\cfrac{N}{p})\times p$;否则$\varphi(N)=\varphi(\cfrac{N}{p})\times (p-1)$。
记$n=p_1^{e_1}\times p_2^{e_2}\times p_3^{e_3}\times …\times p_n^{e_n}$,还有计算公式:($O(k)$)
1 2 3 4 5 6 7 8 9 10 11 int euler_phi (int n) { int m = int (sqrt (n + 0.5 )); int ans = n; for (int i = 2 ; i <= m; i++) if (n % i == 0 ) { ans = ans / i * (i - 1 ); while (n % i == 0 ) n /= i; } if (n > 1 ) ans = ans / n * (n - 1 ); return ans; }
线筛欧拉函数 $O(nloglogn)$ 1 2 3 4 5 6 7 8 9 10 11 int phi[MAXN];void phi_table (int n, int * phi) { for (int i = 2 ; i <= n; i++) phi[i] = 0 ; phi[1 ] = 1 ; for (int i = 2 ; i <= n; i++) if (!phi[i]) for (int j = i; j <= n; j += i) { if (!phi[j]) phi[j] = j; phi[j] = phi[j] / i * (i - 1 ); } }
欧拉定理 $a^{\varphi(m)}\equiv1(\text{mod m})$ 若$gcd(a,m)=1$,则$a^{\varphi(m)}\equiv1(\text{mod m})$ 。
扩展欧拉定理见快速幂部分。
莫比乌斯函数 $\mu(n)$ $O(n)$ 设 $n = p_1^{k_1}p_2^{k_2}\cdot \cdot \cdot p_m^{k_m}$,其中 $p$ 为素数,则莫比乌斯函数为:
可知有平方因子的数的莫比乌斯函数值为 0。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 bool v[MAXN];int mu[MAXN], p[MAXN];void Mobius () { mu[1 ] = 1 ; for (int i = 2 ; i <= 1e7 ; ++i) { if (!v[i]) mu[i] = -1 , p[++tot] = i; for (int j = 1 ; j <= tot && i <= 1e7 / p[j]; ++j) { v[i * p[j]] = 1 ; if (i % p[j] == 0 ) { mu[i * p[j]] = 0 ; break ; } mu[i * p[j]] = -mu[i]; } }
莫比乌斯反演 给定函数 $F(n)$ 和 $f(n)$ 满足:( $d|n$ 表示b遍历 $n$ 的所有因子 $d$)
则莫比乌斯反演为:
莫比乌斯反演一般来说套路比较固定,看到两个求和符号,基本上就逃不开莫比乌斯了。
求和的东西可能不同,但是一般都与 $\gcd$ 有关,枚举 $\gcd$,再使用 $\gcd = k$ 的个数 + 分块一般能 $O(n)$,
将后面的部分提前 + 积性函数预处理 + 分块一般能单组 $O(\sqrt n)$。
莫比乌斯公式 若$d=gcd(i,j)$,则
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 bool v[MAXN];int mu[MAXN], p[MAXN];void Mobius () { mu[1 ] = 1 ; for (int i = 2 ; i <= 1e7 ; ++i) { if (!v[i]) mu[i] = -1 , p[++tot] = i; for (int j = 1 ; j <= tot && i <= 1e7 / p[j]; ++j) { v[i * p[j]] = 1 ; if (i % p[j] == 0 ) { mu[i * p[j]] = 0 ; break ; } mu[i * p[j]] = -mu[i]; } }
求gcd == 1
的个数 代入莫比乌斯公式,
换序(前面两重求和 $i,j$ 可以化简快速计数),
然后就可以分块了。
求gcd == k
的个数
https://www.luogu.com.cn/blog/An-Amazing-Blog/mu-bi-wu-si-fan-yan-ji-ge-ji-miao-di-dong-xi
https://blog.sengxian.com/algorithms/mobius-inversion-formula
莫比乌斯反演得:
1 2 3 4 5 6 7 8 9 10 ll F (int n, int m, int d) { if (n > m) swap(n, m); ll ans = 0 ; n /= d, m /= d; for (int i = 1 , last = 1 ; i <= n; i = last + 1 ) { last = min (n / (n / i), m / (m / i)); ans += (ll)(sum[last] - sum[i - 1 ]) * (n / i) * (m / i); } return ans; }
求gcd == Prime
的个数 给定 $n,m$,统计在$1\leq x\leq n,1\leq y\leq m$ 范围内且 $gcd(x,y)$ 为质数的 $(x,y)$ 有多少对。$T$ 组数据。
设 $f(d)$ 为 $gcd(i,j)=d$ 的个数;$F(n)$ 为 $gcd(i,j)=kd\;\;(k\geq 1)$ 的个数。
若按照gcd == k
的做法,则复杂度将会是$O(P\sqrt n)$。
考虑优化:(令 $k=pd$)
其中,$f(x)=\sum_{p \in \text { prime }, p | x} \mu\left(\frac{x}{p}\right)$,其线筛方程为
总时间复杂度: $O(n+T\sqrt n)$。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 #include <cstdio> #include <algorithm> const int N=1e7 +5 ,lnN=15 ;int n,m,tot,p[N/lnN],mu[N],f[N];bool flg[N];void init () { mu[1 ]=flg[1 ]=1 ; for (int i=2 ;i<N;++i) { if (!flg[i]) p[++tot]=i,mu[i]=-1 ,f[i]=1 ; for (int j=1 ;j<=tot&&i*p[j]<N;++j) { int x=i*p[j]; flg[x]=1 ; if (i%p[j]==0 ) { f[x]=mu[i]; mu[x]=0 ; break ; } else { f[x]=-f[i]+mu[i]; mu[x]=-mu[i]; } } f[i]+=f[i-1 ]; } } int main () { init(); int T; for (scanf ("%d" ,&T);T--;) { scanf ("%d%d" ,&n,&m); long long ans=0 ; if (n>m) n^=m^=n^=m; for (int i=1 ,j;i<=n;i=j+1 ) { j=std ::min (n/(n/i),m/(m/i)); ans+=1L L*(f[j]-f[i-1 ])*(n/i)*(m/i); } printf ("%lld\n" ,ans); } return 0 ;
求gcd(i, j)
的幂 分块化如下:
整除分块法 $O(\sqrt n)$ 计算形如 $\sum_{i=1}^{n}\left\lfloor\cfrac{n}{i}\right\rfloor$ 的式子,复杂度为 $O(\sqrt n)$。
有时候,可能推出来的式子不一定就是一个很裸的整除分块,可能会与某些积性函数相乘,如:μ ,φ …… 这时候,我们就需要对这些函数统计一个前缀和。因为,每当我们使用整除分块跳过一个区间的时候,其所对应的函数值也跳过了一个区间。所以此时,就需要乘上那一个区间的函数值。
约数个数和 注意到:
对于 $[1,r]$ 的整数 $k$ 作为因数在 $[1,r]$ 中出现了 $\left\lfloor\cfrac{r}{k}\right\rfloor$ 次,显然对答案的贡献为 $\left\lfloor\cfrac{r}{k}\right\rfloor$。
1 2 3 4 5 6 7 8 int DC (int n) { int ret = 0 ; for (int L=1 , R; L<=n; L=R+1 ){ R = n/(n/L); ret += (R-L+1 ) * (n/L); } return ret; }
变式,计算:($d(i j)=\sum_{x | i} \sum_{y | j}[g c d(x, y)=1]$)
注意到:(记 $g(n)=\sum_{i=1}^{n}\left\lfloor\frac{n}{i}\right\rfloor$)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 #include <cstdio> #include <algorithm> const int N=5e4 +5 ;int tot,mu[N],p[N];long long s[N];bool flg[N];void init () { mu[1 ]=1 ; for (int i=2 ;i<=5e4 ;++i) { if (!flg[i]) p[++tot]=i,mu[i]=-1 ; for (int j=1 ;j<=tot&&i*p[j]<=5e4 ;++j) { flg[i*p[j]]=1 ; if (i%p[j]==0 ) { mu[i*p[j]]=0 ; break ; } else { mu[i*p[j]]=-mu[i]; } } } for (int i=1 ;i<=5e4 ;++i) mu[i]+=mu[i-1 ]; for (int x=1 ;x<=5e4 ;++x) { long long res=0 ; for (int i=1 ,j;i<=x;i=j+1 ) j=x/(x/i),res+=1L L*(j-i+1 )*(x/i); s[x]=res; } } int main () { init(); int T; for (scanf ("%d" ,&T);T--;) { int n,m; scanf ("%d%d" ,&n,&m); if (n>m) n^=m^=n^=m; long long ans=0 ; for (int i=1 ,j;i<=n;i=j+1 ) { j=std ::min (n/(n/i),m/(m/i)); ans+=1L L*(mu[j]-mu[i-1 ])*s[n/i]*s[m/i]; } printf ("%lld\n" ,ans); } return 0 ; }
杜教筛 $n^{2/3}$ 莫比乌斯函数前缀和 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 #include <bits/stdc++.h> #define For(i, l, r) for(register ll i = (l), _end_ = (ll)(r); i <= _end_; ++i) #define Fordown(i, r, l) for(register ll i = (r), _end_ = (ll)(l); i >= _end_; --i) #define Set(a, v) memset(a, v, sizeof(a)) using namespace std ;typedef long long ll;bool chkmin (int &a, int b) {return b < a ? a = b, 1 : 0 ;}bool chkmax (int &a, int b) {return b > a ? a = b, 1 : 0 ;}inline ll read () { ll x = 0 , fh = 1 ; char ch = getchar(); for (; !isdigit (ch); ch = getchar() ) if (ch == '-' ) fh = -1 ; for (; isdigit (ch); ch = getchar() ) x = (x<<1 ) + (x<<3 ) + (ch ^ '0' ); return x * fh; } void File () {#ifdef zjp_shadow freopen ("1244.in" , "r" , stdin ); freopen ("1244.out" , "w" , stdout ); #endif } const int N = 2e6 + 1e3 ;ll mu[N], Limit, prime[N], cnt; ll summu[N]; bitset <N> is_prime;void Mu_Init (int maxn) { mu[1 ] = 1 ; int res; is_prime.set (); is_prime[0 ] = is_prime[1 ] = false ; For (i, 2 , maxn) { if (is_prime[i]) { prime[++cnt] = i; mu[i] = -1 ; } For (j, 1 , maxn) { res = prime[j] * i; if (res >= maxn) break ; is_prime[res] = false ; if (i % prime[j]) mu[res] = -mu[i]; else { mu[res] = 0 ; break ; } } } For (i, 1 , maxn) summu[i] = summu[i - 1 ] + mu[i]; } ll Size; ll p1[N], p2[N], n; inline ll &App (ll x) { if (x < Size) return p1[x]; return p2[n / x]; }inline ll Mu_Sum_ (ll x) { if (x <= Limit) return summu[x]; ll &htr = App(x); if (htr != p2[0 ]) return htr; ll res = 1 , Nextx; For (i, 2 , x) { Nextx = x / (x / i); res -= 1l l * (Nextx - i + 1 ) * Mu_Sum_(x / i); i = Nextx; } return htr = res; } inline ll Mu_Sum (ll x) { Size = (int )(sqrt (x) + 0.5 ); Set(p1, -2333 ); Set(p2, -2333 ); n = x; return Mu_Sum_(x); } int main () { File (); ll a, b; cin >> a >> b; Limit = N - (int )1e3 ; Mu_Init(Limit); cout << Mu_Sum(b) - Mu_Sum(a - 1 ) << endl ; return 0 ; }
欧拉函数前缀和 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 #include <iostream> #include <string> #include <cstdio> #include <vector> #include <string.h> #include <cstring> #include <set> #include <queue> #include <algorithm> #include <math.h> #include <stdio.h> #include <map> #include <stack> #include <list> #define INF 0x3f3f3f3f #define pii pair<int, int> #define pll pair<LL, LL> #define mp make_pair #define pb push_back #define IOS ios::sync_with_stdio(false);cin.tie(0);cout.tie(0); using namespace std ;typedef long long LL;typedef unsigned long long ULL;const int MOD = 1e9 +7 ;const int HMOD = 999959 ;const int VMOD = 5000000 ;const int MAXN = 1010 ;const int INV = 500000004 ;const double eps = 1e-8 ;int Next[4 ][2 ] = {-1 , 0 , 0 , 1 , 1 , 0 , 0 , -1 };vector <pll> h[HMOD+5 ];LL phi[VMOD+5 ]; bool isp[VMOD+5 ];int pri[VMOD+5 ], cnt = 0 ;LL n; void init () { memset (isp, 0 , sizeof (isp)); phi[1 ] = 1 ; for (int i = 2 ;i <= VMOD;i++) { if (!isp[i]) { pri[++cnt] = i; phi[i] = i-1 ; } for (int j = 1 ;j <= cnt && i*pri[j] <= VMOD;j++) { isp[pri[j]*i] = 1 ; if (i%pri[j] == 0 ) { phi[i*pri[j]] = phi[i]*pri[j]%MOD; break ; } phi[i*pri[j]] = phi[i]*(pri[j]-1 )%MOD; } } for (int i = 2 ;i <= VMOD;i++) phi[i] = (phi[i-1 ]+phi[i])%MOD; } void upd (LL n, LL res) { h[n%HMOD].pb(mp(n, res)); } LL cal (LL n) { if (n <= VMOD) return phi[n]; for (pll pa: h[n%HMOD]) if (pa.first == n) return pa.second; LL ans = 0 ; for (LL l = 2 ;l <= n;) { LL x = n/l; LL r = n/x; ans = (ans + (r-l+1 )%MOD*cal(x)%MOD)%MOD; l = r+1 ; } ans = (((n%MOD*((n+1 )%MOD)%MOD)*INV%MOD-ans)%MOD+MOD)%MOD; upd(n, ans); return ans; } int main () { IOS; init(); scanf ("%lld" , &n); printf ("%lld" , cal(n)); return 0 ; }
双筛代码 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 #include <algorithm> #include <iostream> #include <cstdlib> #include <cstring> #include <cstdio> #include <cmath> #include <tr1/unordered_map> #define re register using namespace std ;using namespace std ::tr1;inline int read () { int X=0 ,w=1 ; char c=getchar(); while (c<'0' ||c>'9' ) { if (c=='-' ) w=-1 ; c=getchar(); } while (c>='0' &&c<='9' ) X=X*10 +c-'0' ,c=getchar(); return X*w; } typedef long long LL;typedef unsigned long long ULL;const int MAXN=5000000 ;int v[MAXN+10 ];int p[MAXN+10 ];LL mu[MAXN+10 ],phi[MAXN+10 ]; inline void init () { v[1 ]=mu[1 ]=phi[1 ]=1 ; re int cnt=0 ; for (re int i=2 ;i<=MAXN;++i) { if (!v[i]) p[++cnt]=i,mu[i]=-1 ,phi[i]=i-1 ; for (re int j=1 ;j<=cnt&&i*p[j]<=MAXN;++j) { v[i*p[j]]=1 ; if (i%p[j]) mu[i*p[j]]=-mu[i],phi[i*p[j]]=phi[i]*phi[p[j]]; else { mu[i*p[j]]=0 ,phi[i*p[j]]=phi[i]*p[j]; break ; } } } for (re int i=1 ;i<=MAXN;++i) mu[i]+=mu[i-1 ],phi[i]+=phi[i-1 ]; } unordered_map <int ,LL> ansmu,ansphi;inline LL S_phi (int n) { if (n<=MAXN) return phi[n]; if (ansphi[n]) return ansphi[n]; LL ans=0 ; for (re int l=2 ,r;r<2147483647 &&l<=n;l=r+1 ) r=n/(n/l),ans+=(r-l+1 )*S_phi(n/l); return ansphi[n]=(ULL)n*(n+1l l)/2l l-ans; } inline LL S_mu (int n) { if (n<=MAXN) return mu[n]; if (ansmu[n]) return ansmu[n]; LL ans=0 ; for (re int l=2 ,r;r<2147483647 &&l<=n;l=r+1 ) r=n/(n/l),ans+=(r-l+1 )*S_mu(n/l); return ansmu[n]=1l l-ans; } int main () { init(); int T=read (); while (T--) { int n=read (); printf ("%lld %lld\n" ,S_phi(n),S_mu(n)); } return 0 ; }
数位分解 将一个数按位分解。如:1234
分解成1 2 3 4
。
除余法 1 2 3 4 5 6 7 8 9 vecotr<int > split (int n) { vector <int > ret; while (n){ int t = n % 10 ; ret.push_back(t); n /= 10 ; } reverse(ret.begin (), ret.end ()); return ret; }
字符法 1 2 3 4 5 6 7 vecotr<int > split (string n) { vector <int > ret; ret.resize(n.size ()); for (int i = 0 ; i < n.size (); ++i) ret[i] = n[i] - '0' ; return ret; }
格雷码 给定一个二进制的位数n
,求出一个$0$到$2^n-1$的排列,使得相邻两项(包括头尾相邻)的二进制表达式中只有恰好一位不同。
Gray序列的第i
位为i xor (i >> 1)
。
1 2 3 4 5 6 vector <int > Gray_Create (int n) { vector <int > res; res.clear (); for (int i = 0 ; i < (1 << n); ++i) res.push_back(i ^ (i >> 1 )); return res; }
高精度整数 Java内置高精度整数。见STL4Java小节。
Python自动支持高精度 。如下:(建议优先使用Python )
1 2 a, b = map(int, input().split(' ' )) print(a + b)
进制转换 进制转换还是Python香。
1 2 3 4 5 6 7 8 9 def convert (A : str, M : int, N : int) : a = int(A, M) base = "0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ" ret = [] while a: ret.append(base[a % N]) a = a // N return '' .join(ret[::-1 ])
C++高精度类
支持加减乘除、进制转换。
进制转换也可以先转成十进制再转成任意进制。(适用于小规模数字)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 struct BigInt { int digit[MAXN]; int SIZE, BASE; BigInt(){BASE = 10 ; SIZE = 0 ;} BigInt(int base){BASE = base; SIZE = 0 ;} BigInt(const string & S, int base = 10 ){ SIZE = S.size (); BASE = base; string baseChar = "0123456789ABCDEF" ; _for(i, 0 , SIZE) digit[i] = baseChar.find (S[i]); } void overflow () { for (int i = 0 ; i < SIZE; ++i) { if (digit[i] >= BASE){ if (i < SIZE - 1 ){ digit[i + 1 ] += digit[i] / BASE; digit[i] %= BASE; } else { digit[SIZE++] = digit[i] / BASE; digit[i] %= BASE; } } } } BigInt operator +(const BigInt &A){ BigInt ret (BASE) ; int max_size = max (SIZE, A.SIZE); ret.SIZE = max_size; for (int i =0 ; i < max_size; ++i){ if (i < SIZE) ret.digit[i] += digit[i]; if (i < A.SIZE) ret.digit[i] += A.digit[i]; } ret.overflow (); return ret; } BigInt operator +(const int x){ BigInt ret = *this ; ret.digit[0 ] += x; ret.overflow (); return ret; } BigInt operator *(const BigInt &A){ BigInt ret (BASE) ; int max_size = SIZE + A.SIZE - 1 ; ret.SIZE = max_size; for (int i = 0 ; i < SIZE; ++i) for (int j = 0 ; j < A.SIZE; ++j) ret.digit[i + j] += digit[i] * A.digit[j]; ret.overflow (); return ret; } BigInt operator *(const int x){ BigInt ret = *this ; for (int i = 0 ; i < SIZE; ++i) ret.digit[i] *= x; ret.overflow (); return ret; } BigInt operator /(const int x){ BigInt ret (BASE) ; ret.SIZE = SIZE; int remainder = 0 ; for (int i = SIZE - 1 ; i >= 0 ; --i){ int t = (remainder * BASE + digit[i]) / x; int r = (remainder * BASE + digit[i]) % x; ret.digit[ret.SIZE++] = t; remainder = r; } reverse(ret.digit, ret.digit + SIZE); while (ret.digit[ret.SIZE - 1 ] == 0 && ret.SIZE > 1 ) ret.SIZE--; return ret; } int operator %(const int x){ int remainder = 0 ; for (int i = SIZE - 1 ; i >= 0 ; --i) remainder = (remainder * BASE + digit[i]) % x; return remainder; } bool operator ==(const BigInt &A){ if (SIZE == A.SIZE){ for (int i = SIZE -1 ; i >= 0 ; --i){ if (digit[i] != A.digit[i]) return false ; } return true ; } else return false ; } bool operator <(const BigInt &A){ if (SIZE == A.SIZE){ for (int i = SIZE -1 ; i >= 0 ; --i){ if (digit[i] < A.digit[i]) return true ; else if (digit[i] > A.digit[i]) return false ; } return false ; } else return SIZE < A.SIZE; } bool operator <=(const BigInt &A){ if (SIZE == A.SIZE){ for (int i = SIZE -1 ; i >= 0 ; --i){ if (digit[i] < A.digit[i]) return true ; else if (digit[i] > A.digit[i]) return false ; } return true ; } else return SIZE < A.SIZE; } }; BigInt convert_to (BigInt A, int base) { BigInt ret (base) ; while (A.SIZE > 1 || A.digit[0 ] > 0 ){ ret.digit[ret.SIZE++] = A % base; A = A / base; } return ret; } string baseChar = "0123456789ABCDEF" ; char S[MAXN];istream& operator >>(istream &in , BigInt &A){ scanf ("%s" , S); A.SIZE = 0 ; A.BASE = 10 ; for (int i = strlen (S) - 1 ; i >= 0 ; --i){ A.digit[A.SIZE++] = baseChar.find (S[i]); } return in; } ostream& operator <<(ostream &out, const BigInt& A){ for (int i = A.SIZE - 1 ; i >= 0 ; --i) printf ("%c" , baseChar[A.digit[i]]); return out; }
压位技巧 关于压位 :为了提高运算效率 ,可以用1个数位来表示多个数位,此时BASE
要设为$10^k$的形式。一般最多压缩4位。(如果digit
数组采用long long
,则可以压缩9位)
只能用来卡常数。
1 2 3 4 5 6 7 8 9 10 11 12 13 BigInt compress4 (BigInt &A) { int tmp = 0 , base = 1 ; BigInt ret (10000 ) ; for (int i = 0 ; i < A.digit.size (); ++i){ tmp += A.digit[i]*base; base *= 10 ; if ((i+1 ) % 4 == 0 ){ ret.digit.push_back(tmp); tmp = 0 ; base = 1 ; } else if (i == A.digit.size () - 1 ){ ret.digit.push_back(tmp); } } return ret; }
高精度模板2号(FFT)
支持FFT乘法。必须在C++11或更新版下使用。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 #include <cstdio> #include <iostream> #include <cmath> #include <string> #include <cstring> #include <vector> #include <algorithm> using namespace std ;const double PI = acos (-1.0 );struct Complex { double x,y; Complex(double _x = 0.0 ,double _y = 0.0 ){ x = _x; y = _y; } Complex operator -(const Complex &b)const { return Complex(x - b.x,y - b.y); } Complex operator +(const Complex &b)const { return Complex(x + b.x,y + b.y); } Complex operator *(const Complex &b)const { return Complex(x*b.x - y*b.y,x*b.y + y*b.x); } }; void change (Complex y[],int len) { int i,j,k; for (int i = 1 ,j = len/2 ;i<len-1 ;i++){ if (i < j) swap(y[i],y[j]); k = len/2 ; while (j >= k){ j = j - k; k = k/2 ; } if (j < k) j+=k; } } void fft (Complex y[],int len,int on) { change(y,len); for (int h = 2 ;h <= len;h<<=1 ){ Complex wn (cos (on*2 *PI/h),sin (on*2 *PI/h)) ; for (int j = 0 ;j < len;j += h){ Complex w (1 ,0 ) ; for (int k = j;k < j + h/2 ;k++){ Complex u = y[k]; Complex t = w*y[k + h/2 ]; y[k] = u + t; y[k + h/2 ] = u - t; w = w*wn; } } } if (on == -1 ){ for (int i = 0 ;i < len;i++){ y[i].x /= len; } } } class BigInt { #define Value(x, nega) ((nega) ? -(x) : (x)) #define At(vec, index) ((index) < vec.size() ? vec[(index)] : 0) static int absComp (const BigInt &lhs, const BigInt &rhs) { if (lhs.size () != rhs.size ()) return lhs.size () < rhs.size () ? -1 : 1 ; for (int i = lhs.size () - 1 ; i >= 0 ; --i) if (lhs[i] != rhs[i]) return lhs[i] < rhs[i] ? -1 : 1 ; return 0 ; } using Long = long long ; const static int Exp = 9 ; const static Long Mod = 1000000000 ; mutable std ::vector <Long> val; mutable bool nega = false ; void trim () const { while (val.size () && val.back() == 0 ) val.pop_back(); if (val.empty()) nega = false ; } int size () const { return val.size (); } Long &operator [](int index) const { return val[index]; } Long &back () const { return val.back(); } BigInt(int size , bool nega) : val(size ), nega(nega) {} BigInt(const std ::vector <Long> &val, bool nega) : val(val), nega(nega) {} public : friend std ::ostream &operator <<(std ::ostream &os, const BigInt &n) { if (n.size ()) { if (n.nega) putchar ('-' ); for (int i = n.size () - 1 ; i >= 0 ; --i) { if (i == n.size () - 1 ) printf ("%lld" , n[i]); else printf ("%0*lld" , n.Exp, n[i]); } } else putchar ('0' ); return os; } friend BigInt operator +(const BigInt &lhs, const BigInt &rhs) { BigInt ret (lhs) ; return ret += rhs; } friend BigInt operator -(const BigInt &lhs, const BigInt &rhs) { BigInt ret (lhs) ; return ret -= rhs; } BigInt(Long x = 0 ) { if (x < 0 ) x = -x, nega = true ; while (x >= Mod) val.push_back(x % Mod), x /= Mod; if (x) val.push_back(x); } BigInt(const char *s) { int bound = 0 , pos; if (s[0 ] == '-' ) nega = true , bound = 1 ; Long cur = 0 , pow = 1 ; for (pos = strlen (s) - 1 ; pos >= Exp + bound - 1 ; pos -= Exp, val.push_back(cur), cur = 0 , pow = 1 ) for (int i = pos; i > pos - Exp; --i) cur += (s[i] - '0' ) * pow , pow *= 10 ; for (cur = 0 , pow = 1 ; pos >= bound; --pos) cur += (s[pos] - '0' ) * pow , pow *= 10 ; if (cur) val.push_back(cur); } BigInt &operator =(const char *s){ BigInt n(s); *this = n; return n; } BigInt &operator =(const Long x){ BigInt n(x); *this = n; return n; } friend std ::istream &operator >>(std ::istream &is, BigInt &n){ string s; is >> s; n=(char *)s.data(); return is; } BigInt &operator +=(const BigInt &rhs) { const int cap = std ::max (size (), rhs.size ()) + 1 ; val.resize(cap); int carry = 0 ; for (int i = 0 ; i < cap - 1 ; ++i) { val[i] = Value(val[i], nega) + Value(At(rhs, i), rhs.nega) + carry, carry = 0 ; if (val[i] >= Mod) val[i] -= Mod, carry = 1 ; else if (val[i] < 0 ) val[i] += Mod, carry = -1 ; } if ((val.back() = carry) == -1 ) { nega = true , val.pop_back(); bool tailZero = true ; for (int i = 0 ; i < cap - 1 ; ++i) { if (tailZero && val[i]) val[i] = Mod - val[i], tailZero = false ; else val[i] = Mod - 1 - val[i]; } } trim(); return *this ; } friend BigInt operator -(const BigInt &rhs) { BigInt ret (rhs) ; ret.nega ^= 1 ; return ret; } BigInt &operator -=(const BigInt &rhs) { rhs.nega ^= 1 ; *this += rhs; rhs.nega ^= 1 ; return *this ; } friend BigInt operator *(const BigInt &lhs, const BigInt &rhs) { int len=1 ; BigInt ll=lhs,rr=rhs; ll.nega = lhs.nega ^ rhs.nega; while (len<2 *lhs.size ()||len<2 *rhs.size ())len<<=1 ; ll.val.resize(len),rr.val.resize(len); Complex x1[len],x2[len]; for (int i=0 ;i<len;i++){ Complex nx(ll[i],0.0),ny(rr[i],0.0); x1[i]=nx; x2[i]=ny; } fft(x1,len,1 ); fft(x2,len,1 ); for (int i = 0 ; i < len; i++) x1[i] = x1[i] * x2[i]; fft( x1 , len , -1 ); for (int i = 0 ; i < len; i++) ll[i] = int ( x1[i].x + 0.5 ); for (int i = 0 ; i < len; i++){ ll[i+1 ]+=ll[i]/Mod; ll[i]%=Mod; } ll.trim(); return ll; } friend BigInt operator *(const BigInt &lhs, const Long &x){ BigInt ret=lhs; bool negat = ( x < 0 ); Long xx = (negat) ? -x : x; ret.nega ^= negat; ret.val.push_back(0 ); ret.val.push_back(0 ); for (int i = 0 ; i < ret.size (); i++) ret[i]*=xx; for (int i = 0 ; i < ret.size (); i++){ ret[i+1 ]+=ret[i]/Mod; ret[i] %= Mod; } ret.trim(); return ret; } BigInt &operator *=(const BigInt &rhs) { return *this = *this * rhs; } BigInt &operator *=(const Long &x) { return *this = *this * x; } friend BigInt operator /(const BigInt &lhs, const BigInt &rhs) { static std ::vector <BigInt> powTwo{BigInt(1 )}; static std ::vector <BigInt> estimate; estimate.clear (); if (absComp(lhs, rhs) < 0 ) return BigInt(); BigInt cur = rhs; int cmp; while ((cmp = absComp(cur, lhs)) <= 0 ) { estimate.push_back(cur), cur += cur; if (estimate.size () >= powTwo.size ()) powTwo.push_back(powTwo.back() + powTwo.back()); } if (cmp == 0 ) return BigInt(powTwo.back().val, lhs.nega ^ rhs.nega); BigInt ret = powTwo[estimate.size () - 1 ]; cur = estimate[estimate.size () - 1 ]; for (int i = estimate.size () - 1 ; i >= 0 && cmp != 0 ; --i) if ((cmp = absComp(cur + estimate[i], lhs)) <= 0 ) cur += estimate[i], ret += powTwo[i]; ret.nega = lhs.nega ^ rhs.nega; return ret; } friend BigInt operator /(const BigInt &num,const Long &x){ bool negat = ( x < 0 ); Long xx = (negat) ? -x : x; BigInt ret; Long k = 0 ; ret.val.resize( num.size () ); ret.nega = (num.nega ^ negat); for (int i = num.size () - 1 ;i >= 0 ; i--){ ret[i] = ( k * Mod + num[i]) / xx; k = ( k * Mod + num[i]) % xx; } ret.trim(); return ret; } bool operator ==(const BigInt &rhs) const { return nega == rhs.nega && val == rhs.val; } bool operator !=(const BigInt &rhs) const { return nega != rhs.nega || val != rhs.val; } bool operator >=(const BigInt &rhs) const { return !(*this < rhs); } bool operator >(const BigInt &rhs) const { return !(*this <= rhs); } bool operator <=(const BigInt &rhs) const { if (nega && !rhs.nega) return true ; if (!nega && rhs.nega) return false ; int cmp = absComp(*this , rhs); return nega ? cmp >= 0 : cmp <= 0 ; } bool operator <(const BigInt &rhs) const { if (nega && !rhs.nega) return true ; if (!nega && rhs.nega) return false ; return (absComp(*this , rhs) < 0 ) ^ nega; } void swap (const BigInt &rhs) const { std ::swap(val, rhs.val); std ::swap(nega, rhs.nega); } }; BigInt ba,bb; int main () { cin >>ba>>bb; cout << ba + bb << '\n' ; cout << ba - bb << '\n' ; cout << ba * bb << '\n' ; BigInt d; cout << (d = ba / bb) << '\n' ; cout << ba - d * bb << '\n' ; return 0 ; }
高精度模板3号
Karatsuba 乘法 $O(n^{log_23})\approx O(n^{1.585})$ 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 int *karatsuba_polymul (int n, int *a, int *b) { if (n <= 32 ) { int *r = new int [n * 2 + 1 ](); for (int i = 0 ; i <= n; ++i) for (int j = 0 ; j <= n; ++j) r[i + j] += a[i] * b[j]; return r; } int m = n / 2 + 1 ; int *r = new int [m * 4 + 1 ](); int *z0, *z1, *z2; z0 = karatsuba_polymul(m - 1 , a, b); z2 = karatsuba_polymul(n - m, a + m, b + m); for (int i = 0 ; i + m <= n; ++i) a[i] += a[i + m]; for (int i = 0 ; i + m <= n; ++i) b[i] += b[i + m]; z1 = karatsuba_polymul(m - 1 , a, b); for (int i = 0 ; i + m <= n; ++i) a[i] -= a[i + m]; for (int i = 0 ; i + m <= n; ++i) b[i] -= b[i + m]; for (int i = 0 ; i <= (m - 1 ) * 2 ; ++i) z1[i] -= z0[i]; for (int i = 0 ; i <= (n - m) * 2 ; ++i) z1[i] -= z2[i]; for (int i = 0 ; i <= (m - 1 ) * 2 ; ++i) r[i] += z0[i]; for (int i = 0 ; i <= (m - 1 ) * 2 ; ++i) r[i + m] += z1[i]; for (int i = 0 ; i <= (n - m) * 2 ; ++i) r[i + m * 2 ] += z2[i]; delete [] z0; delete [] z1; delete [] z2; return r; } void karatsuba_mul (int a[], int b[], int c[]) { int *r = karatsuba_polymul(LEN - 1 , a, b); memcpy (c, r, sizeof (int ) * LEN); for (int i = 0 ; i < LEN - 1 ; ++i) if (c[i] >= 10 ) { c[i + 1 ] += c[i] / 10 ; c[i] %= 10 ; } delete [] r; }
复数 1 2 3 4 5 6 7 8 struct Complex { double r ,i; Complex() {r = 0 , i = 0 ;}; Complex(double real, double imag): r(real), i(imag) {} Complex operator + (Complex A){ return Complex(r + A.r, i + A.i);} Complex operator - (Complex A){ return Complex(r - A.r, i - A.i);} Complex operator * (Complex A){ return Complex(r * A.r - i * A.i, r * A.i + i * A.r);} };
快速傅里叶 FFT $O(nlogn)$ 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 const double PI = acos (-1 );void FFT (Complex *a, int lim) { if (lim == 1 ) return ; Complex a0[lim >> 1 ], a1[lim >> 1 ]; for (int i = 0 ; i < lim; i+=2 ) a0[i >> 1 ] = a[i], a1[i >> 1 ] = a[i + 1 ]; FFT(a0, lim >> 1 ); FFT(a1, lim >> 1 ); Complex wn = Complex(cos (2.0 * PI / lim), sin (2.0 * PI / lim)); Complex w = Complex(1 , 0 ); for (int k = 0 ; k < (lim >> 1 ); ++k){ Complex t = w * a1[k]; a[k] = a0[k] + t; a[k + (lim >> 1 )] = a0[k] - t; w *= wn; } }
迭代版FFT模板 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 const double PI = acos (-1 );int rev[MAXN], len, lim = 1 ;void FFT (Complex *a, int opt = 1 ) { for (int i = 0 ; i < lim; ++i) if (i < rev[i]) swap(a[i], a[rev[i]]); for (int dep = 1 ; dep <= log2(lim); ++dep){ int m = 1 << dep; Complex wn = Complex(cos (2.0 * PI / m), opt * sin (2.0 * PI / m)); for (int k = 0 ; k < lim; k += m){ Complex w = Complex(1 , 0 ); for (int j = 0 ; j < m / 2 ; ++j){ Complex t = w * a[k + j + m/2 ]; Complex u = a[k + j]; a[k + j] = u + t; a[k + j + m/2 ] = u - t; w = w * wn; } } } if (opt == -1 ) for (int i = 0 ; i < lim; ++i) a[i].r /= lim; } int main () { int n; scanf ("%d" , &n); while (lim <= n) lim <<=1 , len++; for (int i =0 ; i < lim; ++i) rev[i] = (rev[i>>1 ] >> 1 ) | ((i&1 ) << (len-1 )); }
多项式 多项式乘法FFT 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 #include <bits/stdc++.h> using namespace std ;#define MAXN 5000005 struct Complex { double r ,i; Complex() {r = 0 , i = 0 ;}; Complex(double real, double imag): r(real), i(imag) {} Complex operator + (Complex A){ return Complex(r + A.r, i + A.i);} Complex operator - (Complex A){ return Complex(r - A.r, i - A.i);} Complex operator * (Complex A){ return Complex(r * A.r - i * A.i, r * A.i + i * A.r);} } F[MAXN], G[MAXN]; const double PI = acos (-1 );int rev[MAXN], len, lim = 1 ;void FFT (Complex *a, int opt = 1 ) { for (int i = 0 ; i < lim; ++i) if (i < rev[i]) swap(a[i], a[rev[i]]); for (int dep = 1 ; dep <= log2(lim); ++dep){ int m = 1 << dep; Complex wn = Complex(cos (2.0 * PI / m), opt * sin (2.0 * PI / m)); for (int k = 0 ; k < lim; k += m){ Complex w = Complex(1 , 0 ); for (int j = 0 ; j < m / 2 ; ++j){ Complex t = w * a[k + j + m/2 ]; Complex u = a[k + j]; a[k + j] = u + t; a[k + j + m/2 ] = u - t; w = w * wn; } } } if (opt == -1 ) for (int i = 0 ; i < lim; ++i) a[i].r /= lim; } int main () { int n, m, t; scanf ("%d %d" , &n, &m); for (int i = 0 ; i <= n; ++i) scanf ("%d" , &t), F[i].r = t; for (int i = 0 ; i <= m; ++i) scanf ("%d" , &t), G[i].r = t; while (lim <= n + m) lim <<=1 , len++; for (int i =0 ; i < lim; ++i) rev[i] = (rev[i>>1 ] >> 1 ) | ((i&1 ) << (len-1 )); FFT(F, 1 ); FFT(G, 1 ); for (int i = 0 ; i <= lim; ++i) F[i] = F[i] * G[i]; FFT(F, -1 ); for (int i = 0 ; i <= n + m; ++i) printf ("%d " , (int )(F[i].r + 0.5 )); return 0 ; }
多项式模乘 MTT
与FFT不同的是,结果需要取模。使用FFT会炸精度,而NTT会因为模数的性质而失去作用。
MTT有2种方法,一种是三模数NTT,然后是拆系数FFT。
其中NTT精度优秀,但常数较大,而FFT则相反。
三模数NTT 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 #include <algorithm> #include <cstdio> #include <cstring> int mod;namespace Math { inline int pw (int base, int p, const int mod) { static int res; for (res = 1 ; p; p >>= 1 , base = static_cast <long long > (base) * base % mod) if (p & 1 ) res = static_cast <long long > (res) * base % mod; return res; } inline int inv (int x, const int mod) { return pw(x, mod - 2 , mod); } } const int mod1 = 998244353 , mod2 = 1004535809 , mod3 = 469762049 , G = 3 ;const long long mod_1_2 = static_cast <long long > (mod1) * mod2;const int inv_1 = Math::inv(mod1, mod2), inv_2 = Math::inv(mod_1_2 % mod3, mod3);struct Int { int A, B, C; explicit inline Int () { } explicit inline Int(int __num) : A(__num), B(__num), C(__num) { } explicit inline Int(int __A, int __B, int __C) : A(__A), B(__B), C(__C) { } static inline Int reduce (const Int &x) { return Int(x.A + (x.A >> 31 & mod1), x.B + (x.B >> 31 & mod2), x.C + (x.C >> 31 & mod3)); } inline friend Int operator + (const Int &lhs, const Int &rhs) { return reduce(Int(lhs.A + rhs.A - mod1, lhs.B + rhs.B - mod2, lhs.C + rhs.C - mod3)); } inline friend Int operator - (const Int &lhs, const Int &rhs) { return reduce(Int(lhs.A - rhs.A, lhs.B - rhs.B, lhs.C - rhs.C)); } inline friend Int operator * (const Int &lhs, const Int &rhs) { return Int(static_cast <long long > (lhs.A) * rhs.A % mod1, static_cast <long long > (lhs.B) * rhs.B % mod2, static_cast <long long > (lhs.C) * rhs.C % mod3); } inline int get () { long long x = static_cast <long long > (B - A + mod2) % mod2 * inv_1 % mod2 * mod1 + A; return (static_cast <long long > (C - x % mod3 + mod3) % mod3 * inv_2 % mod3 * (mod_1_2 % mod) % mod + x) % mod; } } ; #define maxn 100072 namespace Poly {#define N (maxn << 1) int lim, s, rev[N]; Int Wn[N | 1 ]; inline void init (int n) { s = -1 , lim = 1 ; while (lim < n) lim <<= 1 , ++s; for (register int i = 1 ; i < lim; ++i) rev[i] = rev[i >> 1 ] >> 1 | (i & 1 ) << s; const Int t (Math::pw(G, (mod1 - 1 ) / lim, mod1), Math::pw(G, (mod2 - 1 ) / lim, mod2), Math::pw(G, (mod3 - 1 ) / lim, mod3)) ; *Wn = Int(1 ); for (register Int *i = Wn; i != Wn + lim; ++i) *(i + 1 ) = *i * t; } inline void NTT (Int *A, const int op = 1 ) { for (register int i = 1 ; i < lim; ++i) if (i < rev[i]) std ::swap(A[i], A[rev[i]]); for (register int mid = 1 ; mid < lim; mid <<= 1 ) { const int t = lim / mid >> 1 ; for (register int i = 0 ; i < lim; i += mid << 1 ) { for (register int j = 0 ; j < mid; ++j) { const Int W = op ? Wn[t * j] : Wn[lim - t * j]; const Int X = A[i + j], Y = A[i + j + mid] * W; A[i + j] = X + Y, A[i + j + mid] = X - Y; } } } if (!op) { const Int ilim (Math::inv(lim, mod1), Math::inv(lim, mod2), Math::inv(lim, mod3)) ; for (register Int *i = A; i != A + lim; ++i) *i = (*i) * ilim; } } #undef N } int n, m;Int A[maxn << 1 ], B[maxn << 1 ]; int main () { scanf ("%d%d%d" , &n, &m, &mod); ++n, ++m; for (int i = 0 , x; i < n; ++i) scanf ("%d" , &x), A[i] = Int(x % mod); for (int i = 0 , x; i < m; ++i) scanf ("%d" , &x), B[i] = Int(x % mod); Poly::init(n + m); Poly::NTT(A), Poly::NTT(B); for (int i = 0 ; i < Poly::lim; ++i) A[i] = A[i] * B[i]; Poly::NTT(A, 0 ); for (int i = 0 ; i < n + m - 1 ; ++i) { printf ("%d" , A[i].get ()); putchar (i == n + m - 2 ? '\n' : ' ' ); } return 0 ; }
拆系数FFT 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 #include <bits/stdc++.h> typedef long long int64;typedef long double D;int MOD;const int MAXN = 524288 ;const D PI = acos (-1 );struct complex { D real, imag; complex () { real = imag = 0 ; } complex (D x): real(x), imag(0 ) {} complex (D x, D y): real(x), imag(y) {} inline complex conj () { return complex (real, -imag); } inline complex operator +(complex rhs) const { return complex (real + rhs.real, imag + rhs.imag); } inline complex operator -(complex rhs) const { return complex (real - rhs.real, imag - rhs.imag); } inline complex operator *(complex rhs) const { return complex (real * rhs.real - imag * rhs.imag, imag * rhs.real + real * rhs.imag); } inline complex operator *=(complex rhs) { return (*this ) = (*this ) * rhs; } friend inline complex operator *(D x, complex cp) { return complex (x * cp.real, x * cp.imag); } inline complex operator /(D x) const { return complex (real / x, imag / x); } inline complex operator /=(D x) { return (*this ) = (*this ) / x; } friend inline complex operator /(D x, complex cp) { return x * cp.conj() / (cp.real * cp.real - cp.imag * cp.imag); } inline complex operator /(complex rhs) const { return (*this ) * rhs.conj() / (rhs.real * rhs.real - rhs.imag * rhs.imag); } inline complex operator /=(complex rhs) { return (*this ) = (*this ) / rhs; } inline D length () { return sqrt (real * real + imag * imag); } }; inline complex get_omega (int len, bool inv) { return inv ? complex (std ::cos (2 *PI / len), -std ::sin (2 *PI / len)) :complex (std ::cos (2 *PI / len), std ::sin (2 *PI / len)); } inline void fft (int len, complex * A, bool inv = false ) { static int R[MAXN+10 ]; for (int i = 0 ; i < len; i++) R[i] = ((R[i>>1 ]>>1 ) | (len >> (i&1 ))) & (len-1 ); for (int i = 0 ; i < len; i++) if (R[i] > i) std ::swap(A[i], A[R[i]]); for (int step = 1 ; step < len; step <<= 1 ) for (int i = 0 ; i < len; i += (step <<1 )) { complex omega = 1 , t = 0 , rt = get_omega(step <<1 , inv); for (int j = 0 ; j < step ; j++, omega *= rt) { t = A[i+j+step ] * omega; A[i+j+step ] = A[i+j] - t; A[i+j] = A[i+j] + t; } } if (inv) for (int i = 0 ; i < len; i++) A[i] /= len; } void mtt (int deg, int *lhs, int *rhs, int *ret) { static complex A[MAXN+10 ], B[MAXN+10 ], C[MAXN+10 ], D[MAXN+10 ], E[MAXN+10 ], F[MAXN+10 ], G[MAXN+10 ], H[MAXN+10 ]; int len; for (len = 1 ; len <= 2 *deg; len<<=1 ); for (int i = 0 ; i < len; i++) { lhs[i] %= MOD; rhs[i] %= MOD; A[i] = lhs[i] >> 15 ; B[i] = lhs[i] & 0x7fff ; C[i] = rhs[i] >> 15 ; D[i] = rhs[i] & 0x7fff ; } fft(len, A); fft(len, B); fft(len, C); fft(len, D); for (int i = 0 ; i < len; i++) { E[i] = A[i] * C[i]; F[i] = B[i] * C[i]; G[i] = A[i] * D[i]; H[i] = B[i] * D[i]; } fft(len, E, true ); fft(len, F, true ); fft(len, G, true ); fft(len, H, true ); for (int i = 0 ; i < len; i++) ret[i] = (((int64(round(E[i].real)) % MOD)<<30 ) % MOD + ((int64(round(F[i].real)) % MOD)<<15 ) % MOD + ((int64(round(G[i].real)) % MOD)<<15 ) % MOD + int64(round(H[i].real)) % MOD) % MOD; } int deg, n, m, F[MAXN+10 ], G[MAXN+10 ], H[MAXN+10 ];inline int readint () { int f=1 , r=0 ; char c=getchar(); while (!isdigit (c)) { if (c=='-' )f=-1 ; c=getchar(); } while (isdigit (c)) { r=r*10 +c-'0' ; c=getchar(); } return f*r; } int main () { n = readint(); m = readint(); MOD = readint(); for (int i = 0 ; i <= n; i++) F[i] = readint(); for (int i = 0 ; i <= m; i++) G[i] = readint(); mtt(std ::max (n, m), F, G, H); for (int i = 0 ; i <= n + m; i++) printf ("%d " , H[i]); }
高精度整数FFT乘法
溢出须进位,与多项式不同。
FFT的局限性:当任意一个数的位数很少时,朴素乘法的效率高。
Python自带$O(n^{1.585})$的Karatsuba高精度乘法。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 const double PI = acos (-1 );#define MAXN 1000005 #define _for(i,a,b) for(int i =(a); i < (b); ++i) int rev[MAXN], len, lim = 1 ;struct Complex { }; struct BigInt { }; void FFT (Complex *a, int opt = 1 ) { } BigInt operator * (const BigInt &A, const BigInt &B){ int n = A.SIZE - 1 , m = B.SIZE - 1 ; while (lim <= n + m) lim <<=1 , len++; _for(i, 0 , lim) rev[i] = (rev[i>>1 ] >> 1 ) | ((i&1 ) << (len-1 )); _rep(i, 0 , lim) F[i].r = F[i].i = G[i].r = G[i].i = 0 ; _rep(i, 0 , n) F[i].r = A.digit[i]; _rep(i, 0 , m) G[i].r = B.digit[i]; FFT(F, 1 ); FFT(G, 1 ); _rep(i, 0 , lim) F[i] = F[i] * G[i]; FFT(F, -1 ); BigInt ret (A.BASE) ; ret.SIZE = n + m + 1 ; _rep(i, 0 , n + m) ret.digit[i] = (int )(F[i].r + 0.5 ); ret.overflow (); return ret; }
多项式开根
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 #include <bits/stdc++.h> using namespace std ;const int maxn = 1 << 20 , mod = 998244353 ;int a[maxn], b[maxn], g[maxn], gg[maxn];int qpow (int x, int y) { int ans = 1 ; while (y) { if (y & 1 ) { ans = 1L L * ans * x % mod; } x = 1L L * x * x % mod; y >>= 1 ; } return ans; } int inv2 = qpow(2 , mod - 2 );inline void change (int *f, int len) { for (int i = 1 , j = len >> 1 ; i < len - 1 ; i++) { if (i < j) { swap(f[i], f[j]); } int k = len >> 1 ; while (j >= k) { j -= k; k >>= 1 ; } if (j < k) { j += k; } } } inline void NTT (int *f, int len, int type) { change(f, len); for (int q = 2 ; q <= len; q <<= 1 ) { int nxt = qpow(3 , (mod - 1 ) / q); for (int i = 0 ; i < len; i += q) { int w = 1 ; for (int k = i; k < i + (q >> 1 ); k++) { int x = f[k]; int y = 1L L * w * f[k + (q >> 1 )] % mod; f[k] = (x + y) % mod; f[k + (q >> 1 )] = (x - y + mod) % mod; w = 1L L * w * nxt % mod; } } } if (type == -1 ) { reverse(f + 1 , f + len); int iv = qpow(len, mod - 2 ); for (int i = 0 ; i < len; i++) { f[i] = 1L L * f[i] * iv % mod; } } } inline void inv (int deg, int *f, int *h) { if (deg == 1 ) { h[0 ] = qpow(f[0 ], mod - 2 ); return ; } inv(deg + 1 >> 1 , f, h); int len = 1 ; while (len < deg << 1 ) { len <<= 1 ; } copy(f, f + deg, gg); fill (gg + deg, gg + len, 0 ); NTT(gg, len, 1 ); NTT(h, len, 1 ); for (int i = 0 ; i < len; i++) { h[i] = 1L L * (2 - 1L L * gg[i] * h[i] % mod + mod) % mod * h[i] % mod; } NTT(h, len, -1 ); fill (h + deg, h + len, 0 ); } int n, t[maxn];inline void sqrt (int deg, int *f, int *h) { if (deg == 1 ) { h[0 ] = 1 ; return ; } sqrt (deg + 1 >> 1 , f, h); int len = 1 ; while (len < deg << 1 ) { len <<= 1 ; } fill (g, g + len, 0 ); inv(deg, h, g); copy(f, f + deg, t); fill (t + deg, t + len, 0 ); NTT(t, len, 1 ); NTT(g, len, 1 ); NTT(h, len, 1 ); for (int i = 0 ; i < len; i++) { h[i] = 1L L * inv2 * (1L L * h[i] % mod + 1L L * g[i] * t[i] % mod) % mod; } NTT(h, len, -1 ); fill (h + deg, h + len, 0 ); } int main () { cin >> n; for (int i = 0 ; i < n; i++) { scanf ("%d" , &a[i]); } sqrt (n, a, b); for (int i = 0 ; i < n; i++) { printf ("%d " , b[i]); } return 0 ; }
多项式方程求根 牛顿迭代法 开根通常需要结合二分的思想。既可以使用经典二分法,也可以使用二阶收敛的牛顿法。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 struct Polynomial { #define Item double #define MAXORDER 105 int n; Item digit[MAXORDER]; Item& operator [] (const int x){ return digit[x]; } Item operator () (Item x) { return calculate(x); } Polynomial() {memset (digit, 0 , sizeof (digit));} Polynomial(int N){ memset (digit, 0 , sizeof (digit)); n = N + 1 ; digit[N] = 1 ; } Polynomial derivative () { Polynomial ret (n - 2 ) ; _rep(i, 1 , n - 1 ) ret[i - 1 ] = (Item)(i) * digit[i]; return ret; } Item calculate (Item x) { Item ret = 0 , base = 1 ; _rep(i, 0 , n - 1 ){ ret += base * digit[i]; base *= x; } return ret; } Item Newton (Item N) { Item x = N, pre = INF; Polynomial f = *this ; f[0 ] = -N; Polynomial f_1 = f.derivative(); if (f_1(x) < eps) x += eps; while (abs (x - pre) > eps){ pre = x; x = x - f(x) / f_1(x); } return x; } }; Polynomial p (2 ) ;int main () { _rep(i, 0 , 2 ) cout << p.derivative().digit[i] << " " ; cout << endl ; cout << p.calculate(1 ) <<endl ; cout << p(2 ) <<endl ; cout << p.derivative().calculate(3 ) <<endl ; cout << p.Newton(0 ) <<endl ; return 0 ; }
高精度开方-牛顿法 Java 1 2 3 4 5 6 7 8 9 10 11 12 public static BigInteger isqrtNewton (BigInteger n) { BigInteger a = BigInteger.ONE.shiftLeft(n.bitLength() / 2 ); boolean p_dec = false ; for (;;) { BigInteger b = n.divide(a).add(a).shiftRight(1 ); if (a.compareTo(b) == 0 || a.compareTo(b) < 0 && p_dec) break ; p_dec = a.compareTo(b) > 0 ; a = b; } return a; }
指数法开根 1 2 3 4 5 6 7 #define Item long double #define TEST 10 Item solve (Item x, int n) { long long ret = exp (log (x) / n); while (pow (ret + 1 , n) <= x) ret ++; return ret; }
多项式重建 拉格朗日插值 给定n个点,求过这n个点的多项式。
输入n, k,n个点,求 f(k) % 998244353
。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 #include <cstdio> #define N 2005 using namespace std ;const int mod = 998244353 ;int n, k;int x[N], y[N];inline int fpm (int bs, int mi) { int res = 1 ; while (mi) { if (mi & 1 ) res = 1L L * res * bs % mod; bs = 1L L * bs * bs % mod, mi >>= 1 ; } return res; } inline int Sub (int x, int y) { return (x -= y) < 0 ? x + mod : x; }inline int inv (int x) { return fpm(x, mod - 2 ); }int lagrange () { int res = 0 , l; for (int i = 0 ; i < n; ++i) { l = 1 ; for (int j = 0 ; j < n; ++j) { if (i == j) continue ; l = 1L L * l * Sub(k, x[j]) % mod * inv(Sub(x[i], x[j])) % mod; } res = (res + 1L L * l * y[i]) % mod; } return res; } int main () { scanf ("%d%d" , &n, &k); for (int i = 0 ; i < n; ++i) scanf ("%d%d" , &x[i], &y[i]); printf ("%d\n" , lagrange()); return 0 ; }
分数/有理数
只保存最简分数,实现了加减乘除。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 struct Fraction { int num, den; Fraction(int num = 0 , int den = 1 ){ if (den < 0 ){ num = -num; den = -den;} assert(den != 0 ); int g = __gcd(abs (num), den); this ->num = num / g; this ->den = den / g; } Fraction operator +(const Fraction &A) const { return Fraction(num * A.den + den * A.num, den * A.den); } Fraction operator -(const Fraction &A) const { return Fraction(num * A.den - den * A.num, den * A.den); } Fraction operator *(const Fraction &A) const { return Fraction(num * A.num, den * A.den); } Fraction operator /(const Fraction &A) const { return Fraction(num * A.den, den * A.num); } bool operator <(const Fraction &A) const { return num * A.den < den * A.num; } bool operator ==(const Fraction &A) const { return num * A.den == den * A.num; } }; istream& operator >>(istream &in , Fraction &A){ int a, b; cin >> a >> b; A = Fraction(a, b); return in; } ostream& operator <<(ostream &out, Fraction A){ cout << A.num << "/" << A.den; return out; }
Python分数类 1 2 3 4 5 6 7 8 from fractions import Fractiona = Fraction(16 , -10 ) b = Fraction("-8/6" ) c = Fraction('1.414213 \t\n' ) d = a.as_integer_ratio() from math import *e = floor(Fraction(355 , 113 ))
有理数取余 给出一个有理数 $c=\frac{a}{b}$,求 $c \bmod 19260817$ 的值。
费马小定理求逆元 $b^{-1}$,则 $c=a\times b^{-1}$。【核心在于除法变乘法 】
1 2 3 4 5 6 7 8 9 10 11 12 13 14 MOD = 19260817 def p (a,b) : t=a; res=1 while (b!=0 ): if (b%2 ==1 ): res=res*t%MOD t=t*t%MOD b=b//2 return res a=int(input()) b=int(input()) if (b==0 ): print("Angry!" ) else : print((a*p(b,MOD))%MOD)
分数规划 有n个物品,每个物品有两个权值 $a$ 和 $b$,求一组 $w_i\in \{0,1\}$,使得$\cfrac{\sum a_i\times w_i}{\sum b_i\times w_i}$ 最大。
分数规划问题的通用方法是二分。
于是问题转化为如何求 $\sum w_i\times (a_i-mid\times b_i)$ 的最大值。(最小值同理)
我们只需将 $a_i-mid\times b_i$ 作为第 $i$ 个物品的权值,问题便可贪心解决。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 #include <algorithm> #include <cmath> #include <cstdio> #include <cstdlib> #include <cstring> #include <iostream> using namespace std ;inline int read () { int X = 0 , w = 1 ; char c = getchar(); while (c < '0' || c > '9' ) { if (c == '-' ) w = -1 ; c = getchar(); } while (c >= '0' && c <= '9' ) X = X * 10 + c - '0' , c = getchar(); return X * w; } const int N = 100000 + 10 ;const double eps = 1e-6 ;int n;double a[N], b[N];inline bool check (double mid) { double s = 0 ; for (int i = 1 ; i <= n; ++i) if (a[i] - mid * b[i] > 0 ) s += a[i] - mid * b[i]; return s > 0 ; } int main () { n = read (); for (int i = 1 ; i <= n; ++i) a[i] = read (); for (int i = 1 ; i <= n; ++i) b[i] = read (); double L = 0 , R = 1e9 ; while (R - L > eps) { double mid = (L + R) / 2 ; if (check(mid)) L = mid; else R = mid; } printf ("%.6lf\n" , L); return 0 ; }
有约束的分数规划 额外要求:分母 $\sum w_i\times b_i \geq W$。
我们只需将 $b_i$ 视为第 $i$ 个物品的重量,将 $a_i-mid\times b_i$ 作为第 $i$ 个物品的权值,便转化为背包问题。
1 2 3 4 5 6 7 8 9 10 double f[1010 ];inline bool check (double mid) { for (int i = 1 ; i <= W; i++) f[i] = -1e9 ; for (int i = 1 ; i <= n; i++) for (int j = W; j >= 0 ; j--) { int k = min (W, j + b[i]); f[k] = max (f[k], f[j] + a[i] - mid * b[i]); } return f[W] > 0 ; }
矩阵 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 typedef double Item; struct Matrix { int n, m; Item a[MAXN][MAXN]; Matrix (int N = 0 , int M = 0 ){ clear (); n = N; m = M; } void set (int N = 0 , int M = 0 ) { clear (); n = N; m = M; } void clear () { n = m = 0 ; memset (a, 0 , sizeof (a)); } Item* operator [] (int x){ return a[x]; } void set_to_identity (int N = 0 ) { set (N, N); _for(i, 0 , n) _for(j, 0 , m) a[i][j] = (i == j); } void transpose () { _for(i, 0 , n) _for(j, i+1 , m) swap(a[i][j], a[j][i]); swap(n, m); } Matrix operator + (const Matrix &b) const { Matrix tmp (n, m) ; _for(i, 0 , n) _for(j, 0 , m) tmp.a[i][j] = a[i][j] + b.a[i][j]; return tmp; } Matrix operator - (const Matrix &b) const { Matrix tmp (n, m) ; _for(i, 0 , n) _for(j, 0 , m) tmp.a[i][j] = a[i][j] - b.a[i][j]; return tmp; } Matrix operator * (const int &b) const { Matrix tmp (n, m) ; _for(i, 0 , n) _for(j, 0 , m) tmp.a[i][j] = a[i][j] * b; return tmp; } Matrix operator * (const Matrix &b) const { Matrix tmp (n, b.m) ; _for(i, 0 , n) _for(j, 0 , b.m) _for(k, 0 , m) tmp.a[i][j] += a[i][k] * b.a[k][j]; return tmp; } Item gauss () { Item ans = 1 ; for (int i = 0 ; i < n; i++) { int sid = -1 ; for (int j = i; j < n; j++) if (abs (a[j][i]) > eps) { sid = j; break ; } if (sid == -1 ) continue ; if (sid != i) { for (int j = 0 ; j < n; j++) { swap(a[sid][j], a[i][j]); ans = -ans; } } for (int j = i + 1 ; j < n; j++) { Item ratio = a[j][i] / a[i][i]; for (int k = 0 ; k < n; k++) { a[j][k] -= a[i][k] * ratio; } } } for (int i = 0 ; i < n; i++) ans *= a[i][i]; return ans; } }; Matrix qPower (Matrix a, int n) { Matrix ans(n, n), base = a; _for(i, 0 , n) ans.a[i][i] = 1 ; while (n){ if (n & 1 ) ans = ans * base; base = base * base; n >>= 1 ; } return ans; } int readint () { int x; scanf ("%d" , &x); return x; } istream& operator >>(istream &in , Matrix &A){ _for(i, 0 , A.n) _for(j, 0 , A.m){ cin >> A[i][j]; } return in; } ostream& operator <<(ostream &out, Matrix& A){ _for(i, 0 , A.n) _for(j, 0 , A.m) cout << A[i][j] << " \n" [j == A.m - 1 && i != A.n - 1 ]; return out; }
高斯消元/求行列式
高斯消元矩阵需要double
浮点数。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 #include <algorithm> #include <cassert> #include <cmath> #include <cstdio> #include <cstring> #include <iostream> using namespace std ;#define MOD 100000007 #define eps 1e-7 struct matrix { static const int maxn = 20 ; int n, m; double mat[maxn][maxn]; matrix() { memset (mat, 0 , sizeof (mat)); } void print () { cout << "MATRIX " << n << " " << m << endl ; for (int i = 0 ; i < n; i++) { for (int j = 0 ; j < m; j++) { cout << mat[i][j] << "\t" ; } cout << endl ; } } void random (int n) { this ->n = n; this ->m = n; for (int i = 0 ; i < n; i++) for (int j = 0 ; j < n; j++) mat[i][j] = rand() % 100 ; } void initSquare () { this ->n = 4 ; this ->m = 4 ; memset (mat, 0 , sizeof (mat)); mat[0 ][1 ] = mat[0 ][3 ] = 1 ; mat[1 ][0 ] = mat[1 ][2 ] = 1 ; mat[2 ][1 ] = mat[2 ][3 ] = 1 ; mat[3 ][0 ] = mat[3 ][2 ] = 1 ; mat[0 ][0 ] = mat[1 ][1 ] = mat[2 ][2 ] = mat[3 ][3 ] = -2 ; this ->n--; this ->m--; } double gauss () { double ans = 1 ; for (int i = 0 ; i < n; i++) { int sid = -1 ; for (int j = i; j < n; j++) if (abs (mat[j][i]) > eps) { sid = j; break ; } if (sid == -1 ) continue ; if (sid != i) { for (int j = 0 ; j < n; j++) { swap(mat[sid][j], mat[i][j]); ans = -ans; } } for (int j = i + 1 ; j < n; j++) { double ratio = mat[j][i] / mat[i][i]; for (int k = 0 ; k < n; k++) { mat[j][k] -= mat[i][k] * ratio; } } } for (int i = 0 ; i < n; i++) ans *= mat[i][i]; return abs (ans); } }; int main () { srand(1 ); matrix T; T.initSquare(); T.print (); double ans = T.gauss(); T.print (); cout << ans << endl ; }
矩阵求逆 double
型求逆:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 #include <bits/stdc++.h> #define eps 1e-8 using namespace std ;double a[101 ][102 ];int n;bool GUASS () { for (int i=1 ;i<=n;i++){ for (int j=i;j<=n;j++){ if (fabs (a[j][i])>eps){ for (int k=1 ;k<=n;k++){ swap(a[i][k],a[j][k]); } swap(a[i][n+1 ],a[j][n+1 ]); break ; } } for (int j=1 ;j<=n;j++){ if (i==j) continue ; double tmp=a[j][i]/a[i][i]; for (int k=i;k<=n;k++){ a[j][k]-=a[i][k]*tmp; } a[j][n+1 ]-=a[i][n+1 ]*tmp; } } for (int i=1 ;i<=n;i++){ bool lala=0 ; for (int j=1 ;j<=n;j++){ if (a[i][j]!=0 ){ lala=1 ; } } if (lala==0 ){ cout <<"No Solution" ; return 0 ; } } return 1 ; } int main () { cin >>n; for (int i=1 ;i<=n;i++){ for (int j=1 ;j<=n+1 ;j++){ scanf ("%lf" ,&a[i][j]); } } if (GUASS()){ for (int i=1 ;i<=n;i++){ printf ("%.2lf\n" ,a[i][n+1 ]/a[i][i]); } } }
int
型求逆:(模意义下)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 #include <bits/stdc++.h> #define p 1000000007 #define inc(i,a,b) for(register int i=a;i<=b;i++) using namespace std ;int n;int a[1000 ][1000 ];long long KSM (long long a,long long b) { long long res=1 ; while (b){ if (b&1 ) res=res*a%p; a=a*a%p; b/=2 ; } return res%p; } int Gauss () { inc(i,1 ,n){ inc(j,i,n) if (a[j][i]) {swap(a[j],a[i]);break ;} if (a[i][i]==0 ){cout <<"No Solution" ; return 0 ;} long long tmp=KSM(a[i][i],p-2 ); inc(j,i,2 *n) a[i][j]=a[i][j]*tmp%p; inc(j,1 ,n){ if (i==j) continue ; long long tmp2=a[j][i]%p; inc(k,i,2 *n) a[j][k]=(a[j][k]-tmp2*a[i][k]%p+p)%p; } } return 1 ; } int main () { cin >>n; inc(i,1 ,n){ inc(j,1 ,n) scanf ("%d" ,&a[i][j]); a[i][i+n]=1 ; } if (Gauss()){ inc(i,1 ,n){ inc(j,n+1 ,2 *n){ printf ("%d " ,a[i][j]); } printf ("\n" ); } } }
矩阵面积并 扫描线算法 求 $n$ 个矩形的面积并(覆盖)。
分析一下题目,可以用一条直线从左向右扫描所有竖直线段,每次计算这些线段的并集的长度,然后把各个子矩形的面积相加。问题便转化为:
维护一些区间,支持增加和删除区间,每次操作后求这些区间并集的长度。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 #include <stdio.h> #include <iostream> #include <algorithm> #define lson (x << 1) #define rson (x << 1 | 1) using namespace std ;const int MAXN = 1e6 + 10 ;typedef long long ll;int n, cnt = 0 ;ll x1, y1, x2, y2, X[MAXN << 1 ]; struct ScanLine { ll l, r, h; int mark; bool operator < (const ScanLine &rhs) const { return h < rhs.h; } } line [MAXN << 1 ]; struct SegTree { int l, r, sum; ll len; } tree[MAXN << 2 ]; void build_tree (int x, int l, int r) { tree[x].l = l, tree[x].r = r; tree[x].len = 0 ; tree[x].sum = 0 ; if (l == r) return ; int mid = (l + r) >> 1 ; build_tree(lson, l, mid); build_tree(rson, mid + 1 , r); return ; } void pushup (int x) { int l = tree[x].l, r = tree[x].r; if (tree[x].sum ) tree[x].len = X[r + 1 ] - X[l]; else tree[x].len = tree[lson].len + tree[rson].len; } void edit_tree (int x, ll L, ll R, int c) { int l = tree[x].l, r = tree[x].r; if (X[r + 1 ] <= L || R <= X[l]) return ; if (L <= X[l] && X[r + 1 ] <= R) { tree[x].sum += c; pushup(x); return ; } edit_tree(lson, L, R, c); edit_tree(rson, L, R, c); pushup(x); } int main () { scanf ("%d" , &n); for (int i = 1 ; i <= n; i++) { scanf ("%lli %lli %lli %lli" , &x1, &y1, &x2, &y2); X[2 * i - 1 ] = x1, X[2 * i] = x2; line [2 * i - 1 ] = (ScanLine) {x1, x2, y1, 1 }; line [2 * i] = (ScanLine) {x1, x2, y2, -1 }; } n <<= 1 ; sort(line + 1 , line + n + 1 ); sort(X + 1 , X + n + 1 ); int tot = unique(X + 1 , X + n + 1 ) - X - 1 ; build_tree(1 , 1 , tot - 1 ); ll ans = 0 ; for (int i = 1 ; i < n ; i++) { edit_tree(1 , line [i].l, line [i].r, line [i].mark); ans += tree[1 ].len * (line [i + 1 ].h - line [i].h); } printf ("%lli" , ans); return 0 ; }
矩阵周长并 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 #include <iostream> #include <stdio.h> #include <algorithm> #define lson (x << 1) #define rson (x << 1 | 1) using namespace std ;const int MAXN = 2e4 ;int n, X[MAXN << 1 ];int x1, y1, x2, y2, pre = 0 ; struct ScanLine { int l, r, h, mark; if (h == rhs.h) return mark > rhs.mark; return h < rhs.h; } line [MAXN]; struct SegTree { int l, r, sum, len, c; bool lc, rc; } tree[MAXN << 2 ]; void build_tree (int x, int l, int r) { tree[x].l = l, tree[x].r = r; tree[x].lc = tree[x].rc = false ; tree[x].sum = tree[x].len = 0 ; tree[x].c = 0 ; if (l == r) return ; int mid = (l + r) >> 1 ; build_tree(lson, l, mid); build_tree(rson, mid + 1 , r); } void pushup (int x) { int l = tree[x].l, r = tree[x].r; if (tree[x].sum) { tree[x].len = X[r + 1 ] - X[l]; tree[x].lc = tree[x].rc = true ; tree[x].c = 1 ; } else { tree[x].len = tree[lson].len + tree[rson].len; tree[x].lc = tree[lson].lc, tree[x].rc = tree[rson].rc; tree[x].c = tree[lson].c + tree[rson].c; if (tree[lson].rc && tree[rson].lc) tree[x].c -= 1 ; } } void edit_tree (int x, int L, int R, int c) { int l = tree[x].l, r = tree[x].r; if (X[l] >= R || X[r + 1 ] <= L) return ; if (L <= X[l] && X[r + 1 ] <= R) { tree[x].sum += c; pushup(x); return ; } edit_tree(lson, L, R, c); edit_tree(rson, L, R, c); pushup(x); } ScanLine make_line (int l, int r, int h, int mark) { ScanLine res; res.l = l, res.r = r, res.h = h, res.mark = mark; return res; } int main () { scanf ("%d" , &n); for (int i = 1 ; i <= n; i++) { scanf ("%d %d %d %d" , &x1, &y1, &x2, &y2); line [i * 2 - 1 ] = make_line(x1, x2, y1, 1 ); line [i * 2 ] = make_line(x1, x2, y2, -1 ); X[i * 2 - 1 ] = x1, X[i * 2 ] = x2; } n <<= 1 ; sort(line + 1 , line + n + 1 ); sort(X + 1 , X + n + 1 ); int tot = unique(X + 1 , X + n + 1 ) - X - 1 ; build_tree(1 , 1 , tot - 1 ); int res = 0 ; for (int i = 1 ; i < n; i++) { edit_tree(1 , line [i].l, line [i].r, line [i].mark); res += abs (pre - tree[1 ].len); pre = tree[1 ].len; res += 2 * tree[1 ].c * (line [i + 1 ].h - line [i].h); } res += line [n].r - line [n].l; printf ("%d" , res); return 0 ; }
组合数学 组合数 C(n,m)
DP法求组合数 $O(mn)$ 状态转移方程:($C(n,m) = C_n^m$)
初始条件:$C(i,0)=C(i,i)=1$。
1 2 3 4 5 6 7 8 int C[MAXN][MAXN];int cal_c (int n, int m) { for (int i = 0 ; i <= n; ++i) C[i][0 ] = C[i][i] = 1 ; for (int i = 1 ; i <= n; ++i) for (int j = 1 ; j < i; ++j) C[i][j] = C[i-1 ][j-1 ] + C[i-1 ][j]; return C[n][m]; }
记忆化:
1 2 3 4 5 6 #define _c(n,m) (c[n][m] ? c[n][m] : c[n][m] = C(n,m)) int c[MAXN][MAXN];int C (int n, int m) { if (m == 0 || m == n) return 1 ; return _c(n-1 , m) + _c(n-1 , m-1 ); }
线性求组合数/计算二项式系数 $O(n)$ 根据:$C(n,k) = \cfrac{n-k+1}{k}C(n,k-1)$。初始条件:$C(n,0)=1$。
1 2 3 4 5 6 7 int C[MAXN];int cal_c (int n, int m) { C[0 ] = 1 ; for (int i = 1 ; i <= m; ++i) C[i] = C[i-1 ] * (n-i+1 ) / i; return C[m]; }
卢卡斯定理 C(n,m) % p
$O(logn+p)$
适用范围:mod是不大的素数($p<10^5$)的大组合数计算
Lucas 定理用于求解大规模组合数 取模的问题,其中 p 必须为素数 。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 #include <bits/stdc++.h> #define N 100010 using namespace std ;typedef long long ll;ll a[N]; int p;ll pow (ll y,int z,int p) { y%=p;ll ans=1 ; for (int i=z;i;i>>=1 ,y=y*y%p)if (i&1 )ans=ans*y%p; return ans; } ll C (ll n,ll m) { if (m>n)return 0 ; return ((a[n]*pow (a[m],p-2 ,p))%p*pow (a[n-m],p-2 ,p)%p); } ll Lucas (ll n,ll m) { if (!m)return 1 ; return C(n%p,m%p)*Lucas(n/p,m/p)%p; } inline int read () { int f=1 ,x=0 ;char ch; do {ch=getchar();if (ch=='-' )f=-1 ;}while (ch<'0' ||ch>'9' ); do {x=x*10 +ch-'0' ;ch=getchar();}while (ch>='0' &&ch<='9' ); return f*x; } int main () { int T=read (); while (T--){ int n=read (),m=read ();p=read (); a[0 ]=1 ; for (int i=1 ;i<=p;i++)a[i]=(a[i-1 ]*i)%p; cout <<Lucas(n+m,n)<<endl ; } }
扩展卢卡斯定理 C(n,m) % p
$O(log^2n\times p)$ Lucas 定理中对于模数p要求必须为素数,那么对于p不是素数 的情况,就需要用到 exLucas 定理。
适用范围:$p<10^9$。
p可以为任何正整数。
原理:把mod数分解质因数,计算p1^t1 , p2^t2 …. pm^tm,然后用中国剩余定理合并。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 #include <bits/stdc++.h> #define ll long long using namespace std ;inline ll fmi (ll a,ll b,ll p) { ll ans=1 ; a%=p; while (b>0 ) { if (b&1 ) ans=(ans*a)%p; a=(a*a)%p; b>>=1 ; } ans%=p; return ans; } inline ll exgcd (ll a,ll b,ll &x,ll &y) { if (!b) {x=1 ,y=0 ;return a;} ll k=exgcd(b,a%b,x,y); ll z=x;x=y,y=z-a/b*y; return k; } inline ll inv (ll a,ll b) { ll x,y; ll g=exgcd(a,b,x,y); return g==1 ?(x%b+b)%b:-1 ; } inline ll mul (ll n,ll pi,ll pk) { if (!n) return 1 ; ll ans=1 ; if (n/pk) { for (ll i=2 ;i<=pk;i++) if (i%pi) ans=ans*i%pk; ans=fmi(ans,n/pk,pk); } for (ll i=2 ;i<=n%pk;i++) if (i%pi) ans=ans*i%pk; return ans*mul(n/pi,pi,pk)%pk; } inline ll C (ll n,ll m,ll p,ll pi,ll pk) { if (n<m) return 0 ; ll a=mul(n,pi,pk); ll b=mul(m,pi,pk); ll c=mul(n-m,pi,pk),k=0 ,ans; for (ll i=n;i;i/=pi) k+=i/pi; for (ll i=m;i;i/=pi) k-=i/pi; for (ll i=n-m;i;i/=pi) k-=i/pi; ans=a*inv(b,pk)%pk*inv(c,pk)%pk*fmi(pi,k,pk)%pk; ans=ans*(p/pk)%p*inv(p/pk,pk)%p; return ans; } inline ll exLucas (ll n,ll m,ll p) { ll ans=0 ,x=p,t=sqrt (p); for (ll i=2 ;i<=t;i++) if (x%i==0 ) { ll pk=1 ; while (x%i==0 ) x/=i,pk*=i; ans=(ans+C(n,m,p,i,pk))%p; } if (x>1 ) ans=(ans+C(n,m,p,x,x))%p; return ans; } int main () { ll n, m, p; scanf ("%lld%lld%lld" ,&n,&m,&p); printf ("%lld\n" ,exLucas(n,m,p)); return 0 ; }
错排公式 $O(n)$ 所有元素均不在其应在的位置的组合数。$F[n]$表示n个数错排方案的总数。
1 2 3 4 5 6 7 int F[MAXN];int f (int n) { F[1 ] = 0 ; F[2 ] = 1 ; for (int i = 3 ; i <= n; ++i) F[i] = (n-1 ) * F[i-1 ] + (n-1 ) * F[n-2 ]; return F[n]; }
康托展开 $O(nlogn)$ 求 $1\sim N$ 的一个给定全排列在所有 $1\sim N$ 全排列中的排名。(取模)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 #include <bits/stdc++.h> using namespace std ;#define MAXN 1000005 #define rgt register #define mod 998244353 int N, a[MAXN], fac, c[MAXN], ans;char *p;inline void read ( rgt int &x ) { x = 0 ; while ( !isdigit (*p) ) ++p; while ( isdigit (*p) ) x = x * 10 + ( *p & 15 ), ++p; } int main () { scanf ( "%d" , &N ), fac = 1 ; p = new char [N * 8 + 100 ], fread( p, 1 , N * 8 + 100 , stdin ); for ( rgt int i = N; i; --i ) read (a[i]); for ( rgt int i = 1 , s, j; i <= N; ++i ){ for ( s = 0 , j = a[i]; j; j -= j & -j ) s += c[j]; ans = ( ans + 1l l * fac * s ) % mod, fac = 1l l * fac * i % mod; for ( j = a[i]; j <= N; j += j & -j ) ++c[j]; } printf ( "%d\n" , ans + 1 ); return 0 ; }
逆康托展开 $O(nlogn)$
已知排名求排列。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 static const int FAC[] = {1 , 1 , 2 , 6 , 24 , 120 , 720 , 5040 , 40320 , 362880 }; void decantor (int x, int n) { vector <int > v; vector <int > a; for (int i=1 ;i<=n;i++) v.push_back(i); for (int i=m;i>=1 ;i--){ int r = x % FAC[i-1 ]; int t = x / FAC[i-1 ]; x = r; sort(v.begin (),v.end ()); a.push_back(v[t]); v.erase(v.begin ()+t); } }
卡特兰数 Cat(n) $O(n)$ 1 2 3 4 5 6 7 8 int n;long long f[101 ];int Cat (int n) { f[0 ] = 1 ; for (int i = 1 ; i <= n; i++) f[i] = f[i - 1 ] * (4 * i - 2 ) / (i + 1 ); return f[n]; }
斐波那契数 $O(n)$
数列找规律网 (输入前几项,即可找出递推公式)。
1 2 3 4 5 6 7 8 9 10 pair<int, int> fib(int n) { // 返回值是一个二元组(Fn, Fn+1) if (n == 0 ) return {0 , 1 }; auto p = fib(n >> 1 ); int c = p.first * (2 * p.second - p.first); int d = p.first * p.first + p.second * p.second; if (n & 1 ) return {d, c + d}; else return {c, d}; }
2-SAT问题 2-SAT 问题的目标是给每个变量赋值使得所有条件得到满足。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 #include <cstdio> #define min(a,b) (a<b?a:b) const int maxn=1e6 +10 ;struct EDGE { int t,next; }edge[maxn<<2 ]; int head[maxn<<1 ],low[maxn<<1 ],dfn[maxn<<1 ],stack [maxn<<1 ],scc[maxn<<1 ];bool in[maxn<<1 ];int cur,n,tim,top,cnt,m;void add (int u,int v) { edge[++cur].t=v; edge[cur].next=head[u]; head[u]=cur; } void tarjan (int u) { low[u]=dfn[u]=++tim; stack [top++]=u;in[u]=true ; for (int i=head[u]; i; i=edge[i].next) { int v=edge[i].t; if (!dfn[v]) { tarjan(v); low[u]=min (low[u],low[v]); }else if (in[v]) low[u]=min (low[u],dfn[v]); } if (low[u]==dfn[u]) { int v; cnt++; do { v=stack [--top]; in[v]=false ; scc[v]=cnt; }while (u!=v); } } bool two_SAT () { for (int i=1 ; i<=2 *n; i++) if (!dfn[i]) tarjan(i); for (int i=1 ; i<=n; i++) if (scc[i]==scc[i+n]) return false ; return true ; } int main () { scanf ("%d%d" ,&n,&m); for (int i=1 ; i<=m; i++) { int a,b,aval,bval; scanf ("%d%d%d%d" ,&a,&aval,&b,&bval); int nota=aval^1 ,notb=bval^1 ; add(a+nota*n,b+bval*n); add(b+notb*n,a+aval*n); } if (two_SAT()) { printf ("POSSIBLE\n" ); for (int i=1 ; i<=n; i++) printf ("%d " ,scc[i]>scc[i+n]); }else printf ("IMPOSSIBLE" ); return 0 ; }
博弈论 巴什博奕 只有一堆n个物品,两个人轮流从中取物,规定每次最少取一个,最多取m个,最后取光者为胜。
1 2 3 4 5 6 7 8 9 10 #include <iostream> using namespace std ;int main () { int n,m; while (cin >>n>>m) if (n%(m+1 )==0 ) cout <<"后手必胜" <<endl ; else cout <<"先手必胜" <<endl ; return 0 ; }
威佐夫博弈 有两堆各若干的物品,两人轮流从其中一堆取至少一件物品,至多不限,或从两堆中同时取相同件物品,规定最后取完者胜利。
若两堆物品的初始值为(x,y),且x<y,则另z=y-x;
记w=(int)[((sqrt(5)+1)/2)*z ];
若w=x,则先手必败,否则先手必胜。
1 2 3 4 5 6 7 8 9 10 11 ll n,m,w; int jud; int main () { scanf ("%lld%lld" ,&n,&m); if (m>n) swap(n,m); w=n-m; jud=(int )(((sqrt (5.0 )+1.0 )/2.0 )*w); if (jud==m) printf ("0" ); else printf ("1" ); return 0 ; }
Nimm博弈 地上有 $n$ 堆石子(每堆石子数量小于 $10^4$),每人每次可从任意一堆石子里取出任意多枚石子扔掉,可以取完,不能不取。每次只能从一堆里取。最后没石子可取的人就输了。假如甲是先手,且告诉你这 $n$ 堆石子的数量,他想知道是否存在先手必胜的策略。
把每堆物品数全部异或起来,如果得到的值为0,那么先手必败,否则先手必胜。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 #include <cstdio> #include <cmath> #include <iostream> using namespace std ;int main () { int n,ans,temp; while (cin >>n) { temp=0 ; for (int i=0 ;i<n;i++) { cin >>ans; temp^=ans; } if (temp==0 ) cout <<"后手必胜" <<endl ; else cout <<"先手必胜" <<endl ; } return 0 ; }
斐波那契博弈 有一堆物品,两人轮流取物品,先手最少取一个,至多无上限,但不能把物品取完,之后每次取的物品数不能超过上一次取的物品数的二倍且至少为一件,取走最后一件物品的人获胜。
先手胜当且仅当n不是斐波那契数(n为物品总数)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 #include <iostream> #include <string.h> #include <stdio.h> using namespace std ; const int N = 55 ; int f[N]; void Init () { f[0 ] = f[1 ] = 1 ; for (int i=2 ;i<N;i++) f[i] = f[i-1 ] + f[i-2 ]; } int main () { Init(); int n; while (cin >>n) { if (n == 0 ) break ; bool flag = 0 ; for (int i=0 ;i<N;i++) { if (f[i] == n) { flag = 1 ; break ; } } if (flag) puts ("Second win" ); else puts ("First win" ); } return 0 ; }
约瑟夫问题 n 个人标号0,1,…,n。逆时针站一圈,从0号开始,每一次从当前的人逆时针数 k 个,然后让这个人出局。问最后剩下的人是谁。
线性算法 $O(n)$ 1 2 3 4 5 int josephus (int n, int k) { int res = 0 ; for (int i = 1 ; i <= n; ++i) res = (res + k) % i; return res; }
对数算法 $O(klogn)$ 1 2 3 4 5 6 7 8 9 10 11 12 int josephus (int n, int k) { if (n == 1 ) return 0 ; if (k == 1 ) return n - 1 ; if (k > n) return (josephus(n - 1 , k) + k) % n; int res = josephus(n - n / k, k); res -= n % k; if (res < 0 ) res += n; else res += res / (k - 1 ); return res; }
图论 图模板 1 2 3 4 struct edge { int a, b, w; }; vector <int > edge[MAXN];
带权边模板 1 2 3 4 struct Edge { int from, to, dist; Edge(int u, int v, int d): from(u), to(v), dist(d) {} };
(树)节点模板 动态指针型:
1 2 3 4 struct node { int data; node* lc, rc; }
静态数组型:
1 2 3 4 5 int size = 0 ;struct node { int data; int lc, rc; } Node[MAXN];
二叉树遍历 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 void vis (node* root) { cout << root->data; } void preorder (node* root) { if (root == NULL ) return ; vis(root); preorder(root->lc); preorder(root->rc); } void inorder (node* root) { if (root == NULL ) return ; inorder(root->lc); vis(root); inorder(root->rc); } void postorder (node* root) { if (root == NULL ) return ; postorder(root->lc); postorder(root->rc); vis(root); } void layerorder (node* root) { queue <node*> Q; Q.push(root); while (!Q.empty()){ node* now = Q.front(); Q.pop(); vis(now); if (now->lc) Q.push(now->lc); if (now->rc) Q.push(now->rc); } }
二叉排序树 BST 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 struct node { int data; node *lc, *rc; node(int x): data(x), lc(NULL ), rc(NULL ) {} }; struct BST { node* root = NULL ; void insert (int x) { root = Insert(x, root); } node* Insert (int x, node* cur) { if (cur == NULL ) return cur = new node(x); if (x < cur->data) cur->lc = Insert(x, cur->lc); else cur->rc = Insert(x, cur->rc); return cur; } node* find (int x) { node* ret = root; while (ret != NULL && x != ret->data){ if (x < ret->data) ret = ret->lc; else ret = ret->rc; } return ret; } bool erase (int x) { return false ; } }; int main () { BST b; b.insert(7 ); b.insert(5 ); b.insert(6 ); preorder(b.root); cout <<endl ; inorder(b.root); cout <<endl ; postorder(b.root); return 0 ; }
树同构哈希
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 #include <bits/stdc++.h> using namespace std ;const long long maxn=1001 ;long long ans[maxn][maxn],n,m,head[maxn],last[maxn],Next[maxn],tot,x;void add (int x,int y) { last[++tot]=y;Next[tot]=head[x];head[x]=tot; } long long Hash (int x,int f) { long long q[maxn],ans=maxn,top=0 ; for (int i=head[x];i;i=Next[i]) if (last[i]!=f) q[++top]=Hash(last[i],x); sort(q+1 ,q+top+1 ); for (int i=1 ;i<=top;i++) ans=ans*2333 +q[i]; return ans*2333 +maxn+1 ; } int main () { cin >>m; for (int i=1 ;i<=m;i++) { tot=0 ;memset (head,0 ,sizeof (head)); cin >>n; for (int j=1 ;j<=n;j++) { cin >>x; if (x!=0 )add(x,j),add(j,x); } for (int j=1 ;j<=n;j++) ans[i][j]=Hash(j,0 ); sort(ans[i]+1 ,ans[i]+n+1 ); for (int j=1 ,k=0 ;j<=i;j++) { while (k<=n) if (ans[i][++k]!=ans[j][k]) break ; if (k>n){printf ("%d\n" ,j);break ;} } } return 0 ; }
优先级搜索 PFS 广度优先搜索 BFS 邻接矩阵版:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 int n, G{MAXN}[MAXN];bool inq[MAXN] = {false }; void BFS (int u) { queue <int > Q; Q.push(u); inq[u] = true ; while (!Q.empty()){ int u = Q.front(); Q.pop(); EXECUTE(u); for (int v = 0 ; v < n; ++v) if (!inq[v] && G[u][v]){ Q.push(v); inq[v] = true ; } } }
邻接表版:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 int n; vector <int > Adj[MAXN]; bool inq[MAXN] = {false }; void BFS (int u) { queue <int > Q; Q.push(u); inq[u] = true ; while (!Q.empty()){ int u = Q.front(); Q.pop(); EXECUTE(u); for (int i = 0 ; i < Adj[u].size (); ++v) int v = Adj[u][i]; if (!vis[v]){ Q.push(v); inq[v] = true ; } } }
BFS森林:
1 2 3 4 5 void BFSTrave () { for (int u = 0 ; u < n; ++u) if (!inq[u]) BFS(u); }
深度优先搜索 DFS 原理伪代码:
1 2 3 4 5 void DFS (int p) { visited[p] = true ; for (int to : edge[p]) if (!visited[to]) DFS(to); }
邻接矩阵版:
1 2 3 4 5 6 7 8 9 int n, G[MAXN][MAXN]; bool vis[MAXV] = {false }; void DFS (int u, int depth) { vis[u] = true ; EXECUTE(u); for (int v = 0 ; v < n; ++v) if (!vis[v] && G[u][v]) DFS(v, depth + 1 ); }
邻接表版:
1 2 3 4 5 6 7 8 9 10 11 12 int n; vector <int > Adj[MAXN]; bool vis[MAXV] = {false }; void DFS (int u, int depth) { vis[u] = true ; EXECUTE(u); for (int i = 0 ; i < Adj[u].size (); ++v){ int v = Adj[u][i]; if (!vis[v]) DFS(v, depth + 1 ); } }
DFS森林:
1 2 3 4 5 void DFSTrave () { for (int u = 0 ; u < n; ++u) if (!vis[u]) DFS(u, 1 ); }
欧拉路径/欧拉回路* 如果图G中的一个路径包括每个边恰好一次,则该路径称为欧拉路径(Euler path)。
如果一个回路是欧拉路径,则称为欧拉回路(Euler circuit)。
例:构造一个$2^n$位的循环数组 ,满足:
每次从中截取n位,能够遍历所有的n位2进制数。
如:构造0110
数组,将遍历:01
,11
,10
,00
。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 int allOne;vector <bool > vis[2 ];string ans;int twoPow (int x) { return 1 << x; } void dfs (int u, int dex = 0 ) { for (int r = 0 ; r < dex; ++r) cout <<" " ; cout <<"节点" <<bitset <3>(u)<<":" <<endl ; for (int i = 0 ; i <= 1 ; ++i) if (!vis[i][u]) { int v = ((u << 1 ) | i) & allOne; vis[i][u] = 1 ; dfs(v, dex + 1 ); ans.push_back('0' + i); cout <<char ('0' + i)<<endl ; } } string getAnswer (int n) { allOne = twoPow(n - 1 ) - 1 ; ans = "" ; for (int i = 0 ; i < 2 ; ++i) vis[i].resize(twoPow(n - 1 ), 0 ); dfs(0 ); return ans; } int main () { int n; scanf ("%d" , &n); cout << getAnswer(n) << endl ; return 0 ; }
最小生成树 MST Kruskal算法 $O(e\text{log}e)$ 见数据结构部分-并查集uf_set
。Kruskal算法可用于求最小支撑树MST。
提前退出可求解最小生成森林(n-k次连边对应k棵树的k-森林)。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 struct Edge { int from, to, dist; Edge(int u, int v, int d): from(u), to(v), dist(d) {} bool operator < (Edge x) const { return dist > x.dist; } }; int n, m;priority_queue<Edge> PQ; uf_set U (MAXNODE) ; int Kruskal () { int ret = 0 , cnt = 0 ; while (!PQ.empty()){ Edge t = PQ.top(); PQ.pop(); if (U.find (t.from) != U.find (t.to)){ U.Union(t.from, t.to); ret += t.dist; cnt ++; } } if (cnt != n - 1 ) return INF; return ret; } int main () { cin >> n >> m; _for(i, 0 , m){ int x, y, z; cin >> x >> y >> z; PQ.push(Edge(x, y, z)); } int ret = Kruskal(); if (ret == INF) cout << "orz" <<endl ; else cout << ret <<endl ; return 0 ; }
Prim算法 $O(nlogn)$ 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 #include <cstdio> #include <queue> #include <cstring> #include <algorithm> #define R register int using namespace std ;int k,n,m,cnt,sum,ai,bi,ci,head[5005 ],dis[5005 ],vis[5005 ];struct Edge { int v,w,next; }e[400005 ]; void add (int u,int v,int w) { e[++k].v=v; e[k].w=w; e[k].next=head[u]; head[u]=k; } typedef pair <int ,int > pii;priority_queue <pii,vector <pii>,greater<pii> > q; void prim () { dis[1 ]=0 ; q.push(make_pair(0 ,1 )); while (!q.empty()&&cnt<n) { int d=q.top().first,u=q.top().second; q.pop(); if (vis[u]) continue ; cnt++; sum+=d; vis[u]=1 ; for (R i=head[u];i!=-1 ;i=e[i].next) if (e[i].w<dis[e[i].v]) dis[e[i].v]=e[i].w,q.push(make_pair(dis[e[i].v],e[i].v)); } } int main () { memset (dis,127 ,sizeof (dis)); memset (head,-1 ,sizeof (head)); scanf ("%d%d" ,&n,&m); for (R i=1 ;i<=m;i++) { scanf ("%d%d%d" ,&ai,&bi,&ci); add(ai,bi,ci); add(bi,ai,ci); } prim(); if (cnt==n)printf ("%d" ,sum); else printf ("orz" ); }
最小斯坦纳树 最小斯坦纳树,就是要花费最小的代价,连通给定的 k 个关键点 ,这是一个组合优化问题。
给定图中一个点集,求出包含该点集的最小边权和的连通图。
输入格式
第一行:三个整数 n,m,k,表示 G 的结点数、边数和 S 的大小。
接下来 mmm 行:每行三个整数 u,v,w,表示编号为 u,v 的点之间有一条权值为 w 的无向边。
接下来一行:k 个互不相同的正整数,表示 S 的元素。
输出格式
第一行:一个整数,表示 E′ 中边权和的最小值。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 #include <bits/stdc++.h> using namespace std ;const int MAXN=510 ;int n,m,k,x,y,z,eg,p[MAXN],hd[MAXN],ver[2 *MAXN],vis[MAXN],nx[2 *MAXN],edge[2 *MAXN],dp[MAXN][4200 ];priority_queue < pair<int ,int > > q; void add_edge (int x,int y,int z) { ver[++eg]=y; nx[eg]=hd[x],edge[eg]=z; hd[x]=eg; return ; } void dijkstra (int s) { memset (vis,0 ,sizeof (vis)); while (!q.empty()) { pair <int ,int > a=q.top(); q.pop(); if (vis[a.second]) {continue ;} vis[a.second]=1 ; for (int i=hd[a.second];i;i=nx[i]) { if (dp[ver[i]][s]>dp[a.second][s]+edge[i]) { dp[ver[i]][s]=dp[a.second][s]+edge[i]; q.push(make_pair(-dp[ver[i]][s],ver[i])); } } } return ; } int main () { memset (dp,0x3f ,sizeof (dp)); scanf ("%d%d%d" ,&n,&m,&k); for (int i=1 ;i<=m;i++) { scanf ("%d%d%d" ,&x,&y,&z); add_edge(x,y,z),add_edge(y,x,z); } for (int i=1 ;i<=k;i++) { scanf ("%d" ,&p[i]); dp[p[i]][1 <<(i-1 )]=0 ; } for (int s=1 ;s<(1 <<k);s++) { for (int i=1 ;i<=n;i++) { for (int subs=s&(s-1 );subs;subs=s&(subs-1 )) { dp[i][s]=min (dp[i][s],dp[i][subs]+dp[i][s^subs]); } if (dp[i][s]!=0x3f3f3f3f ) {q.push(make_pair(-dp[i][s],i));} } dijkstra(s); } printf ("%d\n" ,dp[p[1 ]][(1 <<k)-1 ]); return 0 ; }
最低公共祖先 LCA LCA(Least Common Ancestors),即最近公共祖先,是指在有根树中,找出某两个结点u和v最近的公共祖先。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 struct Edge { int t, nex; } e[500010 << 1 ]; int head[500010 ], total;void add (int x, int y) { e[++total].t = y; e[total].nex = head[x]; head[x] = total; } int depth[500001 ], fa[500001 ][22 ], lg[500001 ];void dfs (int now, int fath) { fa[now][0 ] = fath; depth[now] = depth[fath] + 1 ; for (int i = 1 ; i <= lg[depth[now]]; ++i) fa[now][i] = fa[fa[now][i-1 ]][i-1 ]; for (int i = head[now]; i; i = e[i].nex) if (e[i].t != fath) dfs(e[i].t, now); } int LCA (int x, int y) { if (depth[x] < depth[y]) swap(x, y); while (depth[x] > depth[y]) x = fa[x][lg[depth[x]-depth[y]] - 1 ]; if (x == y) return x; for (int k = lg[depth[x]] - 1 ; k >= 0 ; --k) if (fa[x][k] != fa[y][k]) x = fa[x][k], y = fa[y][k]; return fa[x][0 ]; } int main () { int n, m, s; scanf ("%d%d%d" , &n, &m, &s); for (int i = 1 ; i <= n-1 ; ++i) { int x, y; scanf ("%d%d" , &x, &y); add(x, y); add(y, x); } for (int i = 1 ; i <= n; ++i) lg[i] = lg[i-1 ] + (1 << lg[i-1 ] == i); dfs(s, 0 ); for (int i = 1 ; i <= m; ++i) { int x, y; scanf ("%d%d" ,&x, &y); printf ("%d\n" , LCA(x, y)); } return 0 ; }
欧拉序列+RMQ法 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 #include <bits/stdc++.h> using namespace std ;inline int read () ;const int M = 1000016 ;int n,m,s;vector <int > edge[M];int st[M], ed[M],deep[M];int euler[M*2 ],cnt;void dfs (int now,int fa=0 ) { euler[++cnt] = now; st[now] = cnt; deep[now] = deep[fa]+1 ; for (auto nxt:edge[now]) if (nxt!=fa) { dfs(nxt,now); euler[++cnt] = now; } ed[now] = cnt; } class SpraseTable { static int lg[M]; int n; function<int (int ,int )> cmp; vector <vector <int >> table; public : SpraseTable(int *arr, int _n, function<int (int ,int )> _cmp = [](int a,int b){return a<b?a:b;} ) : n(_n), cmp(_cmp) { if (!lg[0 ]) {lg[0 ]=-1 ;for (int i=1 ;i<M;i++)lg[i]=lg[i/2 ]+1 ;} table = vector <vector <int >>(lg[n] + 1 , vector <int >(n + 1 )); for (int i = 1 ; i <= n; i++) table[0 ][i] = arr[i]; for (int i = 1 ; i <= lg[n]; i++) for (int j = 1 ; j <= n; j++) if (j + (1 << i) - 1 <= n) table[i][j] = cmp(table[i-1 ][j], table[i-1 ][j+(1 <<(i-1 ))]); } inline int query (int x, int y) { int t = lg[y - x + 1 ]; return cmp(table[t][x], table[t][y - (1 << t) + 1 ]); } };int SpraseTable::lg[M]; int main (void ) { n=read (),m=read (),s=read (); for (int i=1 ,a,b;i<n;i++) { a=read (),b=read (); edge[a].push_back(b); edge[b].push_back(a); } dfs(s); SpraseTable stable (euler,cnt,[](int a,int b){return st[a]<st[b]?a:b;}) ; while (m--) { int a = st[read ()], b=st[read ()]; if (a>b) swap(a,b); printf ("%d\n" ,stable.query(a,b)); } return 0 ; } inline int read () { int x=0 ,f=1 ;char ch=getchar(); while (ch<'0' ||ch>'9' ) {if (ch=='-' )f=-1 ;ch=getchar();} while (ch>='0' &&ch<='9' ){x=x*10 +ch-'0' ;ch=getchar();} return x*f; }
最短路 Dijkstra算法 $O(elogn)$ 单源最短路。
d数组保存计算结果(源点s到各个点的最短距离),p数组保存最短路。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 struct Dijkstra { struct Edge { int from, to, dist; Edge(int u, int v, int d): from(u), to(v), dist(d) {} }; int n; vector <Edge> edges; vector <int > G[MAXN]; bool done[MAXN]; int d[MAXN], p[MAXN]; Dijkstra(int n) { this ->n = n; _for(i, 0 , n) G[i].clear (); edges.clear (); } void AddEdge (int from, int to, int dist) { edges.push_back(Edge(from, to, dist)); G[from].push_back(edges.size () - 1 ); } struct HeapNode { int d, u; bool operator < (const HeapNode& rhs) const { return d > rhs.d;} }; void dijkstra (int s) { priority_queue<HeapNode> PQ; _for(i, 0 , n) d[i] = INF; d[s] = 0 ; memset (done, 0 , sizeof (done)); PQ.push((HeapNode){0 , s}); while (!PQ.empty()){ HeapNode x = PQ.top(); PQ.pop(); int u = x.u; if (done[u]) continue ; done[u] = true ; _for(i, 0 , G[u].size ()){ Edge & e = edges[G[u][i]]; if (d[e.to] > d[u] + e.dist) { d[e.to] = d[u] + e.dist; p[e.to] = G[u][i]; PQ.push((HeapNode){d[e.to], e.to}); } } } } }; int main () { freopen("OJ.in" , "r" , stdin ); int n, m, s; cin >> n >> m >> s; Dijkstra d (n) ; _for(i, 0 , m){ int from = readint() - 1 , to = readint() - 1 , w = readint(); d.AddEdge(from, to, w); } d.dijkstra(s - 1 ); _for(i, 0 , n) cout << d.d[i] << " \n" [i == n-1 ]; return 0 ; }
Bellman-Ford算法 $O(en)$ 单源最短路。允许负边权(无负环)。可检测负环 。
这里采用Bellman-Ford的队列改进版本,即SPFA算法。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 struct Bellman_ford { struct Edge { int from, to, dist; Edge(int u, int v, int d): from(u), to(v), dist(d) {} }; int n; vector <Edge> edges; vector <int > G[MAXN]; bool inq[MAXN]; int cnt[MAXN]; int d[MAXN], p[MAXN]; Bellman_ford(int n) { this ->n = n; _for(i, 0 , n) G[i].clear (); edges.clear (); } void AddEdge (int from, int to, int dist) { edges.push_back(Edge(from, to, dist)); G[from].push_back(edges.size () - 1 ); } bool bellman_ford (int s) { queue <int > Q; memset (inq, 0 , sizeof (inq)); memset (cnt, 0 , sizeof (cnt)); _for(i, 0 , n) d[i] = INF; d[s] = 0 ; inq[s] = true ; Q.push(s); while (!Q.empty()){ int u = Q.front(); Q.pop(); inq[u] = false ; _for(i, 0 , G[u].size ()){ Edge & e = edges[G[u][i]]; if (d[u] < INF && d[e.to] > d[u] + e.dist) { d[e.to] = d[u] + e.dist; p[e.to] = G[u][i]; if (!inq[e.to]){ Q.push(e.to); inq[e.to] = true ; if (++cnt[e.to] > n) return false ; } } } } return true ; } };
差分约束系统* 差分约束系统 是一种特殊的 $n$ 元一次不等式组,它包含 $n$ 个变量 $x_1,x_2,…,x^n$ 以及 $m$ 个约束条件,每个约束条件是由两个其中的变量做差构成的,形如 $x_i-x_j\leq c_k$,其中 $c_k$ 是任意常数(可以是非负数,也可以是负数)。我们要解决的问题是:求一组解 $x_1=a_1,x_2=a_2,…,x_n=a^n$,使得所有的约束条件得到满足,否则判断出无解。
求解方法:
将每个变元 $x_i$ 视作一个节点 $i$,每个约束条件$x_i-x_j\leq c_k$ 对应一条从节点 $j$ 到 节点 $i$ 的权为 $c_k$的有向边。
设 $d[i] = 0$ 并向每个节点连一条边,运行Bellman-Ford算法求最短路/最长路 。
若检测负环,系统无解;否则,$x_i=d[i]$ 为该差分约束系统的一组解。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 struct Bellman_ford { struct Edge { int from, to, dist; Edge(int u, int v, int d): from(u), to(v), dist(d) {} }; int n; vector <Edge> edges; vector <int > G[MAXN]; bool inq[MAXN]; int cnt[MAXN]; int d[MAXN], p[MAXN]; Bellman_ford(int n) { this ->n = n; _for(i, 0 , n) G[i].clear (); edges.clear (); } void AddEdge (int from, int to, int dist) { edges.push_back(Edge(from, to, dist)); G[from].push_back(edges.size () - 1 ); } bool bellman_ford (int s) { queue <int > Q; memset (inq, 0 , sizeof (inq)); memset (cnt, 0 , sizeof (cnt)); _for(i, 0 , n) d[i] = -INF; d[s] = 0 ; inq[s] = true ; Q.push(s); while (!Q.empty()){ int u = Q.front(); Q.pop(); inq[u] = false ; _for(i, 0 , G[u].size ()){ Edge & e = edges[G[u][i]]; if (d[u] < INF && d[e.to] < d[u] + e.dist) { d[e.to] = d[u] + e.dist; p[e.to] = G[u][i]; if (!inq[e.to]){ Q.push(e.to); inq[e.to] = true ; if (++cnt[e.to] > n) return false ; } } } } return true ; } }; int main () { int n, m; cin >> n >> m; Bellman_ford B (n) ; _rep(i, 1 , n) B.AddEdge(0 , i, 0 ); _for(i, 0 , m){ int a, b, c; cin >> a >> b>> c; B.AddEdge(a, b, -c); } if (! B.bellman_ford(0 )) cout << "NO" <<endl ; else { _rep(i, 1 , n) cout << B.d[i] << " \n" [i == n]; } return 0 ; }
Floyd算法 $O(n^3)$ 全源最短路。允许负边权(无负环 )。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 #define INF 1000000007 #define _for(i,a,b) for(int i =(a); i < (b); ++i) struct Edge { int from, to, dist; Edge(int u, int v, int d): from(u), to(v), dist(d) {} }; int d[i][j]; void Floyd_init (int n, int m) { _for(i, 0 , n) _for(j, 0 , n) d[i][j] = (i == j) ? 0 : INF; int from, to, dist; _for(i, 0 , m){ cin >> from >> to >> dist; d[from][to] = min (d[from][to], dist); } } void Floyd (int n) { _for(k, 0 , n) _for(i, 0 , n) _for(j, 0 , n) if (d[i][j] < INF && d[k][j] < INF) d[i][j] = min (d[i][j], d[i][k] + d[k][j]); }
最小环 $O(n^3)$ 给出一个图,由 $n$ 个节点构成,问其中的边权和最小的环 $(n\geq 3)$ 是多大。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 int val[maxn + 1 ][maxn + 1 ]; inline int floyd (const int &n) { static int dis[maxn + 1 ][maxn + 1 ]; for (int i = 1 ; i <= n; ++i) for (int j = 1 ; j <= n; ++j) dis[i][j] = val[i][j]; int ans = inf; for (int k = 1 ; k <= n; ++k) { for (int i = 1 ; i < k; ++i) for (int j = 1 ; j < i; ++j) ans = std ::min (ans, dis[i][j] + val[i][k] + val[k][j]); for (int i = 1 ; i <= n; ++i) for (int j = 1 ; j <= n; ++j) dis[i][j] = std ::min ( dis[i][j], dis[i][k] + dis[k][j]); } return ans; }
Johnson算法 $O(nelogn)$ 全源最短路。将负权图 转化为正权图。
我们新建一个虚拟节点 (在这里我们就设它的编号为 0)。从这个点向其他所有点连一条边权为 0 的边。
接下来用 Bellman-Ford 算法求出从 0 号点到其他所有点的最短路,记为 $h_i$ 。
假如存在一条从 $u$ 点到 $v$ 点,边权为 $w$ 的边,则我们将该边的边权重新设置为 $w+h_u-h_v$ 。
接下来以每个点为起点,跑 n 轮 Dijkstra 算法即可求出任意两点间的最短路了。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 struct Dijkstra { };struct Bellman_ford { };int d[MAXN][MAXN]; void Johnson () { int n, m, s; cin >> n >> m >> s; Bellman_ford B (n) ; _for(i, 0 , m){ int from = readint() - 1 , to = readint() - 1 , w = readint(); B.AddEdge(from, to, w); } B.bellman_ford(s); Dijkstra D (n) ; _for(i, 0 , m){ int & u = B.edges[i].from, v = B.edges[i].to, dist = B.edges[i].dist; D.AddEdge(u, v, dist + B.d[u] - B.d[v]); } _for(i, 0 , n){ D.dijkstra(i); _for(j, 0 , n) d[i][j] = D.d[j]; } }
参考实现:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 #include <cstring> #include <iostream> #include <queue> #define INF 1e9 using namespace std ;struct edge { int v,w,next; }e[10005 ]; struct node { int dis,id; bool operator <(const node&a)const { return dis>a.dis; } node(int d,int x) { dis=d,id=x; } }; int head[5005 ],vis[5005 ],t[5005 ];int cnt,n,m;long long h[5005 ],dis[5005 ];void addedge (int u,int v,int w) { e[++cnt].v=v; e[cnt].w=w; e[cnt].next=head[u]; head[u]=cnt; } bool spfa (int s) { queue <int > q; memset (h,63 ,sizeof (h)); h[s]=0 ,vis[s]=1 ; q.push(s); while (!q.empty()) { int u=q.front(); q.pop(); vis[u]=0 ; for (int i=head[u];i;i=e[i].next) { int v=e[i].v; if (h[v]>h[u]+e[i].w) { h[v]=h[u]+e[i].w; if (!vis[v]) { vis[v]=1 ; q.push(v); t[v]++; if (t[v]==n)return false ; } } } } return true ; } void dijkstra (int s) { priority_queue<node> q; for (int i=1 ;i<=n;i++) dis[i]=INF; memset (vis,0 ,sizeof (vis)); dis[s]=0 ; q.push(node(0 ,s)); while (!q.empty()) { int u=q.top().id; q.pop(); if (vis[u])continue ; vis[u]=1 ; for (int i=head[u];i;i=e[i].next) { int v=e[i].v; if (dis[v]>dis[u]+e[i].w) { dis[v]=dis[u]+e[i].w; if (!vis[v])q.push(node(dis[v],v)); } } } return ; } int main () { ios::sync_with_stdio(false ); cin >>n>>m; for (int i=1 ;i<=m;i++) { int u,v,w; cin >>u>>v>>w; addedge(u,v,w); } for (int i=1 ;i<=n;i++) addedge(0 ,i,0 ); if (!spfa(0 )) { cout <<-1 <<endl ; return 0 ; } for (int u=1 ;u<=n;u++) for (int i=head[u];i;i=e[i].next) e[i].w+=h[u]-h[e[i].v]; for (int i=1 ;i<=n;i++) { dijkstra(i); long long ans=0 ; for (int j=1 ;j<=n;j++) { if (dis[j]==INF)ans+=j*INF; else ans+=j*(dis[j]+h[j]-h[i]); } cout <<ans<<endl ; } return 0 ; }
仙人掌最短路 $O(qlogn)$ 任意一条边至多只出现在一条简单回路的无向连通图称为仙人掌。给你一个有 n 个点和 m 条边的仙人掌图,和 q 组询问,每次询问两个点 u,v,求两点之间的最短路。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 #pragma GCC target("avx,sse2,sse3,sse4,popcnt" ) #pragma GCC optimize("O2,Ofast,inline,unroll-all-loops,-ffast-math" ) #include <bits/stdc++.h> #define N 10005 #define M 80005 using namespace std ;inline void rd (int &X) { X=0 ;char ch=0 ; while (!isdigit (ch))ch=getchar(); while ( isdigit (ch))X=(X<<3 )+(X<<1 )+(ch^48 ),ch=getchar(); } int n,m,q,cnt=1 ,tot,num;int head[N],dis[N],v[N];struct nd {int nxt,to,v;}e[M];int d[N],Dis[N],vis[M],pre[N],sum[M],c[N],f[N][15 ];#define For(x) for(int y,i=head[x];(y=e[i].to);i=e[i].nxt) inline void add (int x,int y,int v) { e[++cnt]={head[x],y,v};head[x]=cnt; e[++cnt]={head[y],x,v};head[y]=cnt; } void SPFA () { queue <int > q;q.push(1 ); memset (dis,0x3f ,sizeof dis);dis[1 ]=dis[0 ]=0 ; while (q.size ()){ int x=q.front();q.pop();v[x]=0 ; For(x) if (dis[y]>dis[x]+e[i].v){ dis[y]=dis[x]+e[i].v; if (!v[y]) v[y]=1 ,q.push(y); } } } void work (int x,int rt) { if (x==rt) return ; add(rt,x,0 );vis[pre[x]]=vis[pre[x]^1 ]=1 ; sum[num]+=e[pre[x]].v;c[x]=num;work(e[pre[x]^1 ].to,rt); } void dfs (int x) { v[x]=++tot; For(x) if (i^pre[x]^1 ) if (!v[y]) pre[y]=i,Dis[y]=Dis[x]+e[i].v,dfs(y); else if (v[y]<v[x]) { sum[++num]=e[i].v; vis[i]=vis[i^1 ]=1 ; work(x,y); } } void dfs2 (int x) { for (int i=1 ;i<=14 ;i++) f[x][i]=f[f[x][i-1 ]][i-1 ]; For(x) if (!vis[i] and y!=f[x][0 ]) d[y]=d[x]+1 ,f[y][0 ]=x,dfs2(y); } int ask (int x=0 ,int y=0 ) { rd(x);rd(y); if (d[x]<d[y]) swap(x,y); int X=x,Y=y; for (int i=14 ;~i;i--) if (d[f[x][i]]>=d[y]) x=f[x][i]; if (x==y) return dis[X]-dis[Y]; for (int i=14 ;~i;i--) if (f[x][i]!=f[y][i]) x=f[x][i],y=f[y][i]; if (c[x] and c[x]==c[y]){ int now=sum[c[x]],L=abs (Dis[x]-Dis[y]); return min (L,now-L)+dis[X]-dis[x]+dis[Y]-dis[y]; }return dis[X]+dis[Y]-2 *dis[f[x][0 ]]; } signed main () { rd(n);rd(m);rd(q); for (int x,y,v,i=1 ;i<=m;i++) rd(x),rd(y),rd(v),add(x,y,v); SPFA();dfs(1 );dfs2(1 ); while (q--) printf ("%d\n" ,ask()); }
k短路 给定一个有 $n$ 个结点,$m$ 条边的有向图,求从 $s$ 到 $t$ 的所有不同路径中的第 $k$ 短路径的长度。
A*算法求k短路 $O(nklogn)$ 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 #include <algorithm> #include <cstdio> #include <cstring> #include <queue> using namespace std ;const int maxn = 5010 ;const int maxm = 400010 ;const int inf = 2e9 ;int n, m, s, t, k, u, v, ww, H[maxn], cnt[maxn];int cur, h[maxn], nxt[maxm], p[maxm], w[maxm];int cur1, h1[maxn], nxt1[maxm], p1[maxm], w1[maxm];bool tf[maxn];void add_edge (int x, int y, double z) { cur++; nxt[cur] = h[x]; h[x] = cur; p[cur] = y; w[cur] = z; } void add_edge1 (int x, int y, double z) { cur1++; nxt1[cur1] = h1[x]; h1[x] = cur1; p1[cur1] = y; w1[cur1] = z; } struct node { int x, v; bool operator <(node a) const { return v + H[x] > a.v + H[a.x]; } }; priority_queue<node> q; struct node2 { int x, v; bool operator <(node2 a) const { return v > a.v; } } x; priority_queue<node2> Q; int main () { scanf ("%d%d%d%d%d" , &n, &m, &s, &t, &k); while (m--) { scanf ("%d%d%d" , &u, &v, &ww); add_edge(u, v, ww); add_edge1(v, u, ww); } for (int i = 1 ; i <= n; i++) H[i] = inf; Q.push({t, 0 }); while (!Q.empty()) { x = Q.top(); Q.pop(); if (tf[x.x]) continue ; tf[x.x] = true ; H[x.x] = x.v; for (int j = h1[x.x]; j; j = nxt1[j]) Q.push({p1[j], x.v + w1[j]}); } q.push({s, 0 }); while (!q.empty()) { node x = q.top(); q.pop(); cnt[x.x]++; if (x.x == t && cnt[x.x] == k) { printf ("%d\n" , x.v); return 0 ; } if (cnt[x.x] > k) continue ; for (int j = h[x.x]; j; j = nxt[j]) q.push({p[j], x.v + w[j]}); } printf ("-1\n" ); return 0 ; }
可持久化可并堆优化 k 短路算法 $O(nlogm)$ 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 #include <algorithm> #include <cstdio> #include <cstring> #include <queue> using namespace std ;const int maxn = 200010 ;int n, m, s, t, k, x, y, ww, cnt, fa[maxn];struct Edge { int cur, h[maxn], nxt[maxn], p[maxn], w[maxn]; void add_edge (int x, int y, int z) { cur++; nxt[cur] = h[x]; h[x] = cur; p[cur] = y; w[cur] = z; } } e1, e2; int dist[maxn];bool tf[maxn], vis[maxn], ontree[maxn];struct node { int x, v; node* operator =(node a) { x = a.x; v = a.v; return this ; } bool operator <(node a) const { return v > a.v; } } a; priority_queue<node> Q; void dfs (int x) { vis[x] = true ; for (int j = e2.h[x]; j; j = e2.nxt[j]) if (!vis[e2.p[j]]) if (dist[e2.p[j]] == dist[x] + e2.w[j]) fa[e2.p[j]] = x, ontree[j] = true , dfs(e2.p[j]); } struct LeftistTree { int cnt, rt[maxn], lc[maxn * 20 ], rc[maxn * 20 ], dist[maxn * 20 ]; node v[maxn * 20 ]; LeftistTree() { dist[0 ] = -1 ; } int newnode (node w) { cnt++; v[cnt] = w; return cnt; } int merge (int x, int y) { if (!x || !y) return x + y; if (v[x] < v[y]) swap(x, y); int p = ++cnt; lc[p] = lc[x]; v[p] = v[x]; rc[p] = merge(rc[x], y); if (dist[lc[p]] < dist[rc[p]]) swap(lc[p], rc[p]); dist[p] = dist[rc[p]] + 1 ; return p; } } st; void dfs2 (int x) { vis[x] = true ; if (fa[x]) st.rt[x] = st.merge(st.rt[x], st.rt[fa[x]]); for (int j = e2.h[x]; j; j = e2.nxt[j]) if (fa[e2.p[j]] == x && !vis[e2.p[j]]) dfs2(e2.p[j]); } int main () { scanf ("%d%d%d%d%d" , &n, &m, &s, &t, &k); for (int i = 1 ; i <= m; i++) scanf ("%d%d%d" , &x, &y, &ww), e1.add_edge(x, y, ww), e2.add_edge(y, x, ww); Q.push({t, 0 }); while (!Q.empty()) { a = Q.top(); Q.pop(); if (tf[a.x]) continue ; tf[a.x] = true ; dist[a.x] = a.v; for (int j = e2.h[a.x]; j; j = e2.nxt[j]) Q.push({e2.p[j], a.v + e2.w[j]}); } if (k == 1 ) { if (tf[s]) printf ("%d\n" , dist[s]); else printf ("-1\n" ); return 0 ; } dfs(t); for (int i = 1 ; i <= n; i++) if (tf[i]) for (int j = e1.h[i]; j; j = e1.nxt[j]) if (!ontree[j]) if (tf[e1.p[j]]) st.rt[i] = st.merge( st.rt[i], st.newnode({e1.p[j], dist[e1.p[j]] + e1.w[j] - dist[i]})); for (int i = 1 ; i <= n; i++) vis[i] = false ; dfs2(t); if (st.rt[s]) Q.push({st.rt[s], dist[s] + st.v[st.rt[s]].v}); while (!Q.empty()) { a = Q.top(); Q.pop(); cnt++; if (cnt == k - 1 ) { printf ("%d\n" , a.v); return 0 ; } if (st.lc[a.x]) Q.push({st.lc[a.x], a.v - st.v[a.x].v + st.v[st.lc[a.x]].v}); if (st.rc[a.x]) Q.push({st.rc[a.x], a.v - st.v[a.x].v + st.v[st.rc[a.x]].v}); x = st.rt[st.v[a.x].x]; if (x) Q.push({x, a.v + st.v[x].v}); } printf ("-1\n" ); return 0 ; }
连通分量 强连通分量 有向图 G 强连通是指,G 中任意两个结点连通。
Tarjan 算法 $O(n+m)$ 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 int dfn[N], low[N], dfncnt, s[N], in_stack[N], tp;int scc[N], sc; int sz[N]; void tarjan (int u) { low[u] = dfn[u] = ++dfncnt, s[++tp] = u, in_stack[u] = 1 ; for (int i = h[u]; i; i = e[i].nex) { const int &v = e[i].t; if (!dfn[v]) { tarjan(v); low[u] = min (low[u], low[v]); } else if (in_stack[v]) { low[u] = min (low[u], dfn[v]); } } if (dfn[u] == low[u]) { ++sc; while (s[tp] != u) { scc[s[tp]] = sc; sz[sc]++; in_stack[s[tp]] = 0 ; --tp; } scc[s[tp]] = sc; sz[sc]++; in_stack[s[tp]] = 0 ; --tp; } }
Kosaraju 算法 $O(n+m)$ Kosaraju 算法依靠两次简单的 DFS 实现。
第一次 DFS,选取任意顶点作为起点,遍历所有未访问过的顶点,并在回溯之前给顶点编号,也就是后序遍历。
第二次 DFS,对于反向后的图,以标号最大的顶点作为起点开始 DFS。这样遍历到的顶点集合就是一个强连通分量。对于所有未访问过的结点,选取标号最大的,重复上述过程。
两次 DFS 结束后,强连通分量就找出来了,Kosaraju 算法的时间复杂度为 $O(n+m)$。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 void dfs1 (int u) { vis[u] = true ; for (int v : g[u]) if (!vis[v]) dfs1(v); s.push_back(u); } void dfs2 (int u) { color[u] = sccCnt; for (int v : g2[u]) if (!color[v]) dfs2(v); } void kosaraju () { sccCnt = 0 ; for (int i = 1 ; i <= n; ++i) if (!vis[i]) dfs1(i); for (int i = n; i >= 1 ; --i) if (!color[s[i]]) { ++sccCnt; dfs2(s[i]); } }
双连通分量 在无向连通图中,如果将其中一个点以及所有连接该点的边去掉,图就不再连通,那么这个点就叫做割点(cut vertex / articulation point)。若不存在割点,称图双连通。
割点
输入格式
第一行输入两个正整数 n,m。
下面 m 行每行输入两个正整数 x,y 表示 x 到 y 有一条边。
输出格式
第一行输出割点个数。
第二行按照节点编号从小到大输出割点,用空格隔开。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 #include <bits/stdc++.h> using namespace std ;int n, m; int num[100001 ], low[100001 ], inde, res;bool vis[100001 ], flag[100001 ]; vector <int > edge[100001 ]; void Tarjan (int u, int father) { vis[u] = true ; low[u] = num[u] = ++inde; int child = 0 ; for (auto v : edge[u]) { if (!vis[v]) { child++; Tarjan(v, u); low[u] = min (low[u], low[v]); if (father != u && low[v] >= num[u] && !flag [u]) { flag[u] = true ; res++; } } else if (v != father) low[u] = min (low[u], num[v]); } if (father == u && child >= 2 && !flag[u]) { flag[u] = true ; res++; } } int main () { cin >> n >> m; for (int i = 1 ; i <= m; i++) { int x, y; cin >> x >> y; edge[x].push_back(y); edge[y].push_back(x); } for (int i = 1 ; i <= n; i++) if (!vis[i]) { inde = 0 ; Tarjan(i, i); } cout << res << endl ; for (int i = 1 ; i <= n; i++) if (flag[i]) cout << i << " " ; return 0 ; }
桥/割边 下面代码实现了求割边,其中,当 isbridge[x]
为真时, (father[x],x)
为一条割边。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 int low[MAXN], dfn[MAXN], iscut[MAXN], dfs_clock;bool isbridge[MAXN];vector <int > G[MAXN];int cnt_bridge;int father[MAXN];void tarjan (int u, int fa) { father[u] = fa; low[u] = dfn[u] = ++dfs_clock; for (int i = 0 ; i < G[u].size (); i++) { int v = G[u][i]; if (!dfn[v]) { tarjan(v, u); low[u] = min (low[u], low[v]); if (low[v] > dfn[u]) { isbridge[v] = true ; ++cnt_bridge; } } else if (dfn[v] < dfn[u] && v != fa) { low[u] = min (low[u], dfn[v]); } } }
拓扑排序 $O(e)$
图必须满足是DAG。
核心思想:不断获取零入度的点。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 struct TopologySort { vector <int > G[MAXN]; int inDegree[MAXN]; bool sortable; vector <int > topology; TopologySort(): sortable(false ) {} void input_graph (int n, int m) { int from, to; while (m--){ cin >> from >> to; G[from - 1 ].push_back(to - 1 ); inDegree[to - 1 ]++; } } bool topologicalSort (int n) { queue <int > Q; topology.clear (); _for(i, 0 , n) if (inDegree[i] == 0 ) Q.push(i); int num = 0 ; while (!Q.empty()){ int u = Q.front(); Q.pop(); num++; topology.push_back(u); _for(i, 0 , G[u].size ()){ int v = G[u][i]; inDegree[v]--; if (inDegree[v] == 0 ) Q.push(v); } } return sortable = (num == n); } }; TopologySort t; int main () { int n, m; cin >> n >> m; t.input_graph(n, m); t.topologicalSort(n); _for(i, 0 , n) cout << t.topology[i]+1 <<" \n" [i == n - 1 ]; return 0 ; }
若入度为0的点唯一 ,还可以确定“冠军 ”(每条边看作一场比赛)。
关键路径 计算每个活动的最早开始时间和最晚开始时间。两者相等的活动称为关键活动。
注意:权值都放在了边上。
保证单源点/单汇点,如不满足则增加超级源点/超级汇点,分别连接所有源点/汇点。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 struct Edge { int from, to, dist; Edge(int u, int v, int d): from(u), to(v), dist(d) {} }; vector <Edge> edges; vector <int > G[MAXN]; int inDegree[MAXN]; void AddEdge (int from, int to, int dist) { edges.push_back(Edge(from, to, dist)); G[from].push_back(edges.size () - 1 ); inDegree[to] ++; } ll earliest[MAXN], latest[MAXN]; void criticalPath (int n, int m) { vector <int > topology; queue <int > Q; _for(i, 0 , n) if (inDegree[i] == 0 ){ Q.push(i); earliest[i] = 0 ; } while (!Q.empty()){ int u = Q.front(); Q.pop(); topology.push_back(u); _for(i, 0 , G[u].size ()){ int v = edges[G[u][i]].to, d = edges[G[u][i]].dist; earliest[v] = max (earliest[v], earliest[u] + d); inDegree[v] --; if (inDegree[v] == 0 ) Q.push(v); } } for (int i = topology.size () - 1 ; i >= 0 ; --i){ int u = topology[i]; if (G[u].size () == 0 ){ latest[u] = earliest[u]; } else { latest[u] = INF; } _for(j, 0 , G[u].size ()){ int v = edges[G[u][j]].to, d = edges[G[u][j]].dist; latest[u] = min (latest[u], latest[v] - d); } } } int main () { int n, m; cin >> n >> m; vector <int > a; _for(i, 0 , n) a.push_back(readint()); _for(i, 0 , m){ int from = readint() - 1 , to = readint() - 1 ; AddEdge(from, to, a[to]); } _for(i, 0 , n){ AddEdge(n, i, a[i]); AddEdge(i, n + 1 , 0 ); } bool inDegree_0[MAXN]; _for(i, 0 , n) inDegree_0[i] = (inDegree[i] == 1 ); criticalPath(n + 2 , m); _for(i, 0 , G[n].size ()) if (inDegree_0[i]){ earliest[edges[G[n][i]].to] -= a[i]; latest[edges[G[n][i]].to] -= a[i]; } ll ret = 1 ; _for(i, 0 , n) cout << earliest[i] << " " << latest[i] << endl ; _for(i, 0 , n){ ret *= (latest[i] - earliest[i] + 1 ); ret %= MOD; } ll max_latest = 0 ; _for(i, 0 , n) max_latest = max (max_latest, latest[i]); cout << max_latest << endl ; cout << ret << endl ; return 0 ; }
欧拉回路/欧拉通路 计算一个无向图的欧拉回路或欧拉通路。
给定一张有 500 个顶点的无向图,求这张图的一条欧拉路或欧拉回路。如果有多组解,输出最小的那一组。
在本题中,欧拉路或欧拉回路不需要经过所有顶点。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 #include <algorithm> #include <cstdio> #include <stack> #include <vector> using namespace std ;struct edge { int to; bool exists ; int revref; bool operator <(const edge& b) const { return to < b.to; } }; vector <edge> beg[505 ];int cnt[505 ];const int dn = 500 ;stack <int > ans;void Hierholzer (int x) { for (int & i = cnt[x]; i < (int )beg[x].size ();) { if (beg[x][i].exists ) { edge e = beg[x][i]; beg[x][i].exists = 0 ; beg[e.to][e.revref].exists = 0 ; ++i; Hierholzer(e.to); } else { ++i; } } ans.push(x); } int deg[505 ];int reftop[505 ];int main () { for (int i = 1 ; i <= dn; ++i) { beg[i].reserve(1050 ); } int m; scanf ("%d" , &m); for (int i = 1 ; i <= m; ++i) { int a, b; scanf ("%d%d" , &a, &b); beg[a].push_back((edge){b, 1 , 0 }); beg[b].push_back((edge){a, 1 , 0 }); ++deg[a]; ++deg[b]; } for (int i = 1 ; i <= dn; ++i) { if (!beg[i].empty()) { sort(beg[i].begin (), beg[i].end ()); } } for (int i = 1 ; i <= dn; ++i) { for (int j = 0 ; j < (int )beg[i].size (); ++j) { beg[i][j].revref = reftop[beg[i][j].to]++; } } int bv = 0 ; for (int i = 1 ; i <= dn; ++i) { if (!deg[bv] && deg[i]) { bv = i; } else if (!(deg[bv] & 1 ) && (deg[i] & 1 )) { bv = i; } } Hierholzer(bv); while (!ans.empty()) { printf ("%d\n" , ans.top()); ans.pop(); } }
网络流 最大流
若节点也有流量限制,可拆点(A
拆为A1 --> A2
,流量限制转移到边上),从而转化为纯最大流。
Edmond-Karp算法 $O(nm^2)$ 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 struct Edge { int from, to, cap, flow; Edge(int u, int v, int c, int f = 0 ): from(u), to(v), cap(c), flow(f) {} }; vector <Edge> edges;vector <int > G[MAXN];int a[MAXN]; int p[MAXN]; struct EdmondKarp { int n, m; EdmondKarp(int N){ n = N; edges.clear (); _for(i, 0 , n) G[i].clear (); } void AddEdge (int from, int to, int cap) { edges.push_back(Edge(from, to, cap)); edges.push_back(Edge(to, from, 0 )); m = edges.size (); G[from].push_back(m - 2 ); G[ to].push_back(m - 1 ); } int maxFlow (int s, int t) { int flow = 0 ; while (1 ){ memset (a, 0 , sizeof (a)); queue <int > Q; Q.push(s); a[s] = INF; while (!Q.empty()){ int x = Q.front(); Q.pop(); _for(i, 0 , G[x].size ()){ Edge& e = edges[G[x][i]]; if (!a[e.to] && e.cap > e.flow){ p[e.to] = G[x][i]; a[e.to] = min (a[x], e.cap - e.flow); Q.push(e.to); } } if (a[t]) break ; } if (!a[t]) break ; for (int u = t; u != s; u = edges[p[u]].from){ edges[p[u]].flow += a[t]; edges[p[u]^1 ].flow -= a[t]; } flow += a[t]; } return flow; } }; int main () { int n, m; cin >> n >> m; int s, t; cin >> s >> t; EdmondKarp e (n) ; _for(i, 0 , m){ int from, to, cap; cin >> from >> to >> cap; e.AddEdge(from, to, cap); } cout << e.maxFlow(s, t)<<endl ; return 0 ; }
Dinic算法 $O(n^2m)$ 所有容量为1时Dinic算法复杂度为$O(\text{min}(n^{2/3}, \sqrt m)\times m)$,对于二分图,则为$O(\sqrt n m)$。
最大流 = 最小割,若每条边容量均为1,则等价于返回割边数量。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 inline int read () { int x=0 ,f=1 ;char ch=getchar(); while (!isdigit (ch)){if (ch=='-' ) f=-1 ;ch=getchar();} while (isdigit (ch)){x=x*10 +ch-48 ;ch=getchar();} return x*f; } struct Edge { int from, to, cap, flow; Edge(int u, int v, int c, int f = 0 ): from(u), to(v), cap(c), flow(f) {} }; vector <Edge> edges;vector <int > G[MAXN];bool vis[MAXN]; int d[MAXN]; int a[MAXN]; int cur[MAXN]; struct Dinic { int n, m, s, t; Dinic(int N){ n = N; edges.clear (); _for(i, 0 , n) G[i].clear (); } int Maxflow (int S, int T) { int flow = 0 ; s = S; t = T; while (BFS()){ memset (cur, 0 , sizeof (cur)); flow += DFS(s, INF); } return flow; } void AddEdge (int from, int to, int cap) { edges.push_back(Edge(from, to, cap)); edges.push_back(Edge(to, from, 0 )); m = edges.size (); G[from].push_back(m - 2 ); G[ to].push_back(m - 1 ); } bool BFS () { memset (vis, 0 , sizeof (vis)); queue <int > Q; Q.push(s); d[s] = 0 ; vis[s] = 1 ; while (!Q.empty()){ int x = Q.front(); Q.pop(); _for(i, 0 , G[x].size ()){ Edge& e = edges[G[x][i]]; if (!vis[e.to] && e.cap > e.flow){ vis[e.to] = 1 ; d[e.to] = d[x] + 1 ; Q.push(e.to); } } } return vis[t]; } int DFS (int x, int a) { if (x == t || a == 0 ) return a; int flow = 0 , f; for (int & i = cur[x]; i < G[x].size (); ++i){ Edge& e = edges[G[x][i]]; if (d[x] + 1 == d[e.to] && (f = DFS(e.to, min (a, e.cap-e.flow))) > 0 ){ e.flow += f; edges[G[x][i] ^ 1 ].flow -= f; flow += f; a -= f; if (a == 0 ) break ; } } return flow; } }; int main () { int n = read (), m = read (); int s = read (), t = read (); Dinic dinic (n) ; _for(i, 0 , m){ int from = read (), to = read (), cap = read (); dinic.AddEdge(from, to, cap); } cout << dinic.Maxflow(s, t); return 0 ; }
预流推进算法 HLPP $O(n\sqrt m)$ 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 #include <bits/stdc++.h> #define up(l,r,i) for(register int i=l;i<=r;++i) #define ergv(u) for(std::vector<edge>::iterator p=head[u].begin();p!=head[u].end();++p) #define ergl(u) for(std::list<int>::iterator p=lst[u].begin();p!=lst[u].end();++p) const int INF =2147483647 ;const int MAXN = 1200 +3 ,MAXM =120000 +3 ;int n,m,s,t,mxflow,hst,nwh;std ::vector <int > lft,gap,ht,q,mlst[MAXN];struct edge { int ver,nxt,flw; edge(int _ver,int _flw,int _nxt): ver(_ver),flw(_flw),nxt(_nxt){} }; std ::vector <edge> head[MAXN];std ::list <int > lst[MAXN];std ::vector <std ::list <int >::iterator> it;const int SIZ =100000 +3 ;char buf1[SIZ];char *p1=buf1,*p2=buf1;inline char readchar () { if (p1==p2)p1=buf1,p2=buf1+fread(buf1,1 ,SIZ,stdin ); return p1==p2?EOF:*p1++; } inline int qread () { int ret,c; while ((c=readchar())> '9' ||c< '0' );ret=c-'0' ; while ((c=readchar())>='0' &&c<='9' )ret=ret*10 +c-'0' ; return ret; } inline void add () { int u=qread(),v=qread(),k=qread(); head[u].push_back(edge(v,k,head[v].size () )); head[v].push_back(edge(u,0 ,head[u].size ()-1 )); } inline void rlb () { ht.assign(n+3 ,n),gap.assign(n+3 ,0 ),ht[t]=0 ; q.clear (),q.resize(n+3 ); int front=0 ,rear=0 ,u; for (q[rear++]=t;front<rear;){ u=q[front++]; ergv(u) if (ht[p->ver]==n&&head[p->ver][p->nxt].flw) ++gap[ht[p->ver]=ht[u]+1 ],q[rear++]=p->ver; } up(1 ,n,i) mlst[i].clear (),lst[i].clear (); up(1 ,n,i) if (ht[i]<n){ it[i]=lst[ht[i]].insert(lst[ht[i]].begin (),i); if (lft[i]>0 ) mlst[ht[i]].push_back(i); } hst=(nwh=ht[q[rear-1 ]]); } inline void psh (int u,edge &e) { int v=e.ver,df=std ::min (lft[u],e.flw); e.flw-=df,head[v][e.nxt].flw+=df,lft[u]-=df,lft[v]+=df; if (lft[v]>0 &&lft[v]<=df) mlst[ht[v]].push_back(v); } inline void psh (int u) { int nh=n,htu=ht[u]; ergv(u)if (p->flw) if (ht[u]==ht[p->ver]+1 ){ psh(u,*p); if (lft[u]<=0 ) return ; } else nh=std ::min (nh,ht[p->ver]+1 ); if (gap[htu]==1 ){ up(htu,hst,i){ ergl(i) --gap[ht[*p]],ht[*p]=n; lst[i].clear (); } hst=htu-1 ; }else { --gap[htu],it[u]=lst[htu].erase(it[u]),ht[u]=nh; if (nh==n) return ; gap[nh]++,it[u]=lst[nh].insert(lst[nh].begin (),u); hst=std ::max (hst,nwh=nh),mlst[nh].push_back(u); } } inline int HLPP () { nwh=hst=0 ,ht.assign(n+3 ,0 ),ht[s]=n,it.resize(n+3 ); up(1 ,n,i) if (i!=s)it[i]=lst[ht[i]].insert(lst[ht[i]].begin (),i); gap.assign(n+3 ,0 ),gap[0 ]=n-1 ,lft.assign(n+3 ,0 ),lft[s]=INF,lft[t]=-INF; ergv(s) psh(s,*p);rlb(); for (int u;nwh;){ if (mlst[nwh].empty()) nwh--; else u=mlst[nwh].back(),mlst[nwh].pop_back(),psh(u); } return lft[t]+INF; } int main () { n=qread(),m=qread(),s=qread(),t=qread(); up(1 ,m,i) add(); std ::printf ("%d\n" ,HLPP()); return 0 ; }
全局最小割 Stoer-Wagner 算法 定义无向图 G 的最小割为:一个去掉后可以使 G 变成两个连通分量,且边权和最小的边集的边权和。
给出一张无向图 G,求其最小割。
两个顶点s、t,要么在全局最小割的同一个集合中,要么在不同的集合中。因此分情况讨论即可:
step1:在图G 中找出任意s −t 最小割cut-of-the-phase step2:合并s 、t ,重复执行step1直到图G只剩下一个顶点 step3:输出最小的cut-of-the-phase 为最终结果
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 #include <bits/stdc++.h> using namespace std ;const int MAXN=610 ,INF=1e9 ;int n,m,x,y,z,s,t,dis[MAXN][MAXN],w[MAXN],dap[MAXN],vis[MAXN],ord[MAXN];int proc (int x) { memset (vis,0 ,sizeof (vis)); memset (w,0 ,sizeof (w)); w[0 ]=-1 ; for (int i=1 ;i<=n-x+1 ;i++) { int mx=0 ; for (int j=1 ;j<=n;j++) { if (!dap[j]&&!vis[j]&&w[j]>w[mx]) {mx=j;} } vis[mx]=1 ,ord[i]=mx; for (int j=1 ;j<=n;j++) { if (!dap[j]&&!vis[j]) {w[j]+=dis[mx][j];} } } s=ord[n-x],t=ord[n-x+1 ]; return w[t]; } int sw () { int res=INF; for (int i=1 ;i<n;i++) { res=min (res,proc(i)); dap[t]=1 ; for (int j=1 ;j<=n;j++) { dis[s][j]+=dis[t][j]; dis[j][s]+=dis[j][t]; } } return res; } int main () { scanf ("%d%d" ,&n,&m); for (int i=1 ;i<=m;i++) { scanf ("%d%d%d" ,&x,&y,&z); dis[x][y]+=z,dis[y][x]+=z; } printf ("%d\n" ,sw()); return 0 ; }
最小费用最大流 SPFA法* 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 inline int read () { int x=0 ,f=1 ;char ch=getchar(); while (!isdigit (ch)){if (ch=='-' ) f=-1 ;ch=getchar();} while (isdigit (ch)){x=x*10 +ch-48 ;ch=getchar();} return x*f; } struct Edge { int from, to, cap, flow, cost; Edge(int u, int v, int c, int f, int w): from(u), to(v), cap(c), flow(f), cost(w) {} }; vector <Edge> edges;vector <int > G[MAXN];bool inq[MAXN]; int d[MAXN]; int a[MAXN]; int p[MAXN]; struct MCMF { int n, m; MCMF(int N){ n = N; edges.clear (); _for(i, 0 , n) G[i].clear (); } int MincostMaxflow (int s, int t, long long & cost) { int flow = 0 ; cost = 0 ; while (spfa(s, t, flow, cost)); return flow; } void AddEdge (int from, int to, int cap, int cost) { edges.push_back(Edge(from, to, cap, 0 , cost)); edges.push_back(Edge(to, from, 0 , 0 ,-cost)); m = edges.size (); G[from].push_back(m - 2 ); G[ to].push_back(m - 1 ); } bool spfa (int s, int t, int & flow, long long &