茶茶白的C++模板

在此留下C++模板貌似是个不错的选择。

虽说程序需求千变万化,但是一些短小精湛的函数块还是很值得整理收藏的。
在此提醒自己:找到一种慢的解法可能是找到最终解法的一步。

优质的代码本身就是一种解释,添加完全不必要的注释只会让人恶心

以前觉得 main 函数 return 0 只是标准写法,现在(2020/5/22)才知道能返回就能提前优雅的结束!

全局变量数组元素自动默认初始化为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
83
84
85
86
//#pragma comment(linker,"/STACK:10240000,10240000") // C++ only
/*
#include <cstring>
#include <cstdio>
#include <cstdlib>
#include <ctime>
#include <cmath>
#include <iostream>
#include <algorithm>
#include <vector>
#include <map>
#include <set>
#include <queue>
#include <stack>
#include <string>
#include <bitset>
__builtin_parity(uint n): 返回n的二进制1的个数奇偶性
__builtin_clz(uint n): 返回n的二进制前置0的个数
__builtin_ctz(uint n): 返回n的二进制后置0的个数
__builtin_ffs(int n): 返回n的二进制从后往前第一次出现1的位置
= __builtin_ctz(uint n)+1
__builtin_popcount(uint n): 返回n的二进制1的个数,以上函数仅在GCC中有
lowbit(uint n): n&(-n) 返回使得最大的2^i|n
*/
#if ( ( _WIN32 || __WIN32__ ) && __cplusplus < 201103L)
#define lld "%I64d"
#else
#define lld "%lld"
#endif
#include <bits/stdc++.h>
using namespace std;
using LL = long long;
//typedef long long LL;
//template<typename T>
//typedef __int128 ll; // for g++
//const int INF = 0x3f3f3f3f;
//1e9+7,1e9+9,1e18+3,1e18+9 are prime
//479<<21|1=1004535809,17<<27|1=2281701377 and g=3
/*------ Welcome to my blog: http://dna049.com ------*/
int main(){
//freopen("in","r",stdin);
std::ios::sync_with_stdio(false);std::cin.tie(nullptr);
time_t now = time(0);

cout<<"Time used: "<<difftime(time(0),now)<<endl;
return 0;
}
/*
_ooOoo_
o8888888o
88" . "88
(| -_- |)
O\ = /O
____/`---'\____
.' \\| |// `.
/ \\||| : |||// \
/ _||||| -:- |||||- \
| | \\\ - /// | |
| \_| ''\---/'' | |
\ .-\__ `-` ___/-. /
___`. .' /--.--\ `. . __
."" '< `.___\_<|>_/___.' >'"".
| | : `- \`.;`\ _ /`;.`/ - ` : | |
\ \ `-. \_ __\ /__ _/ .-` / /
======`-.____`-.___\_____/___.-`____.-'======
`=---='
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
佛祖保佑 永无BUG
*/

/**
* ┌───┐ ┌───┬───┬───┬───┐ ┌───┬───┬───┬───┐ ┌───┬───┬───┬───┐ ┌───┬───┬───┐
* │Esc│ │ F1│ F2│ F3│ F4│ │ F5│ F6│ F7│ F8│ │ F9│F10│F11│F12│ │P/S│S L│P/B│ ┌┐ ┌┐ ┌┐
* └───┘ └───┴───┴───┴───┘ └───┴───┴───┴───┘ └───┴───┴───┴───┘ └───┴───┴───┘ └┘ └┘ └┘
* ┌───┬───┬───┬───┬───┬───┬───┬───┬───┬───┬───┬───┬───┬───────┐ ┌───┬───┬───┐ ┌───┬───┬───┬───┐
* │~ `│! 1│@ 2│# 3│$ 4│% 5│^ 6│& 7│* 8│( 9│) 0│_ -│+ =│ BacSp │ │Ins│Hom│PUp│ │N L│ / │ * │ - │
* ├───┴─┬─┴─┬─┴─┬─┴─┬─┴─┬─┴─┬─┴─┬─┴─┬─┴─┬─┴─┬─┴─┬─┴─┬─┴─┬─────┤ ├───┼───┼───┤ ├───┼───┼───┼───┤
* │ Tab │ Q │ W │ E │ R │ T │ Y │ U │ I │ O │ P │{ [│} ]│ | \ │ │Del│End│PDn│ │ 7 │ 8 │ 9 │ │
* ├─────┴┬──┴┬──┴┬──┴┬──┴┬──┴┬──┴┬──┴┬──┴┬──┴┬──┴┬──┴┬──┴─────┤ └───┴───┴───┘ ├───┼───┼───┤ + │
* │ Caps │ A │ S │ D │ F │ G │ H │ J │ K │ L │: ;│" '│ Enter │ │ 4 │ 5 │ 6 │ │
* ├──────┴─┬─┴─┬─┴─┬─┴─┬─┴─┬─┴─┬─┴─┬─┴─┬─┴─┬─┴─┬─┴─┬─┴────────┤ ┌───┐ ├───┼───┼───┼───┤
* │ Shift │ Z │ X │ C │ V │ B │ N │ M │< ,│> .│? /│ Shift │ │ ↑ │ │ 1 │ 2 │ 3 │ │
* ├─────┬──┴─┬─┴──┬┴───┴───┴───┴───┴───┴──┬┴───┼───┴┬────┬────┤ ┌───┼───┼───┐ ├───┴───┼───┤ E││
* │ Ctrl│ │Alt │ Space │ Alt│ │ │Ctrl│ │ ← │ ↓ │ → │ │ 0 │ . │←─┘│
* └─────┴────┴────┴───────────────────────┴────┴────┴────┴────┘ └───┴───┴───┘ └───────┴───┴───┘
*/

int128 在Linux下使用:

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
#include<iostream>
using namespace std;
inline __int128 read(){
__int128 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;
}

inline void print(__int128 x){
if(x<0){
putchar('-');
x=-x;
}
if(x>9) print(x/10);
putchar(x%10+'0');
}

int main(){
__int128 a = read();
print(a*a);
cout<<endl;
return 0;
}

数论篇

向上取整 $\lceil \frac{a}{n} \rceil$

1
2
3
LL ceil(LL a,LL n){  a >= 0, n>0
return (a+n-1)/n; // 注意一般不等于 (a-1)/n+1 !!!
}

注意到 C/C++ 中 ,整数除法 / 是向 0 取整(int(x)也是向0取整),但是Python(Sagemath) 整数除法// 是向下取整。

Greatest Common divisor

1
2
3
4
5
6
7
8
9
10
11
LL gcd(LL a,LL b){
return b?gcd(b,a%b):a;
}
LL exgcd(LL a,LL b,LL& x,LL& y){ // ax + by = gcd(a,b)
if(b==0){
x=1;y=0;return a;
}
LL d=exgcd(b,a%b,y,x);
y-=a/b*x;
return d;
}

当然可以用直接 GNU 内建函数 __gcd

模乘法逆元

1
2
3
4
LL inv(LL a,LL p){ // 0<a<p and gcd(a,p)=1
if(a==1) return 1;
return (p-p/a)*inv(p%a,p)%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
LL pow(LL x,int n){ 
LL r=1;for(;n;n>>=1,x*=x) if(n&1) r*=x;return r;
}
LL powmod(LL x,LL n,LL p){
LL r=1;
while(n){
if(n&1) r=r*x%p;
n>>=1; x=x*x%p;
}
return r;
}
LL mulmod(LL x,LL n,LL p){
LL r=0;
x%=p; n%=p;
while(n){
if(n&1){
r+=x;
if(r>=m) r-=m;
}
n>>=1;x<<=1;
if(x>=m) x-=m;
}
return r;
}
LL mulmod(LL x,LL n,LL p){
LL r=0;
while(n){
if(n&1) r=(r+x)%p;
n>>=1; x=(x+x)%p;
}
return r;
}
LL powmod(LL x,LL n,LL p){
LL r=1;
while(n){
if(n&1) r=mulmod(r,x,p);
n>>=1; x=mulmod(x,x,p);
}
return r;
}

整数开根号

1
2
3
4
5
6
7
8
9
10
int sqrt2(LL x){
LL r = (LL)sqrt(x-0.1);
while(r*r<=x) ++r;
return int(r-1);
}
int sqrt3(LL x){
LL r = (LL)cbrt(x-0.1);
while(r*r*r<=x) ++r;
return int(r-1);
}

自然数方幂和 $O(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
const int N = 1e6+6;
int sp[N],p[N];
LL inv[N],AP[N],AS[N],f[N];
LL getpowsum(LL n,int k,LL mod){ // mod > k
if(k==0) return n%mod;
if(p[0]!=2) spf();
int nk=k+1;
LL tmp = 1;
for(int i=2;i<=nk;++i) tmp=tmp*i%mod;
inv[nk] = powmod(tmp,mod-2,mod);
for(int i=nk-1;i>=0;--i) inv[i]=inv[i+1]*(i+1)%mod;
f[0]=0;f[1]=1;
for(int i=2;i<=nk;++i){
if(sp[i]==i) f[i]=powmod(i,k,mod);
else f[i]=f[sp[i]]*f[i/sp[i]]%mod;
}
for(int i=1;i<=nk;++i){
f[i]+=f[i-1];
if(f[i]>=mod) f[i]-=mod;
}
if(n<=nk) return f[n];
AP[0]=AS[nk]=1;
for(int i=1;i<=nk;++i) AP[i]=AP[i-1]*(n+1-i)%mod;
for(int i=nk-1;i>=0;--i) AS[i]=AS[i+1]*(n-i-1)%mod;
LL res = 0;
for(int i=1;i<=nk;++i){ // because f(i)=0
LL x = f[i]*AP[i]%mod*AS[i]%mod*inv[i]%mod*inv[nk-i]%mod;
if((nk-i)&1) res-=x; // be careful
else res+=x;
}
return (res%mod+mod)%mod;
}

自然数方幂和精确版

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
#include <boost/multiprecision/cpp_int.hpp>
using namespace boost::multiprecision;
cpp_int f[N];
cpp_int getpowsum(LL n,int k){ // k<1000
if(k==0) return cpp_int(n);
if(p[0]!=2) spf();
int nk=2*k+1;
f[0]=0;f[1]=1;
for(int i=2;i<=nk+1;++i){
if(sp[i]==i) f[i]=pow(cpp_int(i),k);
else f[i]=f[sp[i]]*f[i/sp[i]];
}
for(int i=1;i<=nk;++i) f[i]+=f[i-1];
if(n<=nk) return f[n];
cpp_int res = 0,tl=1,tr=1;
for(int i=nk-1;i>=0;--i) tr=tr*(n-i-1)/(nk-i);
for(int i=0;i<=nk;++i){
if((nk-i)&1) res -= f[i]*tl*tr;
else res += f[i]*tl*tr;
tl = tl*(n-i)/(i+1);
tr = tr*(nk-i)/(n-i-1);
}
return res;
}

需要下载 boost包 类似的包还有NTL,GMP

离散对数

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
LL baby_step_giant_step(LL a,LL b,LL p){ // a^x = b mod p
a%=p,b%=p;
if(a==0) return b%p?-1:1;
if(b==1) return 0;
LL cnt = 0,t=1;
for(LL g=__gcd(a,p);g!=1;g=__gcd(a,p)){
if(b%g) return -1;
p/=g,b/=g,t=t*(a/g)%p;
++cnt;
if(b==t) return cnt;
}
map<LL,LL> hash;
LL m = LL(sqrt(p+0.1) + 1);
LL base = b;
for(LL i=0;i!=m;++i){
hash[base] = i;
base = base*a%p;
}
base = powmod(a,m,p);
for(LL i=1;i<=m+1;++i){
t=t*base%p;
if(hash.count(t)) return i*m-hash[t]+cnt;
}
return -1;
}

模素数开根号

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
LL modsqrt(LL a,LL p){  // find x s.t x*x=a mod p;
a = (p+a%p)%p;
if(a==0) return 0;
if(p==2) return (a&1)?1:0;
LL q=(p-1)>>1;
if(powmod(a,q,p)!=1) return -1;
if(q&1) return powmod(a,(q+1)>>1,p);
LL b,cnt=1;
while(powmod(b=rand()%p,q,p)==1);//find a non quadratic residue
while(!(q&1)) ++cnt,q>>=1;
b=powmod(b,q,p);
LL x=powmod(a,(q+1)>>1,p);
for(LL s=1,t=powmod(a,q,p);t!=1;s=1){ //keep x*x=t*a;t^{2^{cnt-1}}=1
for(LL tt=t*t%p;s<cnt&&tt!=1;++s) tt=tt*tt%p;
LL d=powmod(b,1<<(cnt-s-1),p);
x=(x*d)%p;b=d*d%p;t=t*b%p;cnt=s;
}
return 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
LL mod_sqrt_p(LL a,LL n,LL p,int k){//find x s.t x*x=a mod p^k where (a,p)=1,p>2
LL q = (n/p)*((p-1)>>1);
if(powmod(a,q,n)!=1) return -1;
if(q&1) return powmod(a,(q+1)>>1,n);
LL b,cnt=1;
while(powmod(b=rand()%n,q,n)==1);//find a non quadratic residue
while(!(q&1)) ++cnt,q>>=1;
b=powmod(b,q,n);
LL x=powmod(a,(q+1)>>1,n);
for(LL s=1,t=powmod(a,q,n);t!=1;s=1){// keep x*x=t*a;t^{2^{cnt-1}}=1
for(LL tt=t*t%p;s<cnt&&tt!=1;++s) tt=tt*tt%p;
LL d = powmod(b,1<<(cnt-s-1),n);
x=(x*d)%n;b=d*d%n;t=t*b%n;cnt=s
}
return x;
}
LL pow(LL x,int n){
LL r=1;for(;n;n>>=1,x*=x) if(n&1) r*=x;return r;
}
LL mod_sqrt(LL a,LL p,int k=1){//find smallest x>=0 s.t x*x=a mod p^k
if(a==0) return 0;
int ka=0;
while(a%p==0) a/=p,++ka,--k;
if(k<=0) return 0;
if(ka&1) return -1;
LL n = pow(p,k),x;
if(p==2){
if(k==1||k==2) x = a==1?1:-1;
else if(a%8!=1) x=-1;
else{
x=1;
for(int i=4;i<=k;++i){
if( (x*x)%(1<<i) == a%(1<<i)) continue;
x+=1<<(i-2);
}
}
}else x = mod_sqrt_p(a,n,p,k);
return x==-1?-1:pow(p,ka>>1)*(x<n-x?x:n-x);
}

对于模一般的 $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
LL C(int n,int m){
if(n<m) return 0;
LL ans = 1;
for(int i=n;i+m>n;--i){
ans = ans*i/(n-i+1);
}
return ans;
}
const int N = 12;
void get_C(){
C[0][0]=C[1][0]=C[1][1]=1;
for(int i=2;i<=N;++i){
c[i][0]=c[i][i]=1;
for(int j=1;j<i;++j){
c[i][j]=c[i-1][j]+c[i-1][j-1];
}
}
}
LL f[N]; // can add invf[N] for speed
LL Lucas(LL n, LL m, LL p) { // C(n,m)%p,f[n]=n!%p
LL r = 1;
while(n&&m){
LL np=n%p,mp=m%p;
if(np<mp) return 0;
r = r*f[np]%p*inv(f[mp]*f[np-mp]%p,p)%p;
n/=p,m/=p;
}
return 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
string rational(int n,int m){
if(m==0) return "";
string r;
if(m<0) n=-n,m=-m;
if(n<0){
n=-n;
r.push_back('-');
}
int t = n/m; // get the integer part of n/m
if(t==0) r.push_back('0');
while(t){
r.push_back(t%10+'0');
t/=10;
}
r.reserve();
n%=m;
if(n==0) return r;
r.push_back('.');
t = __gcd(n,m);
n/=t;m/=t;
while(m%10==0){
m/=10;
r.push_back(n/m+'0');
n%=m;
}
while(m%2==0){
m/=2;n*=5;
r.push_back(n/m+'0');
n%=m;
}
while(m%5==0){
m/=5;n*=2;
r.push_back(n/m+'0');
n%=m;
}
if(m==1) return r;
r.push_back('(');
t=1;
do{
r.push_back(10*n/m+'0');
n=10*n%m;
t=10*t%m;
}while(t!=1);
r.push_back(')');
return r;
}

正整数拆分方案数

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
const int N = 1e5+2;
const LL M = 1e9+7;
int a[N];
void init(){
a[0]=1;
for(int i=1;i<N;++i){
for(int j=1;j*(3*j-1)/2<=i;++j){
(a[i]+=(j&1?1:-1)*a[i-j*(3*j-1)/2])%=M;
}
for(int j=-1;j*(3*j-1)/2<=i;--j){
(a[i]+=(j&1?1:-1)*a[i-j*(3*j-1)/2])%=M;
}
if(a[i]<0) a[i]+=M;
}
}

大素数Miller-Rabin概率判别法

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
bool Witness(LL a,LL n,LL m,int t){
LL x=powmod(a,m,n);
if(x==1||x==n-1) return false;
while(t--){
x=mulmod(x,x,n);
if(x==n-1) return false;
}
return true;
}
bool Rabin(LL n){
if(n<2) return false;
if(n==2) return true;
if(!(n&1)) return false;
LL m=n-1;
int t=0,cnt=33;
while(!(m&1)){
++t;m>>=1;
}
while(cnt--){
LL a=rand()%(n-1)+1;
if(Witness(a,n,m,t)) return false;
}
return true;
}

线性素数筛(可以分段来搞节省空间)

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
const int N = 1e9 + 9;
const int NN = 1e8 + 8;
int p[NN];
bool isp[N];
int initprime(){ // 8.6s 推荐,可以做常数优化但没必要!
int cnt=1;p[1]=2;isp[2] = true;
for(int i=3;i<N;i+=2) isp[i]=true;
for(int i=3;i<N;i+=2){
if(isp[i]) p[++cnt] = i;
for(int j = 2,t=(N-1)/i + 1;j <= cnt && p[j] < t; ++j){ // 用除号是防止溢出
isp[i * p[j]] = false;
// 不要下面的一步的话,复杂度 O(nloglogn), 优势在于:不用除法,常数小
if(i%p[j] == 0) break;
}
}
return cnt;
}
bool isprime(int n){
if(p[1]!=2) initprime();
if(n<N) return isp[n];
for(int i=0;p[i]*p[i]<=n;++i){
if(n%p[i]==0) return false;
}
return true;
}

$O(n)$ 最小素因子预处理

1
2
3
4
5
6
7
8
9
10
11
12
13
14
const int N = 2e8; // 再大内存吃不消了 
int sp[N],p[N];
int spf(){ // samllest prime factor i5-3470 2s 执行完
int cnt=0;p[cnt++]=2;
for(int i=2;i<N;i+=2) sp[i]=2;
for(int i=1;i<N;i+=2) sp[i]=i;
for(int i=3;i<N;i+=2){
if(sp[i]==i) p[cnt++] = i;
for(int j = 1;j < cnt && p[j]<=sp[i] && i * p[j] < N; ++j){ //防止乘法溢出
sp[i * p[j]] = p[j]; // 注意到sp只被赋值一次
}
}
return cnt;
}

大整数的最小素因子

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
LL pollardrho(LL n){
LL x=rand()%(n-1)+1;
LL y=x,i=1,k=2,c=rand()%(n-1)+1;
while(true){
x=(mulmod(x,x,n)+c)%n;
LL d=gcd(y-x+n,n);
if(d>1) return d;
if(++i==k){
y=x;
if(k>n) return n;
k<<=1;
}
}
}
LL ans;
void findp(LL n){
if(Rabin(n)){
ans=min(ans,n);
return;
}
LL p=n;
while(p==n) p=pollardrho(n);
findp(p);
findp(n/p);
}

Mobius function

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
const int N = 1e7+2;
int p[N],pi[N],mu[N];
bool isp[N];
int initmu(){
int cnt=1;p[1]=2;mu[1] = 1;isp[2] = true;
for(int i=3;i<N;i+=2) isp[i]=true;
for(int i = 3;i<N;i+=2){
if(isp[i]) mu[i] = -1,p[++cnt] = i;
for(int j=1,t=(N-1)/i + 1;j<=cnt && p[j]<t;++j){
isp[i * p[j]] = false;
if(i%p[j] == 0) break;
mu[i*p[j]] = -mu[i];
}
}
for(int i=2;i<N;i+=4) mu[i]=-mu[i>>1];
return cnt;
}
int getmu(int n){
int r=1;
for(int i=0;p[i]*p[i]<=n;++i){
if(n%p[i]==0){
n/=p[i];
if(n%p[i]==0) return 0;
r=-r;
}
}
return n>1?-r:r;
}
int sumu[N];
void initsumu(){
if(mu[1]!=1) initmu();
for(int i=1;i<N;++i) sumu[i]=sumu[i-1]+mu[i];
}
map<int,int> mp;
int M(int n){ // M(n) = M(n-1) + mu(n)
if(n<N) return sumu[n];
map<int,int>::iterator it = mp.find(n);
if(it!=mp.end()) return it->second;
int r=1;
for(int i=2,j;i<=n;i=j+1){
j=n/(n/i);
r-=(j-i+1)*M(n/i);
}
mp.insert(pair<int,int>(n,r));
return r;
}
LL Q(LL n){ // Q(n) = Q(n-1) + |mu(n)|
LL r=0;
for(LL i=1,t;(t=i*i)<n;++i){
r+=mu[i]*(n/t);
}
return r;
}

Euler totient function

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
const int N = 1e7+2;
int phi[N],p[N];
void init_phi(){
for(int i=1;i<N;++i) phi[i]=i;
for(int i=2,cnt=0;i<N;++i){
if(phi[i]==i){
--phi[i];
p[cnt++]=i;
}
for(int j=0;j<cnt&&i*p[j]<N;++j){
if(i%p[j]==0){
phi[i*p[j]]=phi[i]*p[j];break;
}
phi[i*p[j]]=phi[i]*(p[j]-1);
}
}
}
void init_phi(){
for(int i=1;i<N;i+=2) phi[i]=i;
for(int i=2;i<N;i+=2) phi[i]=i>>1;
for(int i=3;i<N;i+=2){
if(phi[i]!=i) continue;
for(int j=i;j<N;j+=i) phi[j]=phi[j]/i*(i-1);
}
}
LL getphi(LL x){
LL r=x;
for(int i=0;p[i]*p[i]<=x;++i){
if(x%p[i]==0){
r=r/p[i]*(p[i]-1);
while(x%p[i]==0) x/=p[i];
}
}
if(x>1) r=r/x*(x-1);
return r;
}
LL sumphi[N];
void init_sumphi(){
if(phi[2]!=1) init_phi();
for(int i=1;i<N;++i) sumphi[i]=sumphi[i-1]+phi[i];
}
map<int,LL> mp;
LL getsumphi(int n){
if(n<N) return (LL)sumphi[n];
map<int,LL>::iterator it = mp.find(n);
if(it!=mp.end()) return it->second;
LL r=LL(n+1)*n/2;
for(int i=2,j;i<=n;i=j+1){
j=n/(n/i);
r-=(j-i+1)*getsumphi(n/i);
}
mp.insert(pair<int,LL>(n,r));
return r;
}

$\pi(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
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<bits/stdc++.h>
using namespace std;
typedef long long LL;
const int N = 1e8+2;
int p[N],pi[N];
bool isp[N];
int initprime(){
int cnt=1;p[1]=2;isp[2] = true;
for(int i=3;i<N;i+=2) isp[i]=true;
for(int i=3;i<N;i+=2){
if(isp[i]) p[++cnt] = i;
for(int j = 2,t=(N-1)/i + 1;j <= cnt && p[j] < t; ++j){
isp[i * p[j]] = false;
if(i%p[j] == 0) break;
}
}
return cnt;
}
const int M = 7;
const int PM = 2*3*5*7*11*13*17;
int phi[PM+1][M+1],sz[M+1];
void init(){
initprime();
pi[2] = 1;
for(int i=3;i<N;++i){
pi[i]=pi[i-1];
if(isp[i]) ++pi[i];
}
sz[0]=1;
for(int i=0;i<=PM;++i) phi[i][0]=i;
for(int i=1;i<=M;++i){
sz[i]=p[i]*sz[i-1];
for(int j=1;j<=PM;++j){
phi[j][i]=phi[j][i-1]-phi[j/p[i]][i-1];
}
}
}
LL primepi(LL x);
LL primephi(LL x, int s){
if(s <= M) return phi[x%sz[s]][s]+(x/sz[s])*phi[sz[s]][s];
if(x/p[s] <= p[s]) return primepi(x)-s+1;
if(x/p[s]/p[s] <= p[s] && x<N){
int s2x = pi[(int)(sqrt(x+0.2))];
LL ans = pi[x]-(s2x+s-2)*(s2x-s+1)/2;
for(int i=s+1;i<=s2x;++i){
ans+=pi[x/p[i]];
}
return ans;
}
return primephi(x,s-1)-primephi(x/p[s],s-1);
}
LL primepi(LL x){
if(x<N) return pi[x];
int ps2x = primepi(int(sqrt(x+0.2)));
int ps3x = primepi(int(cbrt(x+0.2)));
LL ans = primephi(x,ps3x) + LL(ps2x+ps3x-2)*(ps2x-ps3x+1)/2;
for(int i =ps3x+1,ed = ps2x;i<=ed;++i){
ans -= primepi(x/p[i]);
}
return ans;
}

int main(){
init();
LL n;
while(cin>>n){ // n=98765432109876 = 9.8*10^13 用时 42s, N 大点,耗时会小点
time_t now = time(0);
cout<<primepi(n)<<endl;
cout<<"Time used: "<<difftime(time(0),now)<<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
using namespace std;
typedef long long LL;
const int N = 1e7+2;
LL L[N],R[N];
LL primepi(LL n){
LL rn = (LL)sqrt(n+0.2);
for(LL i=1;i<=rn;++i) R[i]=n/i-1;
LL ln = n/(rn+1)+1;
for(LL i=1;i<=ln;++i) L[i]=i-1;
for(LL p=2;p<=rn;++p){
if(L[p]==L[p-1]) continue;
for(LL i=1,tn=min(n/(p*p),rn);i<=tn;++i){
R[i] -= (i*p<=rn?R[i*p]:L[n/(i*p)])-L[p-1];
}
for(LL i=ln;i>=p*p;--i){
L[i] -= L[i/p]-L[p-1];
}
}
return R[1];
}
int main(){
LL n;
while(cin>>n){ // n=98765432109876 = 9.8*10^13 用时 193s
cout<<primepi(n)<<endl;
}
return 0;
}

求奇素数的一个原根

首先,模 $m$ 有原根的充要条件是:$m=2,4,p^a,2p^a$,其中$p$为奇素数。
对于求模$p$的原根方法:对 $p-1$ 素因子分解:$p-1 = p_1^{a_1} \cdots p_s^{a_s}$ 若恒有

则 $g$ 是 模$p$的原根。对于 $p^a$ 可以用 $p$ 的原根简单构造,而 $2p^a$ 的原根为 $p^a$ 的原根 与 $p^a$ 的原根和 $p^a$的和中奇数者。(证明见P150《数论基础》潘承洞)

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
LL fact[N];
int factor(LL n){
int cnt = 0;
for(int i=0;LL(p[i])*p[i]<=n; ++i){
if(n%p[i]==0){
fact[cnt++] = p[i];
while(n%p[i]==0) n /= p[i];
}
}
if(n > 1) fact[cnt++] = n;
return cnt;
}
LL proot(LL p){ // p must be odd prime
int cnt = factor(p-1);
for(LL i = 1; i<p;++i){
bool flag = true;
for(int j=0;j<cnt;++j){
if(powmod(i,(p-1)/fact[j],p)==1){
flag = false; break;
}
}
if(flag) return i;
}
return -1;
}

求所有原根见

数论函数的Dirichlet乘积

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
class Dirichlet{
public:
const static int N = 1e5+2;
int a[N],n;
Dirichlet(int _n = 0){
n=_n;
for(int i=1;i<=n;++i){
a[i]=0;
}
}
Dirichlet operator*(const Dirichlet& A)const{
Dirichlet R(n);
for(int i=1;i<=n;++i){
for(int j=1;i*j<=n;++j){
R.a[i*j]+=a[i]*A.a[j];
}
}
return R;
}
};
LL Dirichlet(int n,LL a[],LL b[]){
LL ans = 0;int s;
for(s=1;s*s<n;++s){
if(n%s==0) continue;
ans+=a[s]*b[n/s]+b[s]*a[n/s];
}
if(s*s==n) ans+=a[s]*b[s];
return ans;
}

线性同余方程

1
2
3
4
5
6
7
int modeq(LL a,LL b,LL n){ // a*x=b mod n return x
LL x,y,d;
d=exgcd(a,n,x,y);
if(b%d!=0) return -1;
a/=d;n/=d;b/=d;
return ((x*b)%n+n)%n;
}

中国剩余定理

1
2
3
4
5
6
7
8
9
10
11
12
13
LL chineseRemain(int n,LL a[],LL m[]){ //x=a[i] mod m[i], 本质上是处理n=2的情况
LL x,y,a1,a2,m1,m2,d;
a1=a[0],m1=m[0];
for(int i=1;i<n;++i){
a2=a[i],m2=m[i];
d=exgcd(m1,m2,x,y);
if((a2-a1)%d!=0) return -1;
a1+=(((a2-a1)/d*x)%m2+m2)%m2*m1;
m1=m1/d*m2;
a1=(a1%m1+m1)%m1;
}
return a1;
}

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
void change(double *x,double *y,int len,int loglen){
for(int i=0;i<len;++i){
int t=i,k=0;
for(int j=0;j<loglen;++j,t>>=1){
k=(k<<1)|(t&1);
}
if(k<i){
double tmp=x[i];x[i]=x[k];x[k]=tmp;
tmp=y[i];y[i]=y[k];y[k]=tmp;
}
}
}
void fft(double *x,double *y,int len,int loglen,bool isInverse){
change(x,y,len,loglen);
for(int step=2;step<=len;step<<=1){
int half=step>>1;
double theta=PI/half;
double wmx=cos(theta),wmy=sin(theta);
if(isInverse) wmy=-wmy;
for(int i=0;i<len;i+=step){
double wx=1,wy=0;
for(int j=0;j<half;++j){
double cx=x[i+j],cy=y[i+j];
double dx=x[i+j+half],dy=y[i+j+half];
double ex=dx*wx-dy*wy,ey=dx*wy+dy*wx;
x[i+j]=cx+ex;y[i+j]=cy+ey;
x[i+j+half]=cx-ex;y[i+j+half]=cy-ey;
double tmp=wx*wmx-wy*wmy;
wy=wx*wmy+wy*wmx;wx=tmp;
}
}
}
if(isInverse) for(int i=0;i<len;++i){
x[i]/=len;y[i]/=len;
}
}

利用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
void mul(int *a,int *b,int len,int loglen,bool same){
for(int i=0;i!=len;++i){
ax[i]=a[i];bx[i]=b[i];
ay[i]=by[i]=0;
}
fft(ax,ay,len,loglen,0);
if(!same) fft(bx,by,len,loglen,0);
else{
for(int i=0;i!=len;++i){
bx[i]=ax[i];by[i]=ay[i];
}
}
for(int i=0;i!=len;++i){
double cx=ax[i]*bx[i]-ay[i]*by[i];
double cy=ax[i]*by[i]+ay[i]*bx[i];
ax[i]=cx;ay[i]=cy;
}
fft(ax,ay,len,loglen,1);
for(int i=0;i!=len;++i){
a[i]=(int)(ax[i]+0.5);
}
for(int i=q-1;i>=0;--i){
a[i+q-p]=(a[i+q-p]+ta*a[q+i])%M;
a[i]=(a[i]+tb*a[q+i])%M;
a[i+q]=0;
}
}

快速数论变换 NFT

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
void change(LL *x,int len,int loglen){
for(int i=1;i<len;++i){
int t=i,k=0;
for(int j=0;j<loglen;++j,t>>=1){
k=(k<<1)|(t&1);
}
if(k<i) swap(x[i],x[k]);
}
}
const LL FM = 479<<21|1;
void nft(LL *x,int len,int loglen,bool isInverse){
LL g = powmod(3,(FM-1)>>loglen,FM);
if(isInverse){
g=inv(g,FM);
LL invlen = powmod(len,FM-2,FM);
for(int i=0;i<len;++i){
x[i]=x[i]*invlen%FM;
}
}
change(x,len,loglen);
for(int step=2;step<=len;step<<=1){
int half = step>>1;
LL wn = powmod(g,len/step,FM);
for(int i=0;i<len;i+=step){
LL w = 1;
for(int j = i;j<i+half;++j){
LL t=(w*x[j+half])%FM;
x[j+half]=(x[j]-t+FM)%FM;
x[j]=(x[j]+t)%FM;
w = w*wn%FM;
}
}
}
}
void mul(LL *a,LL *b,int len,int loglen,bool same){
nft(a,len,loglen,0);
if(same){
for(int i=0;i<len;++i) a[i]=a[i]*a[i]%FM;
}else{
nft(b,len,loglen,0);
for(int i=0;i<len;++i) a[i]=a[i]*b[i]%FM;
}
nft(a,len,loglen,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
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
class Node{
public:
int c,d;
Node* next;
Node(){}
Node(int _c,int _d):c(_c),d(_d){next=NULL;}
};
class Poly{ // a0+a1*x+...+an*x^n
friend ostream& operator<<(ostream& os,const Poly& A){
if(A.head->c==-1){
if(A.head->d==0) os<<"-1";
if(A.head->d==1) os<<"-x";
if(A.head->d>1) os<<"-x^"<<(A.head->d);
}else{
os<<A.head->c;
if(A.head->d==1) os<<"x";
if(A.head->d>1) os<<"x^"<<(A.head->d);
}
Node* p=A.head->next;
while(p!=NULL){
if(p->c>0) os<<"+";
if(p->c==-1) os<<"-";
else if(p->c!=1) os<<(p->c);
if(p->d>1) os<<"x^"<<(p->d);
if(p->d==1) os<<"x";
p=p->next;
}
return os;
}
private:
int Size;
Node *head;
public:
Poly(){
head=new Node(0,0);
Size=-1;
}
Poly(int x){
head=new Node(x,0);
if(x) Size=0;
else Size=-1;
}
Poly(int* a,int n){
if(n<1) {Poly();return;}
head=NULL;
Node* p;
for(Size=n-1;Size!=-1&&a[Size]==0;--Size);
for(int i=Size;i!=-1;--i){
if(a[i]){
p=new Node(a[i],i);
p->next=head;
head=p;
}
}
if(head==NULL) head=new Node(0,0);
// if(head==NULL) new (this)poly();
}
Poly(const Poly& A){
head=new Node(A.head->c,A.head->d);
Size=A.Size;
Node *Ap=A.head->next,*p=head,*tp;
while(Ap!=NULL){
tp=new Node(Ap->c,Ap->d);
p->next=tp;p=tp;
Ap=Ap->next;
}
}
~Poly(){
del();
}
void del(){
Node* p;
while(head!=NULL){
p=head->next;
delete(head);
head=p;
}
Size=-1;
}
Poly& operator=(const Poly& A){
del();
Size=A.Size;
head=new Node(A.head->c,A.head->d);
Node *Ap=A.head->next,*p=head,*tp;
while(Ap!=NULL){
tp=new Node(Ap->c,Ap->d);
p->next=tp;p=tp;
Ap=Ap->next;
}
return *this;
}
Poly operator+(const Poly& A)const{
if(Size<A.Size) return add(A,*this);
else return add(*this,A);
}
Poly operator-()const{
Poly R(*this);
Node *Rp=R.head;
while(Rp!=NULL){
Rp->c=-Rp->c;
Rp=Rp->next;
}
return R;
}
Poly operator-(const Poly& A)const{
if(Size<A.Size) return add(-A,*this);
else return add(*this,-A);
}
Poly mul(int tc,int td)const{
if(tc==0) return Poly();
if(Size==-1) return *this;
Poly R(*this);
Node *Rp=R.head;
while(Rp!=NULL){
Rp->c*=tc;
Rp->d+=td;
R.Size=Rp->d;
Rp=Rp->next;
}
return R;
}
Poly operator*(Poly& A)const{
Poly R;
Node *Ap=A.head;
while(Ap!=NULL){
R=R+mul(Ap->c,Ap->d);
Ap=Ap->next;

}
return R;
}
bool find(int td)const{
Node *p=head;
while(p!=NULL&&p->d<td) p=p->next;
if(p==NULL||p->d>td) return false;
else return true;
}
int getans(int td)const{
Node *p=head;
int ans=0;
while(p!=NULL&&p->d<=td){
ans+=p->c;
p=p->next;
}
return ans;
}
private:
friend Poly add(const Poly& A,const Poly& B){ //A.Size>=B.Size
Poly R(A);
if(B.Size==-1) return R;
Node *Bp=B.head,*Rp=R.head,*tp;
if(R.head->d>B.head->d){
tp=new Node(B.head->c,B.head->d);
tp->next=R.head;
Rp=R.head=tp;
Bp=Bp->next;
}
while(Bp!=NULL){
if(Rp->d==Bp->d){
Rp->c+=Bp->c;
Bp=Bp->next;
}else{
if(Rp->next->d>Bp->d){
tp=new Node(Bp->c,Bp->d);
tp->next=Rp->next;
Rp->next=tp;Rp=tp;
Bp=Bp->next;
}
else Rp=Rp->next;
}
}
while(R.head->next!=NULL&&R.head->c==0){
Rp=R.head;
R.head=Rp->next;
delete(Rp);
}
if(R.head->c==0) R.Size=-1;
else R.Size=R.head->d;
Rp=R.head;
while(Rp->next!=NULL){
if(Rp->next->c==0){
tp=Rp->next;
Rp->next=tp->next;
delete(tp);
}
Rp=Rp->next;
}
return 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
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
class BigInt{ // without nft and we won't set B=10000 in nft
const static int N=30;
const static int B=10000;
const static int BB=4;
public:
int len,a[N]; // keep the inverse order of int > 0
void clear(){
len=0;
memset(a,0,sizeof(a));
}
BigInt(){
clear();
}
BigInt(int x){ // construct function use c.f. brings mistake in clang
clear();
while(x>=B){
a[len++] = x%B;
x/=B;
}
a[len++]=x;
}
BigInt(LL x){ // construct function use c.f. brings mistake in clang
clear();
while(x>=B){
a[len++] = x%B;
x/=B;
}
a[len++]=(int)x;
}
BigInt(const char *s){
clear();
int slen=(int)strlen(s);
while(slen>0&&a[slen-1]=='0') --slen;
if(slen==0) a[len++] = 0;
else{
while(slen>0){
for(int i=1;i<B;i*=10){
if(--slen>=0) a[len]=i*s[slen]-'0';
}
++len;
}
}
}
BigInt(const BigInt& A):len(A.len){
memset(a,0,sizeof(a));
for(int i=0;i!=len;++i) a[i]=A.a[i];
}
BigInt operator+(const BigInt& A)const{
BigInt R;
R.len=max(len,A.len);
int res=0;
for(int i=0;i!=R.len;++i){
res+=a[i]+A.a[i];
R.a[i]=res%B;res/=B;
}
if(res) R.a[R.len++]=res%B;
while(R.len>1&&R.a[R.len-1]==0) --R.len;
return R;
}
BigInt operator*(const BigInt& A)const{
BigInt R;
for(int i=0;i!=len;++i){
int res = 0;
for(int j=0;j!=A.len;++j){
res+=a[i]*A.a[j]+R.a[i+j];
R.a[i+j]=res%B;res/=B;
}
R.a[i+A.len]+=res;
}
R.len = len + A.len;
while(R.len>1&&R.a[R.len-1]==0) --R.len;
return R;
}
bool operator<(const BigInt& A)const{
if(len<A.len) return true;
if(len>A.len) return false;
int tn = len;
while(--tn>=0){
if(a[tn]<A.a[tn]) return true;
if(a[tn]>A.a[tn]) return false;
}
return false;
}
bool isOdd()const{
return a[0]&1;
}
BigInt gethalf()const{
BigInt R = *this;
bool tmp = false;
for(int i=R.len-1;i>=0;--i){
if(tmp) R.a[i]+=B;
if(R.a[i]&1) tmp=true;
else tmp = false;
R.a[i]>>=1;
}
if(R.a[R.len-1] == 0) --R.len;
return R;
}
void print()const{
printf("%d",a[len-1]);
char ss[10];
sprintf(ss,"%%0%dd",BB);
for(int i=len-2;i>=0;--i){
printf(ss,a[i]);
}
puts("");
}
};

$\mathbb{Z}_p$ 二次拓域

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
LL p,q;
class Zpq{
public:
LL a,b; // a+b sqrt q,0 <= a,b
Zpq(){a=b=0;}
Zpq(LL _a,LL _b):a(_a%p),b(_b%p){}
Zpq operator+(const Zpq& A)const{
return Zpq(a+A.a,b+A.b);
}
Zpq operator*(const Zpq& A)const{
return Zpq(a*A.a+b*A.b%p*q,a*A.b+b*A.a);
}
};
Zpq pow(Zpq A,LL n){
Zpq R(1,0);
while(n){
if(n&1) R=R*A;
n>>=1; A=A*A;
}
return 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
LL mod;
class Matrix{
public:
const static int N = 3;
LL a[N][N];
Matrix(){
all(0);
}
Matrix(int x){ // xIn
all(0);
for(int i=0;i<N;++i){
a[i][i]=x;
}
}
void all(int x){
for(int i=0;i<N;++i){
for(int j=0;j<N;++j){
a[i][j]=x;
}
}
}
Matrix operator+(const Matrix& A)const{
Matrix R;
for(int i=0;i<N;++i){
for(int j=0;j<N;++j){
R.a[i][j]=a[i][j]+A.a[i][j];
if(R.a[i][j]>=mod) R.a[i][j]-=mod;
}
}
return R;
}
Matrix operator*(const Matrix& A)const{
Matrix R;
for(int i=0;i<N;++i){
for(int k=0;k<N;++k){
for(int j=0;j<N;++j){
R.a[i][j] = (R.a[i][j]+a[i][k]*A.a[k][j])%mod;
}
}
}
return R;
}
void print(){
for(int i=0;i<N;++i){
for(int j=0;j<N;++j){
cout<<a[i][j]<<" ";
}
cout<<endl;
}
}
};
Matrix pow(Matrix A,LL n){
Matrix R(1);
while(n){
if(n&1) R=R*A;
n>>=1; A=A*A;
}
return R;
}

矩阵类(版本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
LL mod;
class Matrix{
public:
const static int N = 102;
LL a[N][N];
int n;
Matrix(){}
Matrix(int _n,int x=0):n(_n){ // xIn
all(0);
for(int i=0;i<n;++i){
a[i][i]=x;
}
}
void all(int x){
for(int i=0;i<n;++i){
for(int j=0;j<n;++j){
a[i][j]=x;
}
}
}
Matrix operator+(const Matrix& A)const{
Matrix R(n);
for(int i=0;i<n;++i){
for(int j=0;j<n;++j){
R.a[i][j]=a[i][j]+A.a[i][j];
if(R.a[i][j]>=mod) R.a[i][j]-=mod;
}
}
return R;
}
Matrix operator*(const Matrix& A)const{
Matrix R(n);
for(int i=0;i<n;++i){
for(int k=0;k<n;++k){
for(int j=0;j<n;++j){
R.a[i][j] = (R.a[i][j]+a[i][k]*A.a[k][j])%mod;
}
}
}
return R;
}
void print(){
for(int i=0;i<n;++i){
for(int j=0;j<n;++j){
cout<<a[i][j]<<" ";
}
cout<<endl;
}
}
};
Matrix pow(Matrix A,int n){
Matrix R(A.n,1);
while(n){
if(n&1) R=R*A;
n>>=1; A=A*A;
}
return R;
}

注意到,矩阵乘法一定要写成上面的循环形式,这样利用高速缓存执行时间是原有的1/4

另外有序数组的累和要比无序的快很多,也是因为高速缓存(这个不太懂原理)

数据结构

并查集

1
2
3
4
5
6
7
8
9
10
int find(int x){
int ans = x;
while(ans!=p[ans]) ans=p[ans];
while(x!=ans){
int t = p[x];
p[x] = ans;
x = t;
}
return ans;
}

树状数组

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
struct TreeArray{
LL s[N];
int Size;
TreeArray(){}
TreeArray(int _S):Size(_S){clr(s,0);}
int lowbit(int n){
return n&(-n);
}
void add(int id,int p){
while(id<=Size){
s[id]+=p;
id+=lowbit(id);
}
}
LL sum(int id){
LL r=0;
while(id){
r+=s[id];
id-=lowbit(id);
}
return r;
}
};

若设原始数组为 $a$,设数字 $i$ 的二进制表示为 $d_1 \cdots d_n0\dots 0,d_n = 1$ 本质上树状数组 $s[i]$ 保存的其实是前 $n-1$ 位和 $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
#define lrt rt<<1
#define rrt rt<<1|1
#define lson l,m,lrt
#define rson m+1,r,rrt
const int N=1e5+5;
LL sum[N*3],col[N*3];
void PushUp(int rt){
sum[rt]=sum[lrt]+sum[rrt];
}
void PushDown(int rt,int m){
if(col[rt]){
col[lrt]+=col[rt];
col[rrt]+=col[rt];
sum[lrt]+=(m-(m>>1))*col[rt];
sum[rrt]+=(m>>1)*col[rt];
col[rt]=0;
}
}
void build(int l,int r,int rt){
if(l==r){
scanf("%I64d",&sum[rt]);
return;
}
col[rt]=0;
int m=(l+r)>>1;
build(lson);
build(rson);
PushUp(rt);
}
void update(int L,int R,int p,int l,int r,int rt){
if(L<=l&&R>=r){
sum[rt]+=p*(r-l+1);
col[rt]+=p;
return;
}
PushDown(rt,r-l+1);
int m=(l+r)>>1;
if(L<=m) update(L,R,p,lson);
if(R>m) update(L,R,p,rson);
PushUp(rt);
}
LL query(int L,int R,int l,int r,int rt){
if(L<=l&&R>=r) return sum[rt];
PushDown(rt,r-l+1);
int m=(l+r)>>1;
LL ans=0;
if(L<=m) ans+=query(L,R,lson);
if(R>m) ans+=query(L,R,rson);
return ans;
}

RMQ 求区间最大值

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
class RMQ{
static const int N=100005;
public:
int n;
int a[N][20]; //a[i][j]=max(s[i]---s[i+2^j-1])
RMQ(int* s,int _n):n(_n){
for(int i=0;i!=n;++i) a[i][0]=s[i];
int len=(int)(log(double(n))/log(2.0))+2;
for(int i=1;i!=len;++i){
for(int j=0;j<=n-(1<<i);++j){
a[j][i]=max(a[j][i-1],a[j+(1<<(i-1))][i-1]);
}
}
}
int Max(int l,int r){ // 0 <= l <= r < n
int len=(int)(log(double(r-l+1))/log(2.0));
return max(a[l][len],a[r-(1<<len)+1][len]);
}
};

最长(严格)递增子序列

1
2
3
4
5
6
7
8
9
10
11
12
13
int LIS(int a[],int n){ // longest increasing subsquence
if(n<=0) return 0;
int *b = new int[n];
b[0]=a[0];
int k=0;
for(int i=1;i!=n;++i){
if(a[i]>b[k]) b[++k]=a[i];
else b[lower_bound(b,b+k,a[i])-b]=a[i];
}
return k+1;
}
// lower_bound(first,end,val) 表示在单增 [frist,end) 中首次大于等于 val 的位置
// upper_bound(first,end,val) 表示在单增 [frist,end) 中首次大于 val 的位置

我们经常用二分答案的思想,但是其实二分答案是仅仅知道其单调的情况下的策略,实际上,对于具体的问题,我们完全可以对 $m$ 的值进行不同的处理,而非单纯的 $m=(l+r)>>1$。

最大子序列和

1
2
3
4
5
6
7
8
9
10
int MSCS(int a[],int n){ // maximal sum of continue subsquence,mind overflow
if(n<=0) return 0;
int r=a[0],s=a[0];
for(int i=1;i!=n;++i){
s = max(s,0);
s += a[i];
r = max(r,s);
}
return r;
}

最多区间交

做法:左端点标记 1, 右端点标记 -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
const int MAX=1e5+5;
int r[MAX];
void pack(int cash,int num,int v,int w){
if(num==0||w==0||v==0) return;
if(num==1){ // 0-1背包
for(int i=cash;i>=v;--i)
r[i]=max(r[i],r[i-v]+w);
return;
}
if(num*v>=cash-v+1){ //完全背包
for(int i=v;i<=cash;++i)
r[i]=max(r[i],r[i-v]+w);
return;
}
int q[MAX],s[MAX],head,tail;
for(int j=0;j<v;++j){ //多重背包
q[0]=r[j];s[0]=head=tail=0;
for(int i=1,k=j+v;k<=cash;++i,k+=v){
q[i]=r[k]-i*w;
while(s[head]<i-num) ++head;
while(head<=tail&&q[tail]<q[i]) --tail;
s[++tail]=i;
q[tail]=q[i];
r[k]=q[head]+i*w;
}
}
}

字符串匹配的Sunday算法

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
#include<cstring>
const int ASSIZE=256;
int tp[ASSIZE];
void getTp(char *s){
int sLen=strlen(s);
for(int i=0;i!=ASSIZE;++i){
tp[i]=sLen+1;
}
for(const char *p=s;*p;++p){
tp[*p]=sLen-(p-s);
}
return;
}
int Sunday(char *ps,char *s){
getTp(s);
const char *t,*p,*tx=ps;
int pLen=strlen(ps);
int sLen=strlen(s);
while(tx+sLen<=ps+pLen){
for(t=tx,p=s;*p;++p,++t){
if(*p!=*t) break;
}
if(*p=='\0') return tx-ps;
tx+=tp[tx[sLen]];
}
return -1;
}

字符串匹配的KMP算法

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
void getNext(char *s,int *next){
int j=0,k=-1;
int sLen=strlen(s)-1;
next[0]=-1;
while(j<sLen){
if(k==-1||s[j]==s[k]){
next[++j]=++k;
}else{
k=next[k];
}
}
return;
}
int KMP(char *p,char *s){
int pLen=strlen(p);
int sLen=strlen(s);
int *next=new int[sLen];
getNext(s,next);
while(i<pLen&&j<sLen){
if(j==-1||p[i]==s[j]){
++i;++j;
}else{
j=next[j];
}
}
delete [] next;
if(j==sLen) return i-j;
return -1;
}

整数字典树 Trie

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
struct Trie{
Trie *nt[2];
int cnt;
Trie() { nt[0] = nt[1] = NULL; cnt = 0; }
};
const int Maxb = 31;
void insert(Trie* root, int x, int add = 0){
root->cnt += add;
for(int i=Maxb-1;i>=0;--i){
if(root->nt[(x>>i)&1]==NULL) root->nt[(x>>i)&1] = new Trie();
root = root->nt[(x>>i)&1];
root->cnt += add;
}
}
int getans(Trie* root,int x,int k){ // 与 x 异或后不小于 k 的个数
int ans = 0;
for(int i=Maxb-1;i>=0;--i){
Trie* t = root->nt[((x>>i)^1)&1];
root = root->nt[((k^x)>>i)&1];
if(!((k>>i)&1) && t != NULL) ans += t->cnt;
if(root == NULL) break;
}
if(root != NULL) ans += root->cnt;
return ans;
}
void clear(Trie *root){ // 多case时需要清理内存
if(root->nt[0] != NULL) clear(root->nt[0]);
if(root->nt[1] != NULL) clear(root->nt[1]);
delete root;
}

数组型字典树

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
const int N = 1e4+4;
const int Maxb = 31;
int trie[N*Maxb][2],sum[N*Maxb],tt; // tt need init
void add(int x,int p){
for(int i=Maxb,t=0; i>=0; --i){
int &tmp = trie[t][x>>i &1];
if(!tmp) tmp = ++tt;
t = tmp;
sum[t]+=p;
}
}
bool find(int x,int k){ // if exist y such that y^x >= k
for(int i=Maxb,t=0,tmp; i>=0; --i){
tmp = trie[t][!(x>>i &1)];
if(k>>i &1){
if(!sum[tmp]) return 0;
}else{
if(sum[tmp]) return 1;
tmp = trie[t][x>>i &1];
if(!sum[tmp]) return 0;
}
t = tmp;
}
return 1;
}

字母字典树 Trie

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
struct Trie{
const static int MAX=26;
bool isStr;
Trie *next[MAX];
Trie(){
isStr=false;
for(int i=0;i!=MAX;++i) next[i]=NULL;
}
};
void insert(Trie *root,const char *s){
if(root==NULL||*s=='\0') return;
while(*s!='\0'){
if(root->next[*s-'a']==NULL){
root->next[*s-'a']=new Trie;
}
root=root->next[*s-'a'];
++s;
}
root->isStr=true;
}
bool search(Trie *root,const char *s){
while(root!=NULL&&*s!='\0'){
root=root->next[*s-'a'];
++s;
}
return (root!=NULL&&root->isStr);
}
void clear(Trie *root){
for(int i=0;i!=Trie::MAX;++i){
if(root->next[i]!=NULL) clear(root->next[i]);
}
delete root;
}

单调队列和单调栈

单调队列,单调栈顾名思义就是内部元素是单调的,区别就是队列是先进先出,栈是后进先出。

  1. 单调队列的典型应用就是 $O(n)$复杂度求给定宽度的区间在各处的最大值。
  2. 单调栈的典型应用是 $O(n)$ 复杂度求出以每一点为最小值的最大区间。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
#include <bits/stdc++.h>
const int N = 1e5+5;
int L[N],R[N],stk[N]; //L[i]: min x<i s.t. h[x]>h[i]
void monoStack(int h[], int n){
int top = 0;
for(int j = 0; j < n; ++j) {
while(top && h[stk[top-1]] >= h[j]) --top;
L[j] = (top == 0?-1:stk[top-1]);
stk[top++] = j;
}
top = 0;
for(int j = n-1; j >= 0; --j) {
while(top && h[stk[top-1]] >= h[j]) --top;
R[j] = top == 0?n:stk[top-1];
stk[top++] = j;
}
}

优先队列

可以使用C++ STL 的 priority_queue,查找可用 lower_bound 和 upper_bound。C++ STL 中优先队列是用堆来实现的。用途十分广泛,例如加速最小生成树,拓扑排序,等等。

堆的实现一般是用数组。我们可以用 1 作为树的根,对每一个节点 $x$ ,它的两个节点分别就是 $2x$ 和 $2x+1$ 平时都用 $x<<1,x<<1|1$ 表示。 堆只支持三个操作,

  1. 插入一个节点(我们实现时是插入最尾部,这样保证了是一个完全二叉树) $O(\log n)$
  2. 删除最大键值节点(删除根元素的值)$O(\log n)$
  3. 输出最大键值节点(查看根元素的值)$O(1)$

堆维持了一个性质使其能如此高效:父节点的值不小于儿子节点的值。那么显然根节点的键值最大。
那么我们在尾部插入一个节点,如果它大于其父节点,那么跟父节点互换。因为树的高度为 $O(\log n)$, 因此插入操作复杂度也就 $O(\log n)$, 删除最大节点就是先让其和最尾部的节点互换,然后删除最尾部节点,然后从头到尾的与键值大的子节点互换,直到维护了堆的性质,最后查看最大键值节点就是看根节点,显然是 $O(1)$ 的。
它牺牲了许多操作,带来了高效率。

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
const int M = 200005;
const int N = 100005;
int head[N],sc,d[N],ans[N];
bool v[N];
struct Node{
int ed;
int next;
}e[M];
void addedge(int x,int y){
e[sc].ed=y;
e[sc].next=head[x];
head[x]=sc++;
++d[x];
}
void init(int n){
sc = 1;
for(int i=1;i<=n;++i){
head[i]=-1;
d[i]=0;
v[i] = false;
}
}
void topoSort(int n){ // after read data
priority_queue<int> q;
for(int i=1;i<=n;++i){
if(d[i] == 0) q.push(i);
}
while(!q.empty()){
int u = q.top();
q.pop();
if(v[u]) continue;
v[u] = true;
for(int i=head[u];i!=-1;i=e[i].next){
if(--d[e[i].ed] == 0) q.push(e[i].ed);
}
}
}

红黑树 red-black tree

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
class RBT{
typedef int elemType;
struct Node{
const static bool RED = 0;
const static bool BLACK = 1;
Node *ch[2], *fa;// x->fa->ch[x->rs] = x
int sz;
elemType key;
bool color, rs; // is rightson
};
Node *root; // root has no father
void faSon(Node* x, Node* y, bool rs){
y->fa = x;
y->rs = rs;
x->ch[rs] = y;
}
void newNode(Node* x, elemType val, bool rs){ // x->ch[rs]=null
Node *p = new Node;
p->ch[0] = p->ch[1] = null;
p->sz = 1; p->key = val; p->color = Node::RED;
faSon(x, p, rs);
++x->sz;
}
void rotate(Node* x, bool rs){ // x must not null
Node *y = x->ch[!rs];
if(y == null) return;
if(x == root) root = y;
else faSon(x->fa, y, x->rs);
faSon(x, y->ch[rs], !rs);
faSon(y, x, rs);
y->sz = x->sz;
x->sz = x->ch[0]->sz + x->ch[1]->sz + 1;
}
void insMaintain(Node* x){ // x->color is RED
if(x == root || x->fa->color == Node::BLACK) return;
if(x->fa->fa->ch[!x->fa->rs]->color == Node::BLACK){
if(x->rs ^ x->fa->rs) rotate(x->fa, x->fa->rs);
else x = x->fa;
x->color = Node::BLACK;
x->fa->color = Node::RED;
rotate(x->fa, !x->rs);
}else{
x = x->fa->fa;
x->color = Node::RED;
x->ch[0]->color = x->ch[1]->color = Node::BLACK;
insMaintain(x);
}
}
void delCase1(Node* x, Node* y){
y->color = Node::BLACK;
y->fa->color = Node::RED;
y = y->ch[!y->rs];
rotate(x->fa, x->rs);
delCase2(x, y);
}
void delCase2(Node *x, Node *y){
if(y->ch[y->rs]->color == Node::BLACK){
if(y->ch[!y->rs]->color == Node::BLACK){
y->color=Node::RED;
delMaintain(y->fa);
}else{
y->color = Node::RED;
y->ch[!y->rs]->color = Node::BLACK;
rotate(y, y->rs);
delCase3(y->fa);
}
}else delCase3(y);
}
void delCase3(Node *y){
y->color = y->fa->color;
y->ch[y->rs]->color = y->fa->color = Node::BLACK;
rotate(y->fa,!y->rs);
}
void delMaintain(Node* x){
if(x == root || x == null) return;
if(x->color == Node::RED){
x->color = Node::BLACK; return;
}
Node *y = x->fa->ch[!x->rs];
if(y->color == Node::RED) delCase1(x, y);
else delCase2(x, y);
}
Node* pred(Node* x, elemType val){ // max elem <= val
if(x == null) return null;
if(x->key > val) return pred(x->ch[0], val);
else if(x->ch[1] == null) return x;
return pred(x->ch[1], val);
}
Node* succ(Node* x, elemType val){ // min elem >= val
if(x == null) return null;
if(x->key < val) return succ(x->ch[1], val);
else if(x->ch[0] == null) return x;
return succ(x->ch[0], val);
}
int rank(Node* x, elemType val){ // count elem <= val
if(x == null) return 0;
if(x->key > val) return rank(x->ch[0], val);
return x->ch[0]->sz + 1 + rank(x->ch[1], val);
}
Node* select(Node* x, int k){ // k-th smallest elem
if(x == null || x->sz < k) return null;
int sz = x->ch[0]->sz + 1;
if(sz == k) return x;
if(sz < k) return select(x->ch[1], k-sz);
return select(x->ch[0], k);
}
void clear(Node* x){
if(x->ch[0] != null) clear(x->ch[0]);
if(x->ch[1] != null) clear(x->ch[1]);
if(x != null) delete x;
}
void print(Node* x){
printf("key = %d, sz = %d ", x->key, x->sz);
puts(x->color == Node::RED?"RED":"BLACK");
if(x->ch[0] != null) print(x->ch[0]);
if(x->ch[1] != null) print(x->ch[1]);
}
public:
const static int INF = 0x3f3f3f3f;
Node *null;
RBT(){
null = new Node; // no key, rs, father
null->ch[0] = null->ch[1] = null;
null->sz = 0; null->color = Node::BLACK;
root = null; null->key = INF; // for convenient
}
Node* search(elemType val){
Node *x = root;
while(x != null){
if(val == x->key) return x;
x = x->ch[val >= x->key];
}
return null;
}
void insert(elemType val){
if(root == null){
root = new Node; // no father, rs
root->ch[0] = root->ch[1] = null;
root->sz = 1; root->color = Node::BLACK;
root->key = val;
return;
}
Node *x = root;
while(x->ch[val >= x->key] != null){
++x->sz;
x = x->ch[val >= x->key];
}
newNode(x, val, val >= x->key);
insMaintain(x->ch[val >= x->key]);
root->color = Node::BLACK;
}
void del(elemType val){
Node *x = search(val), *y;
if(x == null) return;
while(x->ch[0] != null || x->ch[1] != null){
bool rs = 0;
if(x->ch[rs] == null) rs = !rs;
y = x->ch[rs];
while(y->ch[!rs] != null) y = y->ch[!rs];
swap(x->key, y->key);
x = y;
if(x->color == Node::RED) break;
}
delMaintain(x);
root->color = Node::BLACK;
y = x;
while(y != root){
y = y->fa;
--y->sz;
}
if(x == root) root = null;
else if(x->ch[0] != null) faSon(x->fa, x->ch[0], x->rs);
else if(x->ch[1] != null) faSon(x->fa, x->ch[1], x->rs);
else x->fa->ch[x->rs] = null;
delete x;
}
elemType pred(elemType val){
return pred(root, val)->key;
}
elemType succ(elemType val){
return succ(root, val)->key;
}
int rank(elemType val){
return rank(root, val);
}
elemType select(int k){
return select(root, k)->key;
}
void clear(){
clear(root);
root = null;
}
int size(){
return root->sz;
}
void print(){
if(root != null) print(root);
}
int getans(int a,int k){ // for particular use
return select(root, rank(root, a) + k)->key;
}
};

图论

知乎上看到YYYYLLL 关于Floyd 算法的解释挺好的,再次记录(稍加修改):

1
2
3
4
DP[k][i][j] 表示只经过 1~k 号节点优化,i 点到 j 点的最短路径长度。
则 DP[k][i][j] = min( DP[k-1][i][j], DP[k-1][i][k]+DP[k-1][k][j] )
= min( DP[k-1][i][j], DP[k][i][k]+DP[k][k][j] )
DP[0][][] 是初始图的邻接矩阵,DP[n][][] 就是最终求得的最短路长度矩阵了

本来一开始是没法做空间优化的,但是第二个等式,就保证了可以做空间优化

C++代码如下:

1
2
3
4
5
6
7
8
const int N = 102;
LL dp[N][N];
void Floyd(int n){
for(int k = 0; k!=n;++k)
for(int i<0; i!=n;++i)
for(int j=0;j!=n;++j)
dp[i][j] = min(dp[i][j],dp[i][k] + dp[k][j]);
}

dfs 序

1
2
3
4
5
6
7
8
ncnt = 0; // init every case
void dfs(int x, int fa){ // dfs order
L[x] = ++ncnt;
for(int i=head[x]; i!=-1; i=e[i].next){
if(e[i].ed != fa) dfs(e[i].ed, x);
}
R[x] = ncnt;
}

最大流

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
struct Node{
int t,w,next;
}a[N*N];
int n,ss,head[N],p[N],flow[N],c[N],h[N],numh[N];
void addedge(int x,int y,int w){
a[ss].t=y;
a[ss].w=w;
a[ss].next=head[x];
head[x]=ss++;
}
int max_flow(int s,int t){
int flow,ans=0,neck,k;
memset(h,0,sizeof(h));
memset(numh,0,sizeof(numh));
memset(p,-1,sizeof(p));
for(int i=1;i<=n;++i) c[i]=head[i];
numh[0]=n;
int u=s;
while(h[s]<n){
if(u==t){
flow=1e9;
for(int i=s;i!=t;i=a[c[i]].t){
if(flow>a[c[i]].w){
neck=i;flow=a[c[i]].w;
}
}
for(int i=s;i!=t;i=a[c[i]].t){
a[c[i]].w-=flow;
a[c[i]^1].w+=flow;
}
ans+=flow;
u=neck;
}
for(k=c[u];k!=-1;k=a[k].next){
if(a[k].w&&h[u]==h[a[k].t]+1) break;
}
if(k!=-1){
c[u]=k;
p[a[k].t]=u;
u=a[k].t;
}else{
if(0==--numh[h[u]]) break;
c[u]=head[u];
k=n;
for(int i=head[u];i!=-1;i=a[i].next){
if(a[i].w) k=min(k,h[a[i].t]);
}
h[u]=k+1;
++numh[h[u]];
if(u!=s) u=p[u];
}
}
return ans;
}

Stoer-Wagner 最小割

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
const int N=102;
int p[N],dis[N],map[N][N];
bool vis[N];
int mincut(int n){
int ret=0x7fffffff;
for(int i=0;i!=n;++i) p[i]=i;
while(n>1){
int t=1,s=0;
for(int i=1;i!=n;++i){
dis[p[i]]=map[p[0]][p[i]];
if(dis[p[i]]>dis[p[t]]) t=i;
}
memset(vis,0,sizeof(vis));
vis[p[0]]=true;
for(int i=1;i<n;++i){
if(i==n-1){
if(ret>dis[p[t]]) ret=dis[p[t]];
for(int j=0;j!=n;++j){
map[p[j]][p[s]]+=map[p[j]][p[t]];
map[p[s]][p[j]]=map[p[j]][p[s]];
}
p[t]=p[--n];
}
vis[p[t]]=true;
s=t;t=-1;
for(int j=1;j!=n;++j){
if(!vis[p[j]]){
dis[p[j]]+=map[p[j]][p[s]];
if(t==-1||dis[p[j]]>dis[p[t]]) t=j;
}
}
}
}
return ret;
}

最短路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
const int INF=0x3f3f3f3f;
const int M=2000005;
const int N=1004;
int head[N],rehead[N],dist[N],sc,v[N];
struct A{
int id;
int g;
A(){}
A(int x,int y):id(x),g(y){}
bool operator<(const A& a)const {
return g+dist[id]>a.g+dist[a.id];
}
};
struct Node{
int ed;
int w;
int next;
}e[M];
void addedge(int x,int y,int z){
e[sc].ed=y;
e[sc].w=z;
e[sc].next=head[x];
head[x]=sc++;
e[sc].ed=x;
e[sc].w=z;
e[sc].next=rehead[y];
rehead[y]=sc++;
}
void spfa(int s){
memset(v,0,sizeof(v));
memset(dist,0x3f,sizeof(dist));
dist[s]=0;
v[s]=true;
queue<int> q;
q.push(s);
while(!q.empty()){
int u=q.front();
for(int i=rehead[u];i!=-1;i=e[i].next){
int ed=e[i].ed;
if(dist[ed]>dist[u]+e[i].w){
dist[ed]=dist[u]+e[i].w;
if(!v[ed]){
v[ed]=true;
q.push(ed);
}
}
}
v[u]=false;
q.pop();
}
}

几何

二维凸包

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
const int N=501;
int vc[N];
struct Node{
int x,y,id;
bool operator!=(const Node& A)const{
return x!=A.x||y!=A.y;
}
bool operator<(const Node& A)const{
if(y==A.y) return x<A.x;
return y<A.y;
}
}p[N],q[N];
bool crossLeft(const Node& op,const Node& sp,const Node& ep){
return (sp.x-op.x)*(ep.y-op.y)>(sp.y-op.y)*(ep.x-op.x);
}
int graham(int n){
sort(p,p+n);
if(n==0) return 0;q[0]=p[0];
if(n==1) return 1;q[1]=p[1];
if(n==2) return 2;q[2]=p[2];
int top=1;
for(int i=2;i!=n;++i){
while(top&&crossLeft(q[top],p[i],q[top-1])) --top;
q[++top]=p[i];
}
int len=top;
q[++top]=p[n-2];
for(int i=n-3;i!=-1;--i){
while(top!=len&&crossLeft(q[top],p[i],q[top-1])) --top;
q[++top]=p[i];
}
return top;
}

几类根号算法

$ s(n) = \sum_{i=1}^{n} \lfloor \frac{n}{i} \rfloor $

1
2
3
4
5
6
7
8
LL getsum(LL n){ 
LL sum = 0;
for(LL i=1,j;i<=n;i=j+1){
j = n/(n/i);
sum += (j-i+1)*(n/i);
}
return sum;
}

$\sum_{i=1}^n \lfloor \frac{n}{i} \rfloor \lfloor \frac{m}{i} \rfloor f(i)$

1
2
3
4
5
6
7
8
LL getsum(int n,int m){ // g[i]=f[i]+g[i-1]
LL sum = 0;
for(int i=1,j;i<=min(n,m);i=j+1){
j=min(n/(n/i),m/(m/i)
sum += LL(n/i)*(m/i)*(g[j]-g[i-1]);
}
return sum;
}

$h(n) = \frac{n(n-1)(n-2)}{3} - \sum_{i=2}^n h(\lfloor \frac{n}{i} \rfloor)$

1
2
3
4
5
6
7
8
9
10
11
12
13
map<int,int> mp;
int getans(int n){
map<int,int>::iterator it = mp.find(n);
if (it != mp.end()) return it->second;
int r = LL(n)*(n-1)%M*(n-2)%M*inv3%M;
for(int i=2,j;i<=n;i=j+1){
j = n/(n/i);
r -= LL(j-i+1)*getans(n/i)%M;
if(r<0) r+=M;
}
mp.insert(pair<int,int>(n,r));
return r;
}

To Be Continue

现代C++ : C++ 17

自2020/5/18 从传统C++ 转 C++17, 有问题多查 C++ Reference

  • auto :让编译器自动类型判断,大大简化了模板的写法,为Lambda函数提供了基础

  • pair, tuple:用 get<i>获取单个的值,也可以用Python形式的解析(std::tie 应该被抛弃)

  • 省略参数运算(计算时可自定义左右结合)

  • if switch 现在也支持局部变量了!

  • constexprif constexpr 让编译器自己去优化

  • Lambda 函数 特别好用!配上STL 真香。

    1
    2
    3
    4
    5
    6
    7
    8
    9
    [capture list] (parameters)
    mutable (optional)
    constexpr (optional)
    exception attr (optional)
    -> return type (optional)
    {
    body
    }
    // capture list : [], [&](引用), [=](copy), [x = 1](加入初始自定义量)
  • STL 更好用了!

没有必要关心你自己不用的特征!用的就一定要搞清楚。

日常表白zly

欧尼酱,人家想喝可乐!