C++ Template dna049

发表于 |

突然想到在此留下模板貌似是个不错的选择也 0.0 长期更新

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

1
Composed by dna049 at 2016-4-10 in FDU

通用代码块

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
//#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>
using namespace std;
typedef long long LL;
#if ( ( _WIN32 || __WIN32__ ) && __cplusplus < 201103L)
#define lld "%I64d"
#else
#define lld "%lld"
#endif
template<typename T>
void upmax(T &a,T b){ if(a<b) a=b;}
//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(){
//#define dna049 1
#ifdef dna049
    freopen("/Users/dna049/Desktop/AC/in","r",stdin);
    LL time = clock();
#endif
//    freopen("/Users/dna049/Desktop/AC/out","w",stdout);
#ifdef dna049
    printf("time: %f\n",1.0*(clock()-time)/CLOCKS_PER_SEC);
#endif
    return 0;
}
/*
                   _ooOoo_
                  o8888888o
                  88" . "88
                  (| -_- |)
                  O\  =  /O
               ____/`---'\____
             .'  \\|     |//  `.
            /  \\|||  :  |||//  \
           /  _||||| -:- |||||-  \
           |   | \\\  -  /// |   |
           | \_|  ''\---/''  |   |
           \  .-\__  `-`  ___/-. /
         ___`. .'  /--.--\  `. . __
      ."" '<  `.___\_<|>_/___.'  >'"".
     | | :  `- \`.;`\ _ /`;.`/ - ` : | |
     \  \ `-.   \_ __\ /__ _/   .-` /  /
======`-.____`-.___\_____/___.-`____.-'======
                   `=---='
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
            佛祖保佑       永无BUG
*/

数论篇

Greatest Common divisor

1
2
3
4
5
6
7
8
9
LL gcd(LL a,LL b){
	return b?gcd(b,a%b):a;
}
LL exgcd(LL a,LL b,LL& x,LL& y){
    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
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
LL pow(LL x,int n){ 
    LL r=1;for(;n;n>>=1,x*=x)  if(n&1) r*=x;return r;
}
LL pow_mod(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 mul_mod(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 mul_mod(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 pow_mod(LL x,LL n,LL p){
    LL r=1;
    while(n){
        if(n&1) r=mul_mod(r,x,p);
        n>>=1;  x=mul_mod(x,x,p);
    }
    return r;
}

对于 pow_sum 可以化为矩阵幂解决，也可以用等比数列公式，取大的模解决。

整数开根号

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] = pow_mod(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]=pow_mod(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;
}

离散对数

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 = pow_mod(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(pow_mod(a,q,p)!=1)   return -1;
    if(q&1) return pow_mod(a,(q+1)>>1,p);
    LL b,cnt=1;
    while(pow_mod(b=rand()%p,q,p)==1);//find a non quadratic residue
    while(!(q&1))   ++cnt,q>>=1;
    b=pow_mod(b,q,p);
    LL x=pow_mod(a,(q+1)>>1,p);
    for(LL s=1,t=pow_mod(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=pow_mod(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(pow_mod(a,q,n)!=1)   return -1;
    if(q&1) return pow_mod(a,(q+1)>>1,n);
    LL b,cnt=1;
    while(pow_mod(b=rand()%n,q,n)==1);//find a non quadratic residue
    while(!(q&1))   ++cnt,q>>=1;
    b=pow_mod(b,q,n);
    LL x=pow_mod(a,(q+1)>>1,n);
    for(LL s=1,t=pow_mod(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 = pow_mod(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;
}

详细分析见我的博文。

N以内素数线性筛 $O(n)$

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
const int N = 1e7+2;
bool ip[N];
int p[N];
int getprime(){
    int cnt=0;ip[2]=true;p[cnt++]=2;
    for(int i=1;i<N;i+=2)   ip[i]=true;
    for(int i = 3; i < N; i+=2) {
        if(ip[i])   p[cnt++] = i;
        for(int j = 1;j < cnt && i * p[j] < N; ++j) {
            ip[i * p[j]] = false;
            if(i % p[j] == 0)	break;
        }
    }
    return cnt;
}
bool isp(int n){
	if(n<N)	return ip[n];
	for(int i=0;p[i]*p[i]<=n;++i){
		if(n%p[i]==0)	return false;	
	}
	return true;
}

正整数拆分

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
bool Witness(LL a,LL n,LL m,int t){
    LL x=pow_mod(a,m,n);
    if(x==1||x==n-1)    return false;
    while(t--){
        x=mul_mod(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
const int N = 1e6+2;
int sp[N],p[N];
int spf(){ // samllest prime factor
    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];
        }
    }
    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 Pollard_rho(LL n){
    LL x=rand()%(n-1)+1;
    LL y=x,i=1,k=2,c=rand()%(n-1)+1;
    while(true){
        x=(mul_mod(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=Pollard_rho(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
const int N = 1e7+2;
bool ip[N];
int mu[N],p[N];
void init_mu(){
    mu[1]=1;ip[2]=true;p[0]=2;
    for(int i=1;i<N;i+=2)   ip[i]=true;
    for(int i = 3,cnt = 1;i<N;i+=2){
        if(ip[i]){
            p[cnt++] = i;
            mu[i] = -1;
        }
        for(int j = 1,t;j<cnt&&(t= i * p[j])<N;++j){
            ip[t] = false;
            if(i % p[j] == 0)   break;
            mu[t] = -mu[i];
        }
    }
    for(int i=2;i<N;i+=4)  mu[i]=-mu[i>>1];
}
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 init_sumu(){
    if(mu[1]!=1)    init_mu();
    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
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
const int N = 1e7+2;
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(){
    getprime();
    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 getphi(LL x,int s){
    if(s == 0)  return x;
    if(s <= M)  return phi[x%sz[s]][s]+(x/sz[s])*phi[sz[s]][s];
    if(x <= p[s]*p[s])   return pi[x]-s+1;
    if(x <= p[s]*p[s]*p[s] && x< N){
        int s2x = pi[sqrt2(x)];
        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 getphi(x,s-1)-getphi(x/p[s],s-1);
}
LL getpi(LL x){
    if(x < N)   return pi[x];
    LL ans = getphi(x,pi[sqrt3(x)])+pi[sqrt3(x)]-1;
    for(int i=pi[sqrt3(x)]+1,ed=pi[sqrt2(x)];i<=ed;++i){
        ans -= getpi(x/p[i])-i+1;
    }
    return ans;
}
LL lehmer_pi(LL x){ // x < N*N
    if(x < N)   return pi[x];
    int a = (int)lehmer_pi(sqrt2(sqrt2(x)));
    int b = (int)lehmer_pi(sqrt2(x));
    int c = (int)lehmer_pi(sqrt3(x));
    LL sum = getphi(x, a) + LL(b + a - 2) * (b - a + 1) / 2;
    for (int i = a + 1; i <= b; i++) {
        LL w = x / p[i];
        sum -= lehmer_pi(w);
        if (i > c) continue;
        LL lim = lehmer_pi(sqrt2(w));
        for (int j = i; j <= lim; j++) {
            sum -= lehmer_pi(w / p[j]) - (j - 1);
        }
    }
    return sum;
}

另一个版本：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
const int N = 1e7+2;
LL L[N],R[N];
LL getans(LL n){ // n*n < 1e14
    LL rn = sqrt2(n);
    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];
        }
    }
    LL ans = L[sqrt3(n)];
    for(LL p=2;p<=rn;++p){
        if(L[p] == L[p-1])  continue;
        ans += R[p]-L[p];
    }
    return ans;
}

求奇素数的一个原根

首先，模 $m$ 有原根的充要条件是：$m=2,4,p^a,2p^a$，其中$p$为奇素数。
对于求模$p$的原根方法：对 $p-1$ 素因子分解：$p-1 = p_1^{a_1} \cdots p_s^{a_s}$ 若恒有
$$ g^{\frac{p-1}{p_i}} \neq 1(mod \; p) $$
则 $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(pow_mod(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
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
LL chinaRemain(int n,LL a[],LL m[]){//x=a[i] mod m[i]
    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 = pow_mod(3,(FM-1)>>loglen,FM);
    if(isInverse){
        g=inv(g,FM);
        LL invlen = pow_mod(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 = pow_mod(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
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
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;
    }
};

图论

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;
}

几类根号算法

1.$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;
}

2.$\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;
}

3.$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