基础模板

常用板子

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
template<class T>inline void MAX(T &x,T y){if(y>x)x=y;}
template<class T>inline void MIN(T &x,T y){if(y<x)x=y;}
template<class T>inline void rd(T &x){x=0;char o,f=1;while(o=getchar(),o<48)if(o==45)f=-f;do x=(x<<3)+(x<<1)+(o^48);while(o=getchar(),o>47);x*=f;
}
template<class T>inline void print(T x,bool op=1){static int top,stk[105];if(x<0)x=-x,putchar('-');if(x==0)putchar('0');while(x)stk[++top]=x%10,x/=10;while(top)putchar(stk[top--]+'0');putchar(op?'\n':' ');
}signed main(){#ifndef ONLINE_JUDGEfreopen("jiedai.in","r",stdin);
#endifreturn (0-0);
}

数学题常用板子

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
template<class T>inline void MAX(T &x,T y){if(y>x)x=y;}
template<class T>inline void MIN(T &x,T y){if(y<x)x=y;}
template<class T>inline void rd(T &x){x=0;char o,f=1;while(o=getchar(),o<48)if(o==45)f=-f;do x=(x<<3)+(x<<1)+(o^48);while(o=getchar(),o>47);x*=f;
}
template<class T>inline void print(T x,bool op=1){static int top,stk[105];if(x<0)x=-x,putchar('-');if(x==0)putchar('0');while(x)stk[++top]=x%10,x/=10;while(top)putchar(stk[top--]+'0');putchar(op?'\n':' ');
}
const int SIZE=1e6+5;
const int P=1e9+7;
//const int P=998244353;
int fac[SIZE],inv[SIZE];
int fast(int a,int b=P-2){int res=1;while(b){if(b&1)res=1ll*res*a%P;a=1ll*a*a%P;b>>=1;}return res;
}
int comb(int n,int m){if(n<0||m<0||n<m)return 0;return 1ll*fac[n]*inv[m]%P*inv[n-m]%P;
}
ll gcd(ll n,ll m){if(m==0)return n;return gcd(m,n%m);
}signed main(){#ifndef ONLINE_JUDGEfreopen("jiedai.in","r",stdin);
#endiffac[0]=1;for(int i=1;i<SIZE;i++)fac[i]=1ll*fac[i-1]*i%P;inv[SIZE-1]=fast(fac[SIZE-1]);for(int i=SIZE-2;i>=0;i--)inv[i]=1ll*inv[i+1]*(i+1)%P;return (0-0);
}

输出挂

可以输出__int128

template<class T>inline void print(T x,bool op=1){static int top,stk[105];if(x<0)x=-x,putchar('-');if(x==0)putchar('0');while(x)stk[++top]=x%10,x/=10;while(top)putchar(stk[top--]+'0');putchar(op?'\n':' ');
}

fread快读

fread函数从文件中读入

fread在读取换行的时候可能会出问题（慎用）

int Tot;
char buf[50000005],*ch=buf-1;
template<class T>inline void rd(T &x){int f=1;for(++ch;*ch<48&&ch<buf+Tot;++ch)if(*ch==45)f=-f;if(ch==buf+Tot)return;for(x=0;*ch>47;++ch)x=(x<<3)+(x<<1)+(*ch^48);x*=f;
}Tot=fread(buf,1,5e7,stdin);

高精度

struct INT{static const int base=1e4;vector<int>v;INT(ll x=0){while(x)v.push_back(x%base),x/=base;if(v.size()==0)v.push_back(0);}inline int size()const{return v.size();}inline void resize(int sz=1){v.clear();v.resize(sz,0);}inline void topzero(){while(v.size()>1&&v.back()==0)v.pop_back();}INT operator +(const INT &A)const{INT T;T.resize(max(size(),A.size())+1);for(int i=0;i<T.size()-1;i++){if(i<size())T.v[i]+=v[i];if(i<A.size())T.v[i]+=A.v[i];if(T.v[i]>=base)T.v[i+1]++,T.v[i]-=base;}T.topzero();return T;}INT operator *(const INT &A)const{INT T;T.resize(size()+A.size()+5);for(int i=0;i<size();i++)for(int j=0;j<A.size();j++){T.v[i+j]+=v[i]*A.v[j];if(T.v[i+j]>=2e9)T.v[i+j+1]+=T.v[i+j]/base,T.v[i+j]%=base;}for(int i=0;i<T.size();i++){if(T.v[i]>=base)T.v[i+1]+=T.v[i]/base,T.v[i]%=base;}T.topzero();return T;}INT operator /(const int x)const{INT res;res.resize(size());ll tmp=0;for(int i=size()-1;i>=0;i--){tmp=tmp*base+v[i];res.v[i]=tmp/x,tmp%=x;}res.topzero();return res;}int operator %(const int mod)const{int res=0;for(int i=size()-1;i>=0;i--)res=(1ll*res*base+v[i])%mod;return res;}friend INT fast(INT a,ll b){INT res(1);while(b){if(b&1)res=res*a;a=a*a;b>>=1;}res.topzero();return res;}void read(){string str;cin>>str;resize();for(int i=0;i<str.size();i++)*this=*this*INT(10)+INT(str[i]-'0');}void print(){printf("%d",v[size()-1]);for(int i=size()-2;i>=0;i--)printf("%04d",v[i]);}
};

分数类

适合打表时使用，运算范围在__int128以内

将define去掉后，运算范围在long long以内

struct frac{#define ll __int128ll a,b;frac(ll _a=0,ll _b=1):a(_a),b(_b){}ll gcd(ll a,ll b){if(b==0)return a;return gcd(b,a%b);}frac work(ll a,ll b){if(b<0)a=-a,b=-b;ll tmp=gcd(a,b);if(tmp<0)tmp=-tmp;return (frac){a/tmp,b/tmp};}frac operator +(frac A){return work(a*A.b+b*A.a,b*A.b);}frac operator -(frac A){return work(a*A.b-b*A.a,b*A.b);}frac operator *(frac A){return work(a*A.a,b*A.b);}frac operator *(ll x){return work(a*x,b);}frac operator /(frac A){return work(a*A.b,b*A.a);}frac operator /(ll x){return work(a,b*x);}void pt(ll x){static int stk[105],top;if(x==0)printf("0");if(x<0)x=-x,printf("-");while(x)stk[++top]=x%10,x/=10;while(top)printf("%d",stk[top--]);}void pt(){pt(a),printf("/"),pt(b),printf("  ");}#undef ll
};

打表压缩

代码一：数值小

/*利用44进制对打表的数组进行压缩，大于等于44的数表示分隔对于int范围以内的数压缩效果较好减号 '-'(45) 表示负号去掉 '"'(34) '/'(47) '?'(63) '\'(92)少了 ' '(32) '~'(126)
*/
void encode(ll *A,char *str,int n){for(int i=1,len=0;i<=n;i++){ll x=A[i];if(x<0)str[len++]='-',x=-x;if(x==0)str[len++]=81;while(x){int t=x%44;x/=44;if(x==0)t+=44;str[len++]=t+33+(t>=1)+(t>=11)+(t>=12)+(t>=27)+(t>=55);}}
}
int decode(char *str,ll *A){int n=0,len=strlen(str);for(int l=0,r=0,f=1;l<len;l=r+1,r=l,f=1){if(str[l]=='-')f=-1,l++,r++;while(r+1<len&&str[r]<81)r++;ll x=0;for(int i=r;i>=l;i--){int t=str[i]-33;t-=(t>=1),t-=(t>=11),t-=(t>=12),t-=(t>=27),t-=(t>=55);x=x*44+(t%44);}A[++n]=f*x;}return n;
}

代码二：数值大

/*利用89进制对打表的数组进行压缩，最大限度利用可显示字符对于long long范围以内的大整数压缩效果较好空格 ' '(32) 表示数字分隔符减号 '-'(45) 表示负号去掉 '"'(34) '/'(47) '?'(63) '\'(92)
*/
void encode(ll *A,char *str,int n){for(int i=1,len=0;i<=n;i++){ll x=A[i];if(x<0)str[len++]='-',x=-x;if(x==0)str[len++]=33;while(x){int t=x%89;x/=89;str[len++]=t+33+(t>=1)+(t>=11)+(t>=12)+(t>=27)+(t>=55);}str[len++]=' ';}
}
int decode(char *str,ll *A){int n=0,len=strlen(str);for(int l=0,r=0,f=1;l<len;l=r+2,r=l,f=1){if(str[l]=='-')f=-1,l++,r++;while(r+1<len&&str[r+1]!=' ')r++;ll x=0;for(int i=r;i>=l;i--){int t=str[i]-33;t-=(t>=1),t-=(t>=11),t-=(t>=12),t-=(t>=27),t-=(t>=55);x=x*89+t;}A[++n]=f*x;}return n;
}

基数排序

给 int 范围内的非负整数排序，建议 n>=1e5

void Sort(int n,int *A){static const int M=1e6+5,S=16,P=(1<<S)-1;static int cnt[1<<S],B[M];memset(cnt,0,sizeof(cnt));for(int i=1;i<=n;i++)cnt[A[i]&P]++;for(int i=1;i<=P;i++)cnt[i]+=cnt[i-1];for(int i=n;i>=1;i--)B[cnt[A[i]&P]--]=A[i];memset(cnt,0,sizeof(cnt));for(int i=1;i<=n;i++)cnt[B[i]>>S]++;for(int i=1;i<=P;i++)cnt[i]+=cnt[i-1];for(int i=n;i>=1;i--)A[cnt[B[i]>>S]--]=B[i];
}

杂项

优先队列小顶堆

priority_queue<int,vector<int>,greater<int> >Q;

二进制数中1的个数

inline int popcount(unsigned int x){return __builtin_popcount(x);
}

inline int popcountll(unsigned long long x){return __builtin_popcountll(x);
}

二进制数中末尾0的个数

__builtin_ctzll(0)的返回值为64

inline int ctzll(unsigned long long x){return __builtin_ctzll(x);
}

区间第k小函数

O(length)选出中点并将数组分成两边，使k左边的元素小等于等A[k]，使k右边的元素大于等于A[k]

 nth_element(A+l,A+k,A+r+1);

归并函数

 merge(A+1,A+n+1,B+1,B+m+1,C+1,cmp);

 inplace_merge(A+l,A+mid+1,A+r+1,cmp);

大随机数

rng()返回一个unsigned int范围内的随机数，搭配取模使用

mt19937 rng(chrono::steady_clock::now().time_since_epoch().count());

随机打乱数组元素顺序

mt19937 rng(chrono::steady_clock::now().time_since_epoch().count());shuffle(A+1,A+n+1,rng);

 srand(time(NULL));random_shuffle(A+1,A+n+1);//n大时打乱效果不佳

程序运行时间

inline double run_time(){return 1.0*clock()/CLOCKS_PER_SEC;
}

上/下一个排列

 prev_permutation(A+1,A+n+1);//上一个排列next_permutation(A+1,A+n+1);//下一个排列

    for(int i=1;i<=n;i++)A[i]=i;do{}while(next_permutation(A+1,A+n+1));//枚举所有排列

string转char

 string str="jiedai";puts(str.data());

数据结构

树状数组

一维树状数组

用法一：单点修改，区间查询

单点修改：add(bit,x,v)

区间查询：query(bit,r)-query(bit,l-1)

用法二：区间修改，单点查询

区间修改：add(bit,l,v),add(bit,r+1,-v)

单点查询：query(bit,x)

int limit;//树状数组下标的上界
void add(ll *bit,int x,ll v){while(x<=limit)bit[x]+=v,x+=x&-x;
}
ll query(ll *bit,int x){ll res=0;while(x)res+=bit[x],x-=x&-x;return res;
}

用法三：区间修改，区间查询

令数组 aaa 的差分数组为 ddd ，即 a[n]=∑i=1nd[i]a[n]=\sum\limits_{i=1}^{n}{d[i]}a[n]=i=1∑n​d[i]

则∑i=1ra[i]=∑i=1r∑j=1id[j]\sum\limits_{i=1}^{r}{a[i]}=\sum\limits_{i=1}^{r}{\sum\limits_{j=1}^{i}{d[j]}}i=1∑r​a[i]=i=1∑r​j=1∑i​d[j] =∑i=1r(r+1−i)∗d[i]=\sum\limits_{i=1}^{r}{(r+1-i)*d[i]}=i=1∑r​(r+1−i)∗d[i] =(r+1)∗∑i=1rd[i]−∑i=1rd[i]∗i=(r+1)*\sum\limits_{i=1}^{r}{d[i]}-\sum\limits_{i=1}^{r}{d[i]*i}=(r+1)∗i=1∑r​d[i]−i=1∑r​d[i]∗i

用 bit1bit1bit1 来维护 d[i]d[i]d[i] 的前缀和，用 bit2bit2bit2 来维护 d[i]∗id[i]*id[i]∗i 的前缀和

void real_add(ll *bit1,ll *bit2,int l,int r,ll v){add(bit1,l,v);add(bit2,l,v*l);add(bit1,r+1,-v);add(bit2,r+1,-v*(r+1));
}
ll real_query(ll *bit1,ll *bit2,int l,int r){return (r+1)*query(bit1,r)-query(bit2,r)-l*query(bit1,l-1)+query(bit2,l-1);
}

struct BIT{static const int SIZE=1e6+5;ll bit1[SIZE],bit2[SIZE];int limit;void init(int n=SIZE-1){limit=n;for(int i=1;i<=n;i++)bit1[i]=bit2[i]=0;}BIT(){init();}void add(ll *bit,int x,ll v){while(x<=limit)bit[x]+=v,x+=x&-x;}ll query(ll *bit,int x){ll res=0;while(x)res+=bit[x],x-=x&-x;return res;}void add(int l,int r,ll v){add(bit1,l,v);add(bit2,l,v*l);add(bit1,r+1,-v);add(bit2,r+1,-v*(r+1));}ll query(int l,int r){return (r+1)*query(bit1,r)-query(bit2,r)-l*query(bit1,l-1)+query(bit2,l-1);}
}bit;

用法四：求第一个前缀和大于等于某个值的下标

int sum_lower_bound(ll *bit,ll v){int pos=0;for(int i=20;i>=0;i--)if(pos+(1<<i)<=limit&&bit[pos+(1<<i)]<v)v-=bit[pos|=1<<i];return pos+1;
}

二维树状数组

用法一：单点修改，区间查询

单点修改：add(x,y,v)

区间查询：query(x2,y2)+query(x1-1,y1-1)-query(x1-1,y2)-query(x2,y1-1)

用法二：区间修改，单点查询

区间修改：add(x1,y1,v),add(x2+1,y2+1,v),add(x1,y2+1,-v),add(x2+1,y1,-v)

单点查询：query(x,y)

const int SIZE=5005;
int limit1=SIZE-1,limit2=SIZE-1;//两个维度的下标的上界
ll bit[SIZE][SIZE];void add(int x,int y,ll v){for(int i=x;i<=limit1;i+=i&-i)for(int j=y;j<=limit2;j+=j&-j)bit[i][j]+=v;
}
ll query(int x,int y){ll res=0;for(int i=x;i;i-=i&-i)for(int j=y;j;j-=j&-j)res+=bit[i][j];return res;
}

用法三：区间修改，区间查询

令数组 aaa 的二维差分数组为 ddd ，即 a[n][m]=∑i=1n∑j=1md[i][j]a[n][m]=\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m}{d[i][j]}a[n][m]=i=1∑n​j=1∑m​d[i][j]

则∑i=1x∑j=1ya[i][j]=∑i=1x∑j=1y∑k=1i∑t=1jd[k][t]\sum\limits_{i=1}^{x}\sum\limits_{j=1}^{y}{a[i][j]}=\sum\limits_{i=1}^{x}\sum\limits_{j=1}^{y}\sum\limits_{k=1}^{i}\sum\limits_{t=1}^{j}{d[k][t]}i=1∑x​j=1∑y​a[i][j]=i=1∑x​j=1∑y​k=1∑i​t=1∑j​d[k][t] =∑i=1x∑j=1y(x+1−i)(y+1−j)d[i][j]=\sum\limits_{i=1}^{x}\sum\limits_{j=1}^{y}{(x+1-i)(y+1-j)d[i][j]}=i=1∑x​j=1∑y​(x+1−i)(y+1−j)d[i][j]

=(x+1)(y+1)∑i=1x∑j=1yd[i][j]=(x+1)(y+1)\sum\limits_{i=1}^{x}\sum\limits_{j=1}^{y}{d[i][j]}=(x+1)(y+1)i=1∑x​j=1∑y​d[i][j] −(y+1)∑i=1x∑j=1yd[i][j]∗i-(y+1)\sum\limits_{i=1}^{x}\sum\limits_{j=1}^{y}{d[i][j]*i}−(y+1)i=1∑x​j=1∑y​d[i][j]∗i −(x+1)∑i=1x∑j=1yd[i][j]∗j-(x+1)\sum\limits_{i=1}^{x}\sum\limits_{j=1}^{y}{d[i][j]*j}−(x+1)i=1∑x​j=1∑y​d[i][j]∗j +∑i=1x∑j=1yd[i][j]∗ij+\sum\limits_{i=1}^{x}\sum\limits_{j=1}^{y}{d[i][j]*ij}+i=1∑x​j=1∑y​d[i][j]∗ij

用四个树状数组分别维护 d[i][j]d[i][j]d[i][j] , d[i][j]∗id[i][j]*id[i][j]∗i , d[i][j]∗jd[i][j]*jd[i][j]∗j , d[i][j]∗ijd[i][j]*ijd[i][j]∗ij 的前缀和

const int SIZE=5005;
int limit1=SIZE-1,limit2=SIZE-1;//两个维度的下标的上界
ll bit1[SIZE][SIZE],bit2[SIZE][SIZE],bit3[SIZE][SIZE],bit4[SIZE][SIZE];void add(int x,int y,ll v){for(int i=x;i<=limit1;i+=i&-i)for(int j=y;j<=limit2;j+=j&-j){bit1[i][j]+=v;bit2[i][j]+=v*x;bit3[i][j]+=v*y;bit4[i][j]+=v*x*y;}
}
ll query(int x,int y){ll res=0;for(int i=x;i;i-=i&-i)for(int j=y;j;j-=j&-j)res+=(x+1)*(y+1)*bit1[i][j]-(y+1)*bit2[i][j]-(x+1)*bit3[i][j]+bit4[i][j];return res;
}
void real_add(int x1,int y1,int x2,int y2,ll v){add(x1,y1,v);add(x1,y2+1,-v);add(x2+1,y1,-v);add(x2+1,y2+1,v);
}
ll real_query(int x1,int y1,int x2,int y2){return query(x2,y2)-query(x1-1,y2)-query(x2,y1-1)+query(x1-1,y1-1);
}

线段树

代码一：区间加减+区间求和

const int SIZE=1e6+5;
int LLL=1,RRR=SIZE-5;ll tree[SIZE<<2],lazy[SIZE<<2];
void build(int l=LLL,int r=RRR,int p=1){lazy[p]=0;if(l==r){tree[p]=0;return;}int mid=l+r>>1;build(l,mid,p<<1);build(mid+1,r,p<<1|1);tree[p]=tree[p<<1]+tree[p<<1|1];
}
void down(int p,int l,int r){if(lazy[p]){int mid=l+r>>1;lazy[p<<1]+=lazy[p];lazy[p<<1|1]+=lazy[p];tree[p<<1]+=1ll*(mid-l+1)*lazy[p];tree[p<<1|1]+=1ll*(r-mid)*lazy[p];lazy[p]=0;}
}
void update(int a,int b,ll v,int l=LLL,int r=RRR,int p=1){if(l>b||r<a)return;if(l>=a&&r<=b){tree[p]+=1ll*(r-l+1)*v;lazy[p]+=v;return;}down(p,l,r);int mid=l+r>>1;update(a,b,v,l,mid,p<<1);update(a,b,v,mid+1,r,p<<1|1);tree[p]=tree[p<<1]+tree[p<<1|1];
}
ll query(int a,int b,int l=LLL,int r=RRR,int p=1){if(l>b||r<a)return 0;if(l>=a&&r<=b)return tree[p];down(p,l,r);int mid=l+r>>1;return query(a,b,l,mid,p<<1)+query(a,b,mid+1,r,p<<1|1);
}

代码二：区间加减+区间最大值

const int SIZE=1e6+5;
int LLL=1,RRR=SIZE-5;ll tree[SIZE<<2],lazy[SIZE<<2];
void build(int l=LLL,int r=RRR,int p=1){lazy[p]=0;if(l==r){tree[p]=0;return;}int mid=l+r>>1;build(l,mid,p<<1);build(mid+1,r,p<<1|1);tree[p]=max(tree[p<<1],tree[p<<1|1]);
}
void down(int p){if(lazy[p]){tree[p<<1]+=lazy[p];lazy[p<<1]+=lazy[p];tree[p<<1|1]+=lazy[p];lazy[p<<1|1]+=lazy[p];lazy[p]=0;}
}
void update(int a,int b,ll v,int l=LLL,int r=RRR,int p=1){if(l>b||r<a)return;if(l>=a&&r<=b){tree[p]+=v;lazy[p]+=v;return;}down(p);int mid=l+r>>1;update(a,b,v,l,mid,p<<1);update(a,b,v,mid+1,r,p<<1|1);tree[p]=max(tree[p<<1],tree[p<<1|1]);
}
ll query(int a,int b,int l=LLL,int r=RRR,int p=1){if(l>b||r<a)return 0;//if(l>=a&&r<=b)return tree[p];down(p);int mid=l+r>>1;return max(query(a,b,l,mid,p<<1),query(a,b,mid+1,r,p<<1|1));
}

代码三：查询区间第一个大于等于某个值的位置

int query_pos(int a,int b,int v,int l=LLL,int r=RRR,int p=1){if(l>b||r<a)return 0;if(l==r)return l;down(p);int mid=l+r>>1,res=0;if(tree[p<<1]>=v)res=query_pos(a,b,v,l,mid,p<<1);if(res==0&&tree[p<<1|1]>=v)res=query_pos(a,b,v,mid+1,r,p<<1|1);return res;
}

代码四：P3373 区间加、区间乘、区间求和

const int M=1e5+5;
int n,m,mod,A[M],tree[M<<2],lazy_mul[M<<2],lazy_add[M<<2];
void build(int l=1,int r=n,int p=1){lazy_mul[p]=1;lazy_add[p]=0;if(l==r){tree[p]=A[l]%mod;return;}int mid=l+r>>1;build(l,mid,p<<1);build(mid+1,r,p<<1|1);tree[p]=(tree[p<<1]+tree[p<<1|1])%mod;
}
void down(int p,int l,int r){if(lazy_mul[p]!=1){tree[p<<1]=1ll*tree[p<<1]*lazy_mul[p]%mod;lazy_mul[p<<1]=1ll*lazy_mul[p<<1]*lazy_mul[p]%mod;lazy_add[p<<1]=1ll*lazy_add[p<<1]*lazy_mul[p]%mod;tree[p<<1|1]=1ll*tree[p<<1|1]*lazy_mul[p]%mod;lazy_mul[p<<1|1]=1ll*lazy_mul[p<<1|1]*lazy_mul[p]%mod;lazy_add[p<<1|1]=1ll*lazy_add[p<<1|1]*lazy_mul[p]%mod;lazy_mul[p]=1;}if(lazy_add[p]){int mid=l+r>>1;tree[p<<1]=(tree[p<<1]+1ll*(mid-l+1)*lazy_add[p])%mod;lazy_add[p<<1]=(lazy_add[p<<1]+lazy_add[p])%mod;tree[p<<1|1]=(tree[p<<1|1]+1ll*(r-mid)*lazy_add[p])%mod;lazy_add[p<<1|1]=(lazy_add[p<<1|1]+lazy_add[p])%mod;lazy_add[p]=0;}
}
void update(int a,int b,int mul,int add,int l=1,int r=n,int p=1){if(l>b||r<a)return;if(l>=a&&r<=b){if(mul!=1){tree[p]=1ll*tree[p]*mul%mod;lazy_mul[p]=1ll*lazy_mul[p]*mul%mod;lazy_add[p]=1ll*lazy_add[p]*mul%mod;}if(add){tree[p]=(tree[p]+1ll*(r-l+1)*add)%mod;lazy_add[p]=(lazy_add[p]+add)%mod;}return;}down(p,l,r);int mid=l+r>>1;update(a,b,mul,add,l,mid,p<<1);update(a,b,mul,add,mid+1,r,p<<1|1);tree[p]=(tree[p<<1]+tree[p<<1|1])%mod;
}
int query(int a,int b,int l=1,int r=n,int p=1){if(l>b||r<a)return 0;if(l>=a&&r<=b)return tree[p];down(p,l,r);int mid=l+r>>1;return (query(a,b,l,mid,p<<1)+query(a,b,mid+1,r,p<<1|1))%mod;
}

主席树

静态区间第k小 P3834

const int SIZE=2e5+5;
const int TREE=SIZE*20;
int n,q,uni,A[SIZE],tmp[SIZE];
int trid,rt[SIZE],ls[TREE],rs[TREE],sum[TREE];
void update(int &now,int lst,int x,int l=1,int r=uni){now=++trid;ls[now]=ls[lst],rs[now]=rs[lst],sum[now]=sum[lst]+1;if(l==r)return;int mid=l+r>>1;if(x<=mid)update(ls[now],ls[lst],x,l,mid);else update(rs[now],rs[lst],x,mid+1,r);
}
int query(int left,int right,int k,int l=1,int r=uni){if(l==r)return l;int mid=l+r>>1,res=sum[ls[right]]-sum[ls[left]];if(k<=res)return query(ls[left],ls[right],k,l,mid);else return query(rs[left],rs[right],k-res,mid+1,r);
}
int main(){#ifndef ONLINE_JUDGEfreopen("jiedai.in","r",stdin);
#endifrd(n),rd(q);for(int i=1;i<=n;i++)rd(A[i]),tmp[++uni]=A[i];sort(tmp+1,tmp+1+uni);uni=unique(tmp+1,tmp+1+uni)-tmp-1;for(int i=1;i<=n;i++)A[i]=lower_bound(tmp+1,tmp+1+uni,A[i])-tmp;for(int i=1;i<=n;i++)update(rt[i],rt[i-1],A[i]);while(q--){int l,r,k;rd(l),rd(r),rd(k);printf("%d\n",tmp[query(rt[l-1],rt[r],k)]);}return (0-0);
}

线段树合并/裂开

线段树合并的复杂度证明：合并两个线段树的复杂度是它们的交集部分，合并时每到达一个点都意味着有两个点变成了一个点，也就是少了一个点。因此线段树无论怎么合并，复杂度都不会超过插入时的新建节点总数，即 O(nlogn)O(nlogn)O(nlogn)

P5494

维护若干可重集合，初始集合数为1，支持以下操作：

0 p x y ：将可重集 p 中大于等于 x 且小于等于 y 的值放入一个新的可重集中

1 p t ：将可重集 t 中的数并入可重集 p，且删除可重集 t

2 p x q ：在可重集 p 中加入 x 个数字 q

3 p x y ：查询可重集 p 中大于等于 x 且小于等于 y 的值的个数

4 p k ：查询可重集 p 中第 k 小的数，不存在输出 -1

const int SIZE=2e5+5;
const int TREE=SIZE*32;
int LLL=1,RRR=SIZE-5;int trid,rub,rt[SIZE],bin[TREE],ls[TREE],rs[TREE];
ll cnt[TREE];
int new_node(){if(rub)return bin[rub--];else return ++trid;
}
void del_node(int p){ls[p]=rs[p]=cnt[p]=0;bin[++rub]=p;
}
void update(int &p,int x,int v,int l=LLL,int r=RRR){if(p==0)p=new_node();cnt[p]+=v;if(l==r)return;int mid=l+r>>1;if(x<=mid)update(ls[p],x,v,l,mid);else update(rs[p],x,v,mid+1,r);
}
ll query_cnt(int p,int a,int b,int l=LLL,int r=RRR){if(p==0||l>b||r<a)return 0;if(l>=a&&r<=b)return cnt[p];int mid=l+r>>1;return query_cnt(ls[p],a,b,l,mid)+query_cnt(rs[p],a,b,mid+1,r);
}
int query_kth(int p,ll k,int l=LLL,int r=RRR){if(l==r)return l;int mid=l+r>>1;if(k<=cnt[ls[p]])return query_kth(ls[p],k,l,mid);else return query_kth(rs[p],k-cnt[ls[p]],mid+1,r);
}
int merge(int x,int y){if(!x||!y)return x|y;cnt[x]+=cnt[y];ls[x]=merge(ls[x],ls[y]);rs[x]=merge(rs[x],rs[y]);del_node(y);return x;
}
void split(int x,int &y,ll k){if(x==0)return;y=new_node();ll v=cnt[ls[x]];if(k<=v)swap(rs[x],rs[y]);if(k<v)split(ls[x],ls[y],k);if(k>v)split(rs[x],rs[y],k-v);cnt[y]=cnt[x]-k;cnt[x]=k;
}int n,q,all;
int main(){#ifndef ONLINE_JUDGEfreopen("jiedai.in","r",stdin);
#endifrd(n),rd(q);rt[all=1]=new_node();for(int i=1;i<=n;i++){int x;rd(x);update(rt[1],i,x);}while(q--){int op,p,x,y;ll k;rd(op);if(op==0){rd(p),rd(x),rd(y);ll k1=query_cnt(rt[p],1,y);ll k2=query_cnt(rt[p],x,y);int tmp;split(rt[p],rt[++all],k1-k2);split(rt[all],tmp,k2);rt[p]=merge(rt[p],tmp);}else if(op==1){rd(p),rd(y);rt[p]=merge(rt[p],rt[y]);}else if(op==2){rd(p),rd(x),rd(y);update(rt[p],y,x);}else if(op==3){rd(p),rd(x),rd(y);printf("%lld\n",query_cnt(rt[p],x,y));}else if(op==4){rd(p),rd(k);if(k<=0||k>query_cnt(rt[p],1,n))puts("-1");else printf("%d\n",query_kth(rt[p],k));}}return (0-0);
}

吉司机线段树

Segment Tree Beats!

BZOJ4695

维护一个序列，支持下面几种操作

1 L R v ：给区间[L,R]加上一个数v

2 L R v ：把区间[L,R]中小于v的数变成v

3 L R v ：把区间[L,R]中大于v的数变成v

4 L R ：求区间[L,R]的和

5 L R ：求区间[L,R]的最大值

6 L R ：求区间[L,R]的最小值

时间复杂度 O(nlog2n)O(nlog^2{n})O(nlog2n)

const int SIZE=5e5+5;
const ll INF=1e18;int n;
ll A[SIZE];int len[SIZE<<2],cmx[SIZE<<2],cmi[SIZE<<2];
ll sum[SIZE<<2],mx[SIZE<<2],mx2[SIZE<<2],mi[SIZE<<2],mi2[SIZE<<2];
ll lazy_mx[SIZE<<2],lazy_mi[SIZE<<2],lazy_add[SIZE<<2];
void up(int p){int ls=(p<<1),rs=(p<<1|1);sum[p]=sum[ls]+sum[rs];if(mx[ls]>mx[rs]){mx[p]=mx[ls];mx2[p]=max(mx2[ls],mx[rs]);cmx[p]=cmx[ls];}else if(mx[ls]<mx[rs]){mx[p]=mx[rs];mx2[p]=max(mx[ls],mx2[rs]);cmx[p]=cmx[rs];}else{mx[p]=mx[ls];mx2[p]=max(mx2[ls],mx2[rs]);cmx[p]=cmx[ls]+cmx[rs];}if(mi[ls]<mi[rs]){mi[p]=mi[ls];mi2[p]=min(mi2[ls],mi[rs]);cmi[p]=cmi[ls];}else if(mi[ls]>mi[rs]){mi[p]=mi[rs];mi2[p]=min(mi[ls],mi2[rs]);cmi[p]=cmi[rs];}else{mi[p]=mi[ls];mi2[p]=min(mi2[ls],mi2[rs]);cmi[p]=cmi[ls]+cmi[rs];}
}
void work_max(int p,ll v){lazy_mx[p]+=v;if(mx[p]==mi[p])mi[p]+=v;else if(cmx[p]+cmi[p]==len[p])mi2[p]+=v;mx[p]+=v;sum[p]+=1ll*v*cmx[p];
}
void work_min(int p,ll v){lazy_mi[p]+=v;if(mi[p]==mx[p])mx[p]+=v;else if(cmi[p]+cmx[p]==len[p])mx2[p]+=v;mi[p]+=v;sum[p]+=1ll*v*cmi[p];
}
void work_add(int p,ll v){lazy_add[p]+=v;mx[p]+=v;mi[p]+=v;if(mx2[p]!=-INF)mx2[p]+=v;if(mi2[p]!=INF)mi2[p]+=v;sum[p]+=1ll*v*len[p];
}
void down(int p){int ls=(p<<1),rs=(p<<1|1);if(lazy_mx[p]){if(mx[ls]>mx[rs])work_max(ls,lazy_mx[p]);else if(mx[ls]<mx[rs])work_max(rs,lazy_mx[p]);else work_max(ls,lazy_mx[p]),work_max(rs,lazy_mx[p]);lazy_mx[p]=0;}if(lazy_mi[p]){if(mi[ls]<mi[rs])work_min(ls,lazy_mi[p]);else if(mi[ls]>mi[rs])work_min(rs,lazy_mi[p]);else work_min(ls,lazy_mi[p]),work_min(rs,lazy_mi[p]);lazy_mi[p]=0;}if(lazy_add[p]){work_add(ls,lazy_add[p]);work_add(rs,lazy_add[p]);lazy_add[p]=0;}
}
void build(int l=1,int r=n,int p=1){len[p]=r-l+1;lazy_mx[p]=lazy_mi[p]=lazy_add[p]=0;if(l==r){mx[p]=mi[p]=sum[p]=A[l];mx2[p]=-INF,mi2[p]=INF;cmx[p]=cmi[p]=1;return;}int mid=l+r>>1;build(l,mid,p<<1);build(mid+1,r,p<<1|1);up(p);
}
void update_max(int a,int b,ll v,int l=1,int r=n,int p=1){if(l>b||r<a)return;if(l>=a&&r<=b){if(v<=mi[p])return;else if(v<mi2[p]){work_min(p,v-mi[p]);return;}}down(p);int mid=l+r>>1;update_max(a,b,v,l,mid,p<<1);update_max(a,b,v,mid+1,r,p<<1|1);up(p);
}
void update_min(int a,int b,ll v,int l=1,int r=n,int p=1){if(l>b||r<a)return;if(l>=a&&r<=b){if(v>=mx[p])return;else if(v>mx2[p]){work_max(p,v-mx[p]);return;}}down(p);int mid=l+r>>1;update_min(a,b,v,l,mid,p<<1);update_min(a,b,v,mid+1,r,p<<1|1);up(p);
}
void update_add(int a,int b,ll v,int l=1,int r=n,int p=1){if(l>b||r<a)return;if(l>=a&&r<=b){work_add(p,v);return;}down(p);int mid=l+r>>1;update_add(a,b,v,l,mid,p<<1);update_add(a,b,v,mid+1,r,p<<1|1);up(p);
}
ll query_max(int a,int b,int l=1,int r=n,int p=1){if(l>b||r<a)return -INF;if(l>=a&&r<=b)return mx[p];down(p);int mid=l+r>>1;return max(query_max(a,b,l,mid,p<<1),query_max(a,b,mid+1,r,p<<1|1));
}
ll query_min(int a,int b,int l=1,int r=n,int p=1){if(l>b||r<a)return INF;if(l>=a&&r<=b)return mi[p];down(p);int mid=l+r>>1;return min(query_min(a,b,l,mid,p<<1),query_min(a,b,mid+1,r,p<<1|1));
}
ll query_sum(int a,int b,int l=1,int r=n,int p=1){if(l>b||r<a)return 0;if(l>=a&&r<=b)return sum[p];down(p);int mid=l+r>>1;return query_sum(a,b,l,mid,p<<1)+query_sum(a,b,mid+1,r,p<<1|1);
}
int main(){#ifndef ONLINE_JUDGEfreopen("jiedai.in","r",stdin);
#endifint q,op,l,r;ll v;rd(n);for(int i=1;i<=n;i++)rd(A[i]);build();rd(q);while(q--){rd(op),rd(l),rd(r);if(0-0);else if(op==1)rd(v),update_add(l,r,v);else if(op==2)rd(v),update_max(l,r,v);else if(op==3)rd(v),update_min(l,r,v);else if(op==4)printf("%lld\n",query_sum(l,r));else if(op==5)printf("%lld\n",query_max(l,r));else if(op==6)printf("%lld\n",query_min(l,r));}return (0-0);
}

李超线段树*

可撤销并查集

struct disjoint_set{static const int SIZE=1e6+5;int top,stk[SIZE],fa[SIZE],sz[SIZE];inline void init(int n){top=0;for(int i=1;i<=n;i++)fa[i]=i,sz[i]=1;}inline int getfa(int x){while(x!=fa[x])x=fa[x];return x;}inline bool link(int x,int y){x=getfa(x),y=getfa(y);if(x!=y){if(sz[x]>sz[y])swap(x,y);stk[++top]=x;fa[x]=y;sz[y]+=sz[x];return 1;}else return 0;}inline void pop_to(int tar){while(top>tar){sz[fa[stk[top]]]-=sz[stk[top]];fa[stk[top]]=stk[top];top--;}}
}U;

ST表

静态区间最大值 P3865

const int M=1e5+5;
int n,q,A[M],st[M][20];
inline int ST(int l,int r){int k=log2(r-l+1);return max(st[l][k],st[r-(1<<k)+1][k]);
}
int main(){#ifndef ONLINE_JUDGEfreopen("jiedai.in","r",stdin);
#endifrd(n),rd(q);for(int i=1;i<=n;i++)rd(A[i]);for(int i=1;i<=n;i++)st[i][0]=A[i];for(int i=n;i>=1;i--)for(int j=1;i+(1<<j)-1<=n;j++)st[i][j]=max(st[i][j-1],st[i+(1<<j-1)][j-1]);while(q--){int l,r;rd(l),rd(r);printf("%d\n",ST(l,r));}return (0-0);
}

分块

代码一： O((n))O(\sqrt(n))O((​n))单点修改 + O(1)O(1)O(1)区间查询

struct BLOCK1{static const int M=2e5+5;static const int S=400;static const int K=M/S+5;int bel[M],L[K],R[K];ll big[K],small[M];BLOCK1(){for(int i=1;i<M;i++)bel[i]=(i-1)/S+1;for(int i=1;i<K;i++)L[i]=1+(i-1)*S,R[i]=i*S;}void update(int x,ll v){for(int i=1;i<=bel[x];i++)big[i]+=v;for(int i=L[bel[x]];i<=x;i++)small[i]+=v;}ll query(int l,int r){return big[bel[l]+1]+small[l]-big[bel[r+1]+1]-small[r+1];}
}block1;

代码二： O(1)O(1)O(1)单点修改 + O((n))O(\sqrt(n))O((​n))区间查询

struct BLOCK2{static const int M=2e5+5;static const int S=400;static const int K=M/S+5;int bel[M],L[K],R[K];ll big[K],small[M];BLOCK2(){for(int i=1;i<M;i++)bel[i]=(i-1)/S+1;for(int i=1;i<K;i++)L[i]=1+(i-1)*S,R[i]=i*S;}void update(int x,ll v){big[bel[x]]+=v;small[x]+=v;}ll query(int l,int r){ll res=0;if(bel[l]!=bel[r]){for(int i=l;i<=R[bel[l]];i++)res+=small[i];for(int i=bel[l]+1;i<bel[r];i++)res+=big[i];for(int i=L[bel[r]];i<=r;i++)res+=small[i];}else for(int i=l;i<=r;i++)res+=small[i];return res;}
}block2;

代码三： O((n))O(\sqrt(n))O((​n))区间修改 + O(1)O(1)O(1)单点查询

struct BLOCK3{static const int M=2e5+5;static const int S=400;static const int K=M/S+5;int bel[M],L[K],R[K];ll big[K],small[M];BLOCK3(){for(int i=1;i<M;i++)bel[i]=(i-1)/S+1;for(int i=1;i<K;i++)L[i]=1+(i-1)*S,R[i]=i*S;}void update(int l,int r,ll v){if(bel[l]!=bel[r]){for(int i=l;i<=R[bel[l]];i++)small[i]+=v;for(int i=bel[l]+1;i<bel[r];i++)big[i]+=v;for(int i=L[bel[r]];i<=r;i++)small[i]+=v;}else for(int i=l;i<=r;i++)small[i]+=v;}ll query(int x){return big[bel[x]]+small[x];}
}block3;

代码四： O(1)O(1)O(1)区间修改 + O((n))O(\sqrt(n))O((​n))单点查询

struct BLOCK4{static const int M=2e5+5;static const int S=400;static const int K=M/S+5;int bel[M],L[K],R[K];ll big[K],small[M];BLOCK4(){for(int i=1;i<M;i++)bel[i]=(i-1)/S+1;for(int i=1;i<K;i++)L[i]=1+(i-1)*S,R[i]=i*S;}void update(int l,int r,ll v){big[bel[l]]+=v,small[l]+=v;big[bel[r+1]]-=v,small[r+1]-=v;}ll query(int x){ll res=0;for(int i=0;i<bel[x];i++)res+=big[i];for(int i=L[bel[x]];i<=x;i++)res+=small[i];return res;}
}block4;

莫队

普通莫队

给定长度为n的序列，m次询问，求区间 [l,r][l,r][l,r] 之间有多少种不同的数字 P1972

const int M=1e6+5;
const int S=sqrt(M);
struct node{int l,r,id;bool operator <(const node &A)const{if(l/S!=A.l/S)return l/S<A.l/S;if(l/S&1)return r<A.r;else return r>A.r;}
}Q[M];
int n,m,res,A[M],ans[M],cnt[M];
void insert(int x){cnt[x]++;if(cnt[x]==1)res++;
}
void erase(int x){cnt[x]--;if(cnt[x]==0)res--;
}
int main(){#ifndef ONLINE_JUDGEfreopen("jiedai.in","r",stdin);
#endifrd(n);for(int i=1;i<=n;i++)rd(A[i]);rd(m);for(int i=1;i<=m;i++)rd(Q[i].l),rd(Q[i].r),Q[i].id=i;sort(Q+1,Q+1+m);int l=1,r=0;for(int i=1;i<=m;i++){while(r<Q[i].r)insert(A[++r]);while(l>Q[i].l)insert(A[--l]);while(r>Q[i].r)erase(A[r--]);while(l<Q[i].l)erase(A[l++]);ans[Q[i].id]=res;}for(int i=1;i<=m;i++)printf("%d\n",ans[i]);return (0-0);
}

带修莫队

1. 修改某个位置的数字

2. 求区间 [l,r][l,r][l,r] 之间有多少种不同的数字

这是一道带修莫队的模板题。

每个查询操作除了有左端点l和右端点r，还有一个时间t，表示这个查询在第t个修改操作之后。

将这些查询(l,r,t)进行排序，以l/S为第一关键字，r/S为第二关键字，t为第三关键字。其中S为分块的大小。

与普通莫队相比，当从上一个询问跳到当前的询问时，还要维护目前的修改数now。

当now小于查询的修改数时，就把没有进行的修改操作执行。

当now大于查询的修改数时，就将多出的已执行的修改进行回撤。

需要注意的是，如果执行或回撤的修改操作在目前的区间内，就要更新答案。

在修改数now与查询的修改数相等之后，再进行左右端点的移动。

复杂度分析：n为序列长度，c为修改操作数，q为查询操作数，S为分块的大小。

1. 修改数的移动：

左端点分成n/S块，右端点分成n/S块，每个左右端点块最多移动c。因此移动次数为O(c∗n2/S2)O(c*n^2/S^2)O(c∗n2/S2)

1. 左端点的移动：

在同一左端点块中，每次最多移动S；到下一左端点块时，最多移动2S。因此移动次数为O(qS)O(qS)O(qS)

1. 右端点的移动：

在同一右端点块中，每次最多移动S；到下一右端点块时，每次最多移动2S；

当左端点改变时，每次最多移动n。因此移动次数为O(qS+n2/S)O(qS+n^2/S)O(qS+n2/S)

于是总移动次数为O(qS+c∗n2/S2+n2/S)O(qS+c*n^2/S^2+n^2/S)O(qS+c∗n2/S2+n2/S)。

由于不知道c和q的值，那么c和q统一用m代替，且假设n=m，那么总移动次数就成了O(qS+n3/S2+n2/S)O(qS+n^3/S^2+n^2/S)O(qS+n3/S2+n2/S)。

当S=n2/3S=n^{2/3}S=n2/3时，最优复杂度为n5/3n^{5/3}n5/3。

const int M=140005;
const int S=2500;
const int K=1e6+5;
struct node{int l,r,t,id;bool operator <(const node &A)const{if(l/S!=A.l/S)return l/S<A.l/S;if(r/S!=A.r/S){if(l/S&1)return r/S<A.r/S;else return r/S>A.r/S;}if(r/S&1)return t<A.t;else return t>A.t;}
}Q[M];
char str[M][5];
int n,m,A[M],pos[M],X[M],Y[M],ans[M];
int res,cnt[K];
void insert(int x){cnt[x]++;if(cnt[x]==1)res++;
}
void erase(int x){cnt[x]--;if(cnt[x]==0)res--;
}
int main(){#ifndef ONLINE_JUDGEfreopen("jiedai.in","r",stdin);
#endifrd(n),rd(m);for(int i=1;i<=n;i++)rd(A[i]);int l=1,r=0,now=0,tot=0;for(int i=1;i<=m;i++){scanf("%s",str[i]);int a,b;rd(a),rd(b);if(str[i][0]=='R'){now++;pos[now]=a;X[now]=A[a];Y[now]=b;A[a]=b;}else Q[++tot]=(node){a,b,now,i};}sort(Q+1,Q+1+tot);for(int i=1;i<=tot;i++){while(now<Q[i].t){now++;A[pos[now]]=Y[now];if(pos[now]>=l&&pos[now]<=r)erase(X[now]),insert(Y[now]);}while(now>Q[i].t){A[pos[now]]=X[now];if(pos[now]>=l&&pos[now]<=r)erase(Y[now]),insert(X[now]);now--;}while(r<Q[i].r)insert(A[++r]);while(l>Q[i].l)insert(A[--l]);while(r>Q[i].r)erase(A[r--]);while(l<Q[i].l)erase(A[l++]);ans[Q[i].id]=res;}for(int i=1;i<=m;i++)if(str[i][0]=='Q')printf("%d\n",ans[i]);return (0-0);
}

树上带修莫队

nnn个节点的树，mmm种颜色，qqq个询问。
树上每个节点有一个颜色，每种颜色有一个权值。
第iii次遍历到某种颜色的点，该点的价值乘上对应的数值。
000 xxx yyy : 把xxx点的颜色改为yyy
111 xxx yyy : 查询遍历xxx到yyy路径的总价值。

思路：用树的括号序将查询的一条链转换成括号序列上的一个区间，再在括号序列上利用带修莫队求解。

将链转换成区间

处理出树的括号序，在对树进行 dfsdfsdfs 的时候，进入和离开某个点的时候都要打上标记，记xxx点进入和离开的时间戳为 idL[x]idL[x]idL[x] 和 idR[x]idR[x]idR[x]。
将一条链转换成一个括号区间，要保证链上的点在区间中只出现一个括号，不在链上的点在区间上没有出现，或者出现的两个括号相互抵消。

对于一条端点为 xxx 和 yyy 的链，假设 idL[x]≤idL[y]idL[x]≤idL[y]idL[x]≤idL[y] ，算出 x 和 y 的 LCALCALCA。

1. 如果 xxx 或 yyy 就是 lcalcalca ，那么这条链对应的括号序区间为idL[x]idL[x]idL[x]到idL[y]idL[y]idL[y]。
2. 如果 xxx 和 yyy 都不是 lcalcalca ，那么对应的区间为 idR[x]idR[x]idR[x] 到 idL[y]idL[y]idL[y]。不过这个区间中并没有包含 lcalcalca 这个点，这需要在处理询问的时候特别加上。

带修莫队

莫队的区间中第二次加入一个点时，改为删除该点，撤销该操作就是加入该点，可以标记该点的奇偶性来进行判断。
加上修改操作后同理，利用奇偶性判断是否修改。

const int M=1e5+5;
const int S=2000;
int n,m,q,A[M],B[M],C[M],tot,head[M],to[M<<1],nxt[M<<1];
inline void add_edge(int a,int b){to[++tot]=b;nxt[tot]=head[a];head[a]=tot;
}
int fa[M],son[M],sz[M],deep[M],top[M];
int reid[M<<1],idL[M],idR[M],dfsid;
void dfs(int x,int f){fa[x]=f;sz[x]=1;deep[x]=deep[f]+1;reid[idL[x]=++dfsid]=x;for(int i=head[x];i;i=nxt[i]){int y=to[i];if(y==f)continue;dfs(y,x);sz[x]+=sz[y];if(sz[y]>sz[son[x]])son[x]=y;}reid[idR[x]=++dfsid]=x;
}
void chain(int x,int k){top[x]=k;if(son[x])chain(son[x],k);for(int i=head[x];i;i=nxt[i]){int y=to[i];if(y==fa[x]||y==son[x])continue;chain(y,y);}
}
int LCA(int a,int b){while(top[a]!=top[b]){if(deep[top[a]]<deep[top[b]])swap(a,b);a=fa[top[a]];}return deep[a]<deep[b]?a:b;
}
struct node{int l,r,t,lca,id;bool operator <(const node &A)const{if(l/S!=A.l/S)return l/S<A.l/S;if(r/S!=A.r/S){if(l/S%2)return r/S<A.r/S;else return r/S>A.r/S;}if(r/S%2)return t<A.t;else return t>A.t;}
}Q[M];
int op0,op1,X[M],Y[M],Z[M],mark[M],cnt[M];
ll res,ans[M];
void add(int x){mark[x]^=1;if(mark[x])res+=1ll*A[C[x]]*B[++cnt[C[x]]];else res-=1ll*A[C[x]]*B[cnt[C[x]]--];
}
int main(){#ifndef ONLINE_JUDGEfreopen("jiedai.in","r",stdin);
#endifrd(n),rd(m),rd(q);for(int i=1;i<=m;i++)rd(A[i]);for(int i=1;i<=n;i++)rd(B[i]);for(int i=1;i<n;i++){int a,b;rd(a),rd(b);add_edge(a,b),add_edge(b,a);}for(int i=1;i<=n;i++)rd(C[i]);dfs(1,0);chain(1,1);int now=0;for(int i=1;i<=q;i++){int op,x,y;rd(op),rd(x),rd(y);if(op==0){X[++op0]=x,Y[op0]=y,Z[op0]=C[x];C[x]=y,now++;}else{int lca=LCA(x,y);if(idL[x]>idL[y])swap(x,y);if(x==lca||y==lca)Q[++op1]=(node){idL[x],idL[y],op0,0,i};else Q[++op1]=(node){idR[x],idL[y],op0,lca,i};}}sort(Q+1,Q+1+op1);int L=Q[1].l,R=L-1;for(int i=1;i<=op1;i++){int l=Q[i].l,r=Q[i].r,t=Q[i].t,lca=Q[i].lca,id=Q[i].id;while(L>l)add(reid[--L]);while(R<r)add(reid[++R]);while(L<l)add(reid[L++]);while(R>r)add(reid[R--]);while(now<t){now++;if(mark[X[now]]){add(X[now]);C[X[now]]=Y[now];add(X[now]);}else C[X[now]]=Y[now];}while(now>t){if(mark[X[now]]){add(X[now]);C[X[now]]=Z[now];add(X[now]);}else C[X[now]]=Z[now];now--;}if(lca)add(lca);ans[id]=res;if(lca)add(lca);}for(int i=1;i<=q;i++)if(ans[i])printf("%lld\n",ans[i]);return (0-0);
}

回滚莫队

loj 2874

求[l,r][l,r][l,r]中一种数权值与出现次数之积的最大值

const int M=2e5+5;
const int S=350;
struct node{int l,r,id;bool operator <(const node &A)const{if(l/S!=A.l/S)return l/S<A.l/S;return r<A.r;}
}Q[M];
int n,m,A[M],B[M],cnt[M],cnt2[M];
ll res,ans[M];int main(){#ifndef ONLINE_JUDGEfreopen("jiedai.in","r",stdin);
#endifrd(n),rd(m);for(int i=1;i<=n;i++)rd(A[i]),B[i]=A[i];sort(B+1,B+1+n);int uni=unique(B+1,B+1+n)-B-1;for(int i=1;i<=n;i++)A[i]=lower_bound(B+1,B+1+uni,A[i])-B;for(int i=1;i<=m;i++){int l,r;rd(l),rd(r);Q[i]=(node){l,r,i};}sort(Q+1,Q+1+m);int L=0,R=0;for(int i=1;i<=m;i++){if(i==1||Q[i].l/S!=Q[i-1].l/S){L=min(n+1,Q[i].l/S*S+S);R=L-1;res=0;for(int j=1;j<=uni;j++)cnt[j]=0;}int l=Q[i].l,r=Q[i].r,id=Q[i].id;if(l/S==r/S){ll tmp=0;for(int j=l;j<=r;j++)MAX(tmp,1ll*B[A[j]]*(++cnt2[A[j]]));for(int j=l;j<=r;j++)cnt2[A[j]]=0;ans[id]=tmp;}else{while(R<r)++R,MAX(res,1ll*B[A[R]]*(++cnt[A[R]]));ll tmp=res;for(int j=l;j<L;j++)MAX(tmp,1ll*B[A[j]]*(++cnt[A[j]]));for(int j=l;j<L;j++)cnt[A[j]]--;ans[id]=tmp;}}for(int i=1;i<=m;i++)printf("%lld\n",ans[i]);return (0-0);
}

二次离线莫队

二次离线莫队算法可以计算区间中满足某种条件的点对的数量，例如P5047

题意：给定一个长度为n的序列 m次查询 求区间[L,R]中逆序对的数量

分析：如果直接用普通莫队做这道题的话，每次移动都要在树状数组上修改和查询，复杂度为 O(nnlogn)O(n\sqrt{n}logn)O(nn​logn)

记 ([l1,r1],[l2,r2])([l1,r1],[l2,r2])([l1,r1],[l2,r2]) 表示一个点在 [l1,r1][l1,r1][l1,r1] 中另外一点在 [l2,r2][l2,r2][l2,r2] 中的逆序对数

在做莫队的过程中，假如要将区间 [l,r][l,r][l,r] ，通过右端点右移，变成区间 [l,r′][l,r'][l,r′] 。那么一次从 [l,i−1][l,i-1][l,i−1] 变到 [l,i][l,i][l,i] 的移动，要让当前维护的答案加上 ([l,i−1],[i,i])([l,i-1],[i,i])([l,i−1],[i,i]) ，转化一下也就是 ([1,i−1],[i,i])−([1,l−1],[i,i])([1,i-1],[i,i])-([1,l-1],[i,i])([1,i−1],[i,i])−([1,l−1],[i,i])

因此在所有移动之后，当前答案应该加上：∑i=r+1r′([1,i−1],[i,i])−([1,l−1],[r+1,r′])\sum\limits_{i=r+1}^{r'}{([1,i-1],[i,i])}-([1,l-1],[r+1,r'])i=r+1∑r′​([1,i−1],[i,i])−([1,l−1],[r+1,r′])

前者只有n种情况，用树状数组 O(nlogn)O(nlogn)O(nlogn) 进行预处理，然后利用前缀和进行计算

后者是一个前缀和一个区间的形式，我们可以离线将每个区间保存在其前缀上，那么可以从左到右处理每一个前缀以及挂在这个前缀上的区间 [r+1,r′][r+1,r'][r+1,r′]

所有 [r+1,r′][r+1,r'][r+1,r′] 的总长，由于莫队的性质，是 O(nn)O(n\sqrt{n})O(nn​) 级别的，而前缀的数量是 O(n)O(n)O(n) 的，我们可以用一个 O(n)O(\sqrt{n})O(n​) 单点修改 O(1)O(1)O(1) 查询的分块。每次前缀加入一个新数，就在容器中 O(n)O(\sqrt{n})O(n​) 进行修改。查询一个位置的元素与前缀形成的逆序对数是 O(1)O(1)O(1) 的

总复杂度就是 O(nlogn+nn)O(nlogn+n\sqrt{n})O(nlogn+nn​)

对于左端点左移，左端点右移，右端点左移也类似，要注意的是左端点移动要用后缀

const int M=1e5+5;
const int S=sqrt(M);
const int K=M/S+5;
struct BLOCK{int bel[M],L[K],R[K],big[K],small[M];BLOCK(){for(int i=1;i<M;i++)bel[i]=(i-1)/S+1;for(int i=1;i<K;i++)L[i]=1+(i-1)*S,R[i]=i*S;}void clear(){memset(big,0,sizeof(big));memset(small,0,sizeof(small));}void update(int x,int v){for(int i=1;i<=bel[x];i++)big[i]+=v;for(int i=L[bel[x]];i<=x;i++)small[i]+=v;}int query(int l,int r){return big[bel[l]+1]+small[l]-big[bel[r+1]+1]-small[r+1];}
}block;
struct node{int l,r,id;bool operator <(const node &A)const{if(l/S!=A.l/S)return l/S<A.l/S;if(l/S&1)return r<A.r;else return r>A.r;}
}Q[M];
int n,m,uni,A[M],tmp[M],bit[M];
ll ans[M],sum1[M],sum2[M],res1[M],res2[M];
vector<node>vi1[M],vi2[M];
void add(int x,int v){while(x<=uni)bit[x]+=v,x+=x&-x;
}
int query(int x){int res=0;while(x)res+=bit[x],x-=x&-x;return res;
}
int main(){#ifndef ONLINE_JUDGEfreopen("jiedai.in","r",stdin);
#endifrd(n),rd(m);for(int i=1;i<=n;i++)rd(A[i]),tmp[i]=A[i];for(int i=1;i<=m;i++)rd(Q[i].l),rd(Q[i].r),Q[i].id=i;sort(tmp+1,tmp+1+n);uni=unique(tmp+1,tmp+1+n)-tmp-1;for(int i=1;i<=n;i++)A[i]=lower_bound(tmp+1,tmp+1+uni,A[i])-tmp;sort(Q+1,Q+1+m);for(int i=1,l=1,r=0;i<=m;i++){if(Q[i].r>r)vi1[l-1].push_back((node){r+1,Q[i].r,i}),r=Q[i].r;if(Q[i].l<l)vi2[r+1].push_back((node){Q[i].l,l-1,i}),l=Q[i].l;if(Q[i].r<r)vi1[l-1].push_back((node){Q[i].r+1,r,i}),r=Q[i].r;if(Q[i].l>l)vi2[r+1].push_back((node){l,Q[i].l-1,i}),l=Q[i].l;}memset(bit,0,sizeof(bit));for(int i=1;i<=n;i++){sum1[i]=sum1[i-1]+i-1-query(A[i]);add(A[i],1);}memset(bit,0,sizeof(bit));for(int i=n;i>=1;i--){sum2[i]=sum2[i+1]+query(A[i]-1);add(A[i],1);}block.clear();for(int i=0;i<n;i++){for(int j=0;j<vi1[i].size();j++){int l=vi1[i][j].l,r=vi1[i][j].r,id=vi1[i][j].id;for(int k=l;k<=r;k++)res1[id]+=block.query(A[k]+1,uni);}block.update(A[i+1],1);}block.clear();for(int i=n+1;i>1;i--){for(int j=0;j<vi2[i].size();j++){int l=vi2[i][j].l,r=vi2[i][j].r,id=vi2[i][j].id;for(int k=l;k<=r;k++)res2[id]+=block.query(1,A[k]-1);}block.update(A[i-1],1);}ll res=0;for(int i=1,l=1,r=0;i<=m;i++){if(Q[i].r>r)res+=sum1[Q[i].r]-sum1[r]-res1[i],r=Q[i].r;if(Q[i].l<l)res+=sum2[Q[i].l]-sum2[l]-res2[i],l=Q[i].l;if(Q[i].r<r)res-=sum1[r]-sum1[Q[i].r]-res1[i],r=Q[i].r;if(Q[i].l>l)res-=sum2[l]-sum2[Q[i].l]-res2[i],l=Q[i].l;ans[Q[i].id]=res;}for(int i=1;i<=m;i++)printf("%lld\n",ans[i]);return (0-0);
}

线段树分治

P5787

题意：n个节点的图中，在k时间内会有m条边先出现后消失，问每一时刻这个图是否是二分图

int ans=0;
struct disjoint_set{static const int SIZE=1e6+5;int top,stk[SIZE],fa[SIZE],sz[SIZE];int len[SIZE],dis[SIZE];void init(int n){top=0;for(int i=1;i<=n;i++)fa[i]=i,sz[i]=1;}int getfa(int x){if(fa[x]==x){dis[x]=0;return x;}int f=getfa(fa[x]);dis[x]=dis[fa[x]]+len[x];return f;}void link(int x,int y){int fx=getfa(x),fy=getfa(y);if(fx!=fy){if(sz[fx]>sz[fy])swap(x,y),swap(fx,fy);stk[++top]=fx;fa[fx]=fy;sz[fy]+=sz[fx];len[fx]=(dis[x]+dis[y]+1)%2;}else if(dis[x]%2==dis[y]%2)ans=0;}void pop_to(int tar){while(top>tar){sz[fa[stk[top]]]-=sz[stk[top]];fa[stk[top]]=stk[top];top--;}}
}U;const int M=1e5+5;
int n,m,k;
vector<pair<int,int> >vi[M<<2];
void update(int a,int b,int x,int y,int l=1,int r=k,int p=1){if(l>b||r<a)return;if(l>=a&&r<=b){vi[p].push_back(make_pair(x,y));return;}int mid=l+r>>1;update(a,b,x,y,l,mid,p<<1);update(a,b,x,y,mid+1,r,p<<1|1);
}
void visit(int l=1,int r=k,int p=1){int lst_top=U.top,lst_ans=ans;for(int i=0;i<vi[p].size();i++)U.link(vi[p][i].first,vi[p][i].second);if(l==r)puts(ans?"Yes":"No");else{int mid=l+r>>1;visit(l,mid,p<<1);visit(mid+1,r,p<<1|1);}U.pop_to(lst_top);ans=lst_ans;
}
int main(){#ifndef ONLINE_JUDGEfreopen("jiedai.in","r",stdin);
#endifrd(n),rd(m),rd(k);for(int i=1;i<=m;i++){int x,y,l,r;rd(x),rd(y),rd(l),rd(r);if(l<r)update(l+1,r,x,y);}ans=1;U.init(n);visit();return (0-0);
}

CDQ分治

P3810

题意： nnn 个三元组 (a,b,c)(a,b,c)(a,b,c) ，f(i)f(i)f(i) 表示三维都小于等于第 iii 个三元组的数量，对于每个 ddd ，求 f(i)=df(i)=df(i)=d 的数量

对于完全相同的三元组，需要进行去重，记录相同三元组的数量

三维偏序

const int M=1e5+5;
int n,k,X[M],Y[M],Z[M],rk[M],cnt[M],res[M],ans[M];
bool cmp(int a,int b){if(X[a]!=X[b])return X[a]<X[b];if(Y[a]!=Y[b])return Y[a]<Y[b];return Z[a]<Z[b];
}
struct node{int y,z,id;bool operator <(const node &A)const{if(y!=A.y)return y<A.y;return id<A.id;}
}Q[M];
int bit[M<<1];
void add(int x,int v){while(x<=k)bit[x]+=v,x+=x&-x;
}
int query(int x){int res=0;while(x)res+=bit[x],x-=x&-x;return res;
}
void CDQ(int l,int r){if(l==r)return;int mid=l+r>>1;CDQ(l,mid),CDQ(mid+1,r);sort(Q+l,Q+r+1);//可以用归并排序代替for(int i=l;i<=r;i++){if(Q[i].id<=mid)add(Q[i].z,cnt[Q[i].id]);else res[Q[i].id]+=query(Q[i].z);}for(int i=l;i<=r;i++)if(Q[i].id<=mid)add(Q[i].z,-cnt[Q[i].id]);
}
int main(){#ifndef ONLINE_JUDGEfreopen("jiedai.in","r",stdin);
#endifrd(n),rd(k);for(int i=1;i<=n;i++)rd(X[i]),rd(Y[i]),rd(Z[i]),rk[i]=i;sort(rk+1,rk+1+n,cmp);int uni=0;for(int i=1;i<=n;i++){if(X[rk[i]]!=X[rk[i-1]]||Y[rk[i]]!=Y[rk[i-1]]||Z[rk[i]]!=Z[rk[i-1]])uni++,Q[uni]=(node){Y[rk[i]],Z[rk[i]],uni};cnt[uni]++;}CDQ(1,uni);for(int i=1;i<=uni;i++)ans[res[i]+cnt[i]-1]+=cnt[i];for(int i=0;i<n;i++)printf("%d\n",ans[i]);return (0-0);
}

四维偏序

HDOJ 5126

1. 向集合中加入三元组(a,b,c)

2. 查询集合中有多少个三元组(a,b,c)，满足a1<=a<=a2, b1<=b<=b2, c1<=c<=c2

对查询操作进行差分，分成8个只查询前缀的操作。

1. 首先对第一维（即时间维）进行CDQ分治。对于分治区间 [L,R][L,R][L,R] ，用左半区间的插入操作更新右半区间的查询操作。

2. 在这个区间内，再对第二维（即 aaa 维）进行排序，排序的第二关键字为第一维（即时间维）。

   再对这个区间内的第二维（即 aaa 维）进行CDQ分治，用排序后 aaa 维前一半小的操作去跟新 aaa 维后一半的查询。

3. 在第二层的CDQ中，再对第三维（即 bbb 维）进行排序，排序的第二关键字为第一维（即时间维），从小到大扫描每一个操作：

4. 对于修改操作，只有当这个操作满足第一维在第一层CDQ的左半区间且第二维在第二层CDQ的左半区间，
   才将当前的第四维（即 ccc 维）加入树状数组中。（第四维在预处理时需要离散）

5. 对于查询操作，只有当这个操作满足第一维在第一层CDQ的右半区间且第二维在第二层CDQ的右半区间，
   才将当前的第四维在树状数组中查询，更新答案。

const int M=5e4+5;
const int LEFT=1;
const int RIGHT=2;
int cas,n,q,uni,OP[M],ans[M],A[M<<1],bit[M<<1];
void add(int x,int v){while(x<=uni)bit[x]+=v,x+=x&-x;
}
int query(int x){int res=0;while(x)res+=bit[x],x-=x&-x;return res;
}
struct node{int id,x,y,z,op,part1,part2;
}Q[M<<3],tmp[M<<3];
bool cmp1(node &A,node &B){if(A.x!=B.x)return A.x<B.x;return A.id<B.id;
}
bool cmp2(node &A,node &B){if(A.y!=B.y)return A.y<B.y;return A.id<B.id;
}
void cdq(int l,int r){if(l==r)return;int mid=l+r>>1;cdq(l,mid),cdq(mid+1,r);for(int i=l;i<=mid;i++)Q[i].part2=LEFT;for(int i=mid+1;i<=r;i++)Q[i].part2=RIGHT;int a=l,b=mid+1,now=l;//使用归并排序提高效率while(a<=mid&&b<=r){if(cmp2(Q[a],Q[b]))tmp[now++]=Q[a++];else tmp[now++]=Q[b++];}while(a<=mid)tmp[now++]=Q[a++];while(b<=r)tmp[now++]=Q[b++];for(int i=l;i<=r;i++)Q[i]=tmp[i];//sort(Q+l,Q+r+1,cmp2);//用归并排序代替for(int i=l;i<=r;i++){if(Q[i].op==0&&Q[i].part1==LEFT&&Q[i].part2==LEFT)add(Q[i].z,1);if(Q[i].op!=0&&Q[i].part1==RIGHT&&Q[i].part2==RIGHT)ans[Q[i].id]+=Q[i].op*query(Q[i].z);}for(int i=l;i<=r;i++)if(Q[i].op==0&&Q[i].part1==LEFT&&Q[i].part2==LEFT)add(Q[i].z,-1);
}
void CDQ(int l,int r){if(l==r)return;int mid=l+r>>1;CDQ(l,mid),CDQ(mid+1,r);for(int i=l;i<=mid;i++)Q[i].part1=LEFT;for(int i=mid+1;i<=r;i++)Q[i].part1=RIGHT;sort(Q+l,Q+r+1,cmp1);cdq(l,r);
}
int main(){#ifndef ONLINE_JUDGEfreopen("jiedai.in","r",stdin);
#endifrd(cas);while(cas--){n=uni=0;memset(ans,0,sizeof(ans));rd(q);for(int i=1;i<=q;i++){int x,y,z,x2,y2,z2;rd(OP[i]),rd(x),rd(y),rd(z);if(OP[i]==1){Q[++n]=(node){i,x,y,z,0,0,0};A[++uni]=z;}else{rd(x2),rd(y2),rd(z2);A[++uni]=z2;A[++uni]=z-1;Q[++n]=(node){i,x2,y2,z2,1,0,0};Q[++n]=(node){i,x-1,y2,z2,-1,0,0};Q[++n]=(node){i,x2,y-1,z2,-1,0,0};Q[++n]=(node){i,x2,y2,z-1,-1,0,0};Q[++n]=(node){i,x-1,y-1,z2,1,0,0};Q[++n]=(node){i,x-1,y2,z-1,1,0,0};Q[++n]=(node){i,x2,y-1,z-1,1,0,0};Q[++n]=(node){i,x-1,y-1,z-1,-1,0,0};}}sort(A+1,A+1+uni);uni=unique(A+1,A+1+uni)-A-1;for(int i=1;i<=n;i++)Q[i].z=lower_bound(A+1,A+1+uni,Q[i].z)-A;CDQ(1,n);for(int i=1;i<=q;i++)if(OP[i]==2)printf("%d\n",ans[i]);}return (0-0);
}

整体二分*

Splay Tree

代码一：权值建树 P3369

1. 插入 v
2. 删除 v (若有多个相同的数，只删除一个)
3. 查询 v 的排名 (比 v 小的数的个数 +1 )
4. 查询排名为 k 的数
5. 求小于 v 的最大的数
6. 求大于 v 的最小的数

struct Splay_Tree{static const int SIZE=1e6+5;int trid,root,par[SIZE],ch[SIZE][2],val[SIZE],cnt[SIZE],sz[SIZE];void up(int x){sz[x]=sz[ch[x][0]]+sz[ch[x][1]]+cnt[x];}void UP(int x){while(x)up(x),x=par[x];}void clean_node(int x){par[x]=ch[x][0]=ch[x][1]=val[x]=cnt[x]=sz[x]=0;}void new_node(int &x,int f,int v){x=++trid;par[x]=f;val[x]=v;cnt[x]=sz[x]=1;ch[x][0]=ch[x][1]=0;UP(x);//单点插入时加上，建树时不加}void rotate(int x){int y=par[x],s=ch[y][0]==x;if(ch[x][s])par[ch[x][s]]=y;ch[y][!s]=ch[x][s];par[x]=par[y];if(par[y])ch[par[y]][ch[par[y]][1]==y]=x;par[y]=x;ch[x][s]=y;up(y),up(x);//order is important}void splay(int x,int goal=0){while(par[x]!=goal){int y=par[x];if(par[y]!=goal){if(ch[par[y]][0]==y^ch[y][0]==x)rotate(x);else rotate(y);}rotate(x);}if(goal==0)root=x;}void insert(int v){if(root==0){new_node(root,0,v);return;}for(int x=root;;x=ch[x][v>val[x]]){if(v==val[x]){cnt[x]++;UP(x);splay(x);return;}if(ch[x][v>val[x]]==0){new_node(ch[x][v>val[x]],x,v);splay(ch[x][v>val[x]]);return;}}}int get_pre(){int x=ch[root][0];if(x==0)return 0;while(ch[x][1])x=ch[x][1];return x;}int Find(int v){for(int x=root;x;x=ch[x][v>val[x]])if(v==val[x])return x;return 0;}void erase(int v){int x=Find(v);if(x==0)return;splay(x);cnt[x]--,sz[x]--;if(cnt[x]>0)return;int y=get_pre();if(y==0)par[root=ch[x][1]]=0;else{splay(y,x);par[ch[x][1]]=y;ch[y][1]=ch[x][1];up(y);par[root=y]=0;}clean_node(x);}int query_kth(int k){int x=root;while(1){if(k<=sz[ch[x][0]])x=ch[x][0];else if(k<=sz[ch[x][0]]+cnt[x]){splay(x);return val[x];}//将查询的点旋转到根，保证复杂度else k-=sz[ch[x][0]]+cnt[x],x=ch[x][1];}}int query_rank(int v){int res=0,s=root;for(int x=root;x;s=x,x=ch[x][v>val[x]]){if(v==val[x]){res+=sz[ch[x][0]];break;}if(v>val[x])res+=sz[ch[x][0]]+cnt[x];}splay(s);//将查询的点旋转到根，保证复杂度return res+1;}int query_less(int v){int res=-2e9,s=root;for(int x=root;x;s=x,x=ch[x][v>val[x]])if(val[x]<v)MAX(res,val[x]);splay(s);//将查询的点旋转到根，保证复杂度return res;}int query_greater(int v){int res=2e9,s=root;for(int x=root;x;s=x,x=ch[x][v>=val[x]])if(val[x]>v)MIN(res,val[x]);splay(s);//将查询的点旋转到根，保证复杂度return res;}
}splay;

代码二：下标建树 P3391

struct Splay_Tree{static const int SIZE=1e6+5;int trid,root,par[SIZE],ch[SIZE][2],flip[SIZE],val[SIZE],sz[SIZE];void up(int x){sz[x]=sz[ch[x][0]]+sz[ch[x][1]]+1;}void UP(int x){while(x)up(x),x=par[x];}void down(int x){if(flip[x]){swap(ch[x][0],ch[x][1]);if(ch[x][0])flip[ch[x][0]]^=1;if(ch[x][1])flip[ch[x][1]]^=1;flip[x]=0;}}void DOWN(int x){static int top,stk[SIZE];while(x)stk[++top]=x,x=par[x];while(top)down(stk[top--]);}void clean_node(int x){par[x]=ch[x][0]=ch[x][1]=flip[x]=val[x]=sz[x]=0;}void new_node(int &x,int f,int v){x=++trid;par[x]=f;val[x]=v;sz[x]=1;ch[x][0]=ch[x][1]=flip[x]=0;
//      UP(x);//单点插入时加上，建树时不加}void build(int &x,int f,int l,int r){int mid=l+r>>1;new_node(x,f,mid);if(l<mid)build(ch[x][0],x,l,mid-1);if(mid<r)build(ch[x][1],x,mid+1,r);up(x);}void rotate(int x){int y=par[x],s=ch[y][0]==x;if(ch[x][s])par[ch[x][s]]=y;ch[y][!s]=ch[x][s];par[x]=par[y];if(par[y])ch[par[y]][ch[par[y]][1]==y]=x;par[y]=x;ch[x][s]=y;up(y),up(x);//order is important}void splay(int x,int goal=0){while(par[x]!=goal){int y=par[x];if(par[y]!=goal){if(ch[par[y]][0]==y^ch[y][0]==x)rotate(x);else rotate(y);}rotate(x);}if(goal==0)root=x;}int Find_kth(int k){int x=root;while(1){down(x);if(k<=sz[ch[x][0]])x=ch[x][0];else if(k<=sz[ch[x][0]]+1)return x;else k-=sz[ch[x][0]]+1,x=ch[x][1];}}void reverse(int l,int r){if(l>=r)return;int x=Find_kth(l-1),y=Find_kth(r+1);splay(x),splay(y,x);flip[ch[y][0]]^=1;}int query_kth(int k){int x=root;while(1){down(x);if(k<=sz[ch[x][0]])x=ch[x][0];else if(k<=sz[ch[x][0]]+1){splay(x);return val[x];}//将查询的点旋转到根，保证复杂度else k-=sz[ch[x][0]]+1,x=ch[x][1];}}
}splay;
int main(){#ifndef ONLINE_JUDGEfreopen("jiedai.in","r",stdin);
#endifint n,q,l,r;rd(n),rd(q);splay.build(splay.root,0,0,n+1);//将0和n+1插入方便提取区间while(q--){rd(l),rd(r);splay.reverse(l+1,r+1);}for(int i=1;i<=n;i++)printf("%d%c",splay.query_kth(i+1),i==n?'\n':' ');return (0-0);
}

Link Cut Tree

代码一：动态树 P3690

struct Link_Cut_Tree{static const int SIZE=1e5+5;int par[SIZE],ch[SIZE][2],flip[SIZE],val[SIZE],sum[SIZE];bool is_root(int x){//splay's rootreturn ch[par[x]][0]!=x&&ch[par[x]][1]!=x;}void up(int x){sum[x]=val[x]^sum[ch[x][0]]^sum[ch[x][1]];}void down(int x){if(flip[x]){swap(ch[x][0],ch[x][1]);if(ch[x][0])flip[ch[x][0]]^=1;if(ch[x][1])flip[ch[x][1]]^=1;flip[x]=0;}}void DOWN(int x){static int top,stk[SIZE];for(stk[++top]=x;!is_root(x);x=par[x])stk[++top]=par[x];while(top)down(stk[top--]);}void rotate(int x){int y=par[x],s=ch[y][0]==x;if(ch[x][s])par[ch[x][s]]=y;ch[y][!s]=ch[x][s];par[x]=par[y];if(!is_root(y))ch[par[y]][ch[par[y]][1]==y]=x;par[y]=x;ch[x][s]=y;up(y),up(x);//order is important}void splay(int x){DOWN(x);while(!is_root(x)){int y=par[x];if(!is_root(y)){if(ch[par[y]][0]==y^ch[y][0]==x)rotate(x);else rotate(y);}rotate(x);}up(x);}void Access(int x){int s=0;while(x){splay(x);ch[x][1]=s;up(x);x=par[s=x];}}void make_root(int x){Access(x);splay(x);flip[x]^=1;}void split(int x,int y){make_root(x);Access(y);splay(y);}int find_root(int x){//tree's rootAccess(x);splay(x);while(down(x),ch[x][0])x=ch[x][0];splay(x);return x;}void link(int x,int y){make_root(x);if(find_root(y)==x)return;par[x]=y;splay(x);}void cut(int x,int y){if(find_root(x)!=find_root(y))return;split(x,y);if(par[x]!=y||ch[x][1])return;par[x]=ch[y][0]=0;up(y);}void update(int x,int v){splay(x);val[x]=v;up(x);}int query(int x,int y){split(x,y);return sum[y];}
}lct;

代码二：动态图连通性*

KD Tree

代码一：单点插入+矩阵求和 P4148

操作1：在(x,y)上加v

操作2：查询(x1,y1)到(x2,y2)的矩阵内的数字和

强制在线

1. K-D树：

   以2_D树为例，树上的每个节点都“管理”着一个矩形区域，每个节点上存一个点的坐标，这个坐标一定位于这个节点“管理”的区域内，

   接下来按照这个点的某一维（x或y）上的值将矩阵划分为两个子矩阵（对应由两个子节点“管理”）。

   这样就将多维空间不断地通过某一维度的中位数点进行当前空间的分割，于是平衡时树的高度是 O(log2n)O(log_2{n})O(log2​n) 的。

   KD树的查询其实是暴力+剪枝，复杂度可以当作 O(n))O(\sqrt{n}))O(n​))。

2. 需要维护的基础信息：

   1. treetreetree：存储当前节点的坐标以及权值
   2. lson,rsonlson,rsonlson,rson：当前节点的左右儿子
   3. sumsumsum：当前子树中所有节点的权值和
   4. sizesizesize：当前子树的节点数量
   5. mi[i]mi[i]mi[i]：当前子树中第iii维的最小值
   6. mx[i]mx[i]mx[i]：当前子树中第iii维的最大值

3. 常用基础函数：

   1. nth_element(tmp+l,tmp+mid,tmp+r+1);用O(length)的时间选出中点并将数组分成两边
   2. void rebuild(int &x,int l,int r);对子树进行重构
   3. void to_array(int x);将子树中的所有节点储存入临时数组中
   4. void check(int &x);检查当前子树是否平衡，平衡系数0.75
   5. bool outside(int a1,int a2,int b1,int b2,int x1,int x2,int y1,int y2);判断当前区域(x1,y1)-(x2,y2)是否在查询区域(a1,b1)-(a2,b2)之外
   6. bool inside(int a1,int a2,int b1,int b2,int x1,int x2,int y1,int y2);判断当前区域(x1,y1)-(x2,y2)是否在查询区域(a1,b1)-(a2,b2)之内

const int SIZE=2e5+5;
int TYPE,jiedai;
struct kd_node{int dim[2],val;bool operator <(const kd_node &A)const{return dim[TYPE]<A.dim[TYPE];}
}tree[SIZE],tmp[SIZE];
int root,trid,dfsid,rub,bin[SIZE];
int type[SIZE],ls[SIZE],rs[SIZE],mx[SIZE][2],mi[SIZE][2],sz[SIZE];
ll sum[SIZE];
void up(int x){for(int i=0;i<2;i++){mx[x][i]=mi[x][i]=tree[x].dim[i];if(ls[x])MAX(mx[x][i],mx[ls[x]][i]),MIN(mi[x][i],mi[ls[x]][i]);if(rs[x])MAX(mx[x][i],mx[rs[x]][i]),MIN(mi[x][i],mi[rs[x]][i]);}sum[x]=sum[ls[x]]+sum[rs[x]]+tree[x].val;sz[x]=sz[ls[x]]+sz[rs[x]]+1;
}
void new_node(int &x,int ty,const kd_node &rhs){if(rub)x=bin[rub--];else x=++trid;type[x]=ty;tree[x]=rhs;ls[x]=rs[x]=0;up(x);
}
void clean_node(int x){bin[++rub]=x;ls[x]=rs[x]=sz[x]=sum[x]=0;
}
void rebuild(int &x,int l,int r){int mid=l+r>>1;TYPE=(++jiedai)&1;nth_element(tmp+l,tmp+mid,tmp+r+1);new_node(x,TYPE,tmp[mid]);if(l<mid)rebuild(ls[x],l,mid-1);if(mid<r)rebuild(rs[x],mid+1,r);up(x);
}
void to_array(int x){if(ls[x])to_array(ls[x]);if(rs[x])to_array(rs[x]);tmp[++dfsid]=tree[x];clean_node(x);
}
void insert(int &x,const kd_node &rhs){if(x==0){new_node(x,(++jiedai)&1,rhs);return;}if(rhs.dim[type[x]]<=tree[x].dim[type[x]])insert(ls[x],rhs);else insert(rs[x],rhs);up(x);if(sz[ls[x]]>sz[x]*0.75||sz[rs[x]]>sz[x]*0.75){dfsid=0;to_array(x);rebuild(x,1,dfsid);}
}
ll query(int x,int a1,int a2,int b1,int b2){//a1<=x<=a2 b1<=y<=b2if(x==0)return 0;if(mi[x][0]>a2||mx[x][0]<a1||mi[x][1]>b2||mx[x][1]<b1)return 0;if(mi[x][0]>=a1&&mx[x][0]<=a2&&mi[x][1]>=b1&&mx[x][1]<=b2)return sum[x];if(tree[x].dim[0]>=a1&&tree[x].dim[0]<=a2&&tree[x].dim[1]>=b1&&tree[x].dim[1]<=b2)return query(ls[x],a1,a2,b1,b2)+query(rs[x],a1,a2,b1,b2)+tree[x].val;return query(ls[x],a1,a2,b1,b2)+query(rs[x],a1,a2,b1,b2);
}
int main(){#ifndef ONLINE_JUDGEfreopen("jiedai.in","r",stdin);
#endifint n,op,x,y,v,x1,x2,y1,y2,lstans=0;rd(n);while(1){rd(op);if(op==1){rd(x),rd(y),rd(v);x^=lstans,y^=lstans,v^=lstans;insert(root,(kd_node){x,y,v});}else if(op==2){rd(x1),rd(y1),rd(x2),rd(y2);x1^=lstans,y1^=lstans,x2^=lstans,y2^=lstans;printf("%lld\n",lstans=query(root,x1,x2,y1,y2));}else if(op==3)break;}return (0-0);
}

代码二：k远点对 P4357

题意：
给定二维平面上的NNN个点，求第KKK远的无序点对。

思路：

1. 别问我为什么想到用K-D Tree的，因为是看了题解的。
2. 本题没有插入、删除等高级操作，仅仅建树和查询，代码简洁。
3. 进入正题：考虑暴力，暴力遍历对于每个点而言能形成的所有点对，显然复杂度为O(n2)O(n^2)O(n2)，不可行，接下来考虑剪枝。
4. 首先，K=min(100,n∗(n+1)2)K=min(100,\frac{n*(n+1)}{2})K=min(100,2n∗(n+1)​)，所以先在小顶堆中插入2∗K2*K2∗K个000，在后续暴力搜索前2∗K2*K2∗K大点对的过程中逐渐把它们poppoppop掉。
5. 对这NNN个点的每个点而言，都从K-D Tree的根节点往下遍历，每到一个节点，先计算当前节点与这个点的距离，并更新小顶堆，然后进入到剪枝的关键步骤。
6. 我们考虑每个节点的左右儿子，分别利用左右儿子的每个维度最大最小边界来计算可能的最远点，若当前子空间最远的点都无法对小顶堆进行更新，则不需要进入这个空间了！
7. 另一点剪枝：先遍历左右子空间中最远可能点更远的子空间，这样也许就不用再遍历另外一个空间啦。
8. 复杂度的话。。。还不会算，据说是O(n32)O(n^{\frac{3}{2}})O(n23​)。

const int SIZE=1e5+5;
int TYPE,jiedai;
struct kd_node{int dim[2];bool operator <(const kd_node &A)const{return dim[TYPE]<A.dim[TYPE];}ll calc(kd_node &A){return 1ll*(dim[0]-A.dim[0])*(dim[0]-A.dim[0])+1ll*(dim[1]-A.dim[1])*(dim[1]-A.dim[1]);}ll calc(int mx[2],int mi[2]){int d1=max(abs(dim[0]-mx[0]),abs(dim[0]-mi[0]));int d2=max(abs(dim[1]-mx[1]),abs(dim[1]-mi[1]));return 1ll*d1*d1+1ll*d2*d2;}
}tree[SIZE],tmp[SIZE];
int root,trid;
int type[SIZE],ls[SIZE],rs[SIZE],mx[SIZE][2],mi[SIZE][2];
void up(int x){for(int i=0;i<2;i++){mx[x][i]=mi[x][i]=tree[x].dim[i];if(ls[x])MAX(mx[x][i],mx[ls[x]][i]),MIN(mi[x][i],mi[ls[x]][i]);if(rs[x])MAX(mx[x][i],mx[rs[x]][i]),MIN(mi[x][i],mi[rs[x]][i]);}
}
void new_node(int &x,int ty,const kd_node &rhs){x=++trid;type[x]=ty;tree[x]=rhs;ls[x]=rs[x]=0;up(x);
}
void build(int &x,int l,int r){int mid=l+r>>1;TYPE=(++jiedai)&1;nth_element(tmp+l,tmp+mid,tmp+r+1);new_node(x,TYPE,tmp[mid]);if(l<mid)build(ls[x],l,mid-1);if(mid<r)build(rs[x],mid+1,r);up(x);
}
priority_queue<ll,vector<ll>,greater<ll>>Q;
void query(int x,kd_node &rhs){ll t=rhs.calc(tree[x]),t1=0,t2=0;if(t>Q.top())Q.pop(),Q.push(t);if(ls[x])t1=rhs.calc(mx[ls[x]],mi[ls[x]]);if(rs[x])t2=rhs.calc(mx[rs[x]],mi[rs[x]]);if(t1>t2){if(t1>Q.top())query(ls[x],rhs);if(t2>Q.top())query(rs[x],rhs);}else{if(t2>Q.top())query(rs[x],rhs);if(t1>Q.top())query(ls[x],rhs);}
}
int n,k;
int main(){#ifndef ONLINE_JUDGEfreopen("jiedai.in","r",stdin);
#endifrd(n),rd(k);for(int i=1;i<=n;i++)rd(tmp[i].dim[0]),rd(tmp[i].dim[1]);build(root,1,n);for(int i=1;i<=k+k;i++)Q.push(0);for(int i=1;i<=n;i++)query(root,tree[i]);printf("%lld\n",Q.top());return (0-0);
}

左偏树*

双栈模拟队列

loj6515

const int M=5e4+5;
struct Stack{int st,top;pair<int,int>stk[M];void insert(pair<int,int>x){stk[++top]=x;}void pop(){top--;}bool empty(){return st==top;}
}stk1,stk2;
int casid,m,p;
ll dp1[M][505],dp2[M][505],Q[M],A[M],B[M];
void DP(Stack &stk,ll dp[M][505],int i){memset(dp[i],-1,sizeof(dp[i]));for(int j=0;j<p;j++)if(~dp[i-1][j]){MAX(dp[i][j],dp[i-1][j]);MAX(dp[i][(j+stk.stk[i].first)%p],dp[i-1][j]+stk.stk[i].second);}
}
void all_DP(Stack &stk,ll dp[M][505]){memset(dp[stk.st],-1,sizeof(dp[stk.st]));dp[stk.st][0]=0;for(int i=stk.st+1;i<=stk.top;i++)DP(stk,dp,i);
}
void rebuild(Stack &stk1,Stack &stk2){int mid=(stk1.st+stk1.top+1)/2;for(int i=mid;i>stk1.st;i--)stk2.insert(stk1.stk[i]);stk1.st=mid;
}
int main(){#ifndef ONLINE_JUDGEfreopen("jiedai.in","r",stdin);
#endifmemset(dp1,-1,sizeof(dp1));memset(dp2,-1,sizeof(dp2));rd(casid),rd(m),rd(p);dp1[0][0]=dp2[0][0]=0;for(int k=1;k<=m;k++){char str[15];int a,b;scanf("%s",str);if(0);else if(str[0]=='I'&&str[1]=='F'){rd(a),rd(b);stk1.insert(make_pair(a,b));DP(stk1,dp1,stk1.top);}else if(str[0]=='I'&&str[1]=='G'){rd(a),rd(b);stk2.insert(make_pair(a,b));DP(stk2,dp2,stk2.top);}else if(str[0]=='D'&&str[1]=='F'){if(stk1.empty())rebuild(stk2,stk1),all_DP(stk1,dp1),all_DP(stk2,dp2);stk1.pop();}else if(str[0]=='D'&&str[1]=='G'){if(stk2.empty())rebuild(stk1,stk2),all_DP(stk1,dp1),all_DP(stk2,dp2);stk2.pop();}else if(str[0]=='Q'&&str[1]=='U'){rd(a),rd(b);for(int i=0;i<p;i++)A[i]=dp1[stk1.top][i];for(int i=0;i<p;i++)B[i]=dp2[stk2.top][i];int L=1,R=0;ll ans=-1;for(int i=a;i<=b;i++){while(L<=R&&A[i]>A[Q[R]])R--;Q[++R]=i;}for(int i=p-1;i>=0;i--){while(L<=R&&A[(p+b-i)%p]>=A[Q[R]%p])R--;Q[++R]=p+b-i;while(L<=R&&Q[L]<p+a-i)L++;if(L<=R&&~A[Q[L]%p]&&~B[i])MAX(ans,A[Q[L]%p]+B[i]);}printf("%lld\n",ans);}}return (0-0);
}

哈希表实现map

代码一：STL版本

unordered_map<int,int>mp;

代码二：手写版本（int->int）

struct hash_map{static const int SIZE=1e6+5;int tot,head[SIZE],key[SIZE*1],value[SIZE*1],nxt[SIZE*1];int Hash(int x){if(x<0)x=-x;return x%(SIZE-1)+1;}int& operator [](int x){int pos=Hash(x);for(int i=head[pos];i;i=nxt[i])if(key[i]==x)return value[i];key[++tot]=x,value[tot]=0,nxt[tot]=head[pos],head[pos]=tot;return value[tot];}void clear(){tot=0;memset(head,0,sizeof(head));}hash_map(){clear();}
};

代码三：手写类模板版本（整数->所有类型）

template<class KEY,class VALUE>struct hash_map{static const int SIZE=1e6+5;int tot,head[SIZE],nxt[SIZE*1];KEY key[SIZE*1];VALUE value[SIZE*1];int Hash(KEY x){if(x<0)x=-x;return x%(SIZE-1)+1;}VALUE& operator [](KEY x){int pos=Hash(x);for(int i=head[pos];i;i=nxt[i])if(key[i]==x)return value[i];key[++tot]=x,value[tot]=0,nxt[tot]=head[pos],head[pos]=tot;return value[tot];}void clear(){tot=0;memset(head,0,sizeof(head));}hash_map(){clear();}
};

支持任意删除的大顶堆

struct Multiset{priority_queue<int>A,B;void clear(){while(!A.empty())A.pop();while(!B.empty())B.pop();}bool empty(){return A.empty();}int top(){if(A.empty())return -1e9;else return A.top();}void pop(){if(A.empty())return;A.pop();while(!A.empty()&&!B.empty()&&A.top()==B.top())A.pop(),B.pop();}void push(int x){A.push(x);}void erase(int x){B.push(x);while(!A.empty()&&!B.empty()&&A.top()==B.top())A.pop(),B.pop();}
};

树论

树的直径

const int SIZE=1e6+5;
int n,tot,head[SIZE],to[SIZE<<1],co[SIZE<<1],nxt[SIZE<<1];
inline void add_edge(int a,int b,int c){to[++tot]=b;co[tot]=c;nxt[tot]=head[a];head[a]=tot;
}
int mxID;
ll mxdis;
void dfs(int x,int f,ll d){if(d>mxdis)mxdis=d,mxID=x;for(int i=head[x];i;i=nxt[i]){int y=to[i];if(y==f)continue;dfs(y,x,d+co[i]);}
}
int main(){#ifndef ONLINE_JUDGEfreopen("jiedai.in","r",stdin);
#endifrd(n);for(int i=1;i<n;i++){int a,b,c;rd(a),rd(b),rd(c);add_edge(a,b,c),add_edge(b,a,c);}int s=0,t=0;mxdis=-1,dfs(1,0,0),s=mxID;mxdis=-1,dfs(s,0,0),t=mxID;printf("%d %d %lld\n",s,t,mxdis);return (0-0);
}

树的重心

1. 删除重心后所得的所有子树，节点数不超过原树的1/2，一棵树最多有两个重心，且它们相邻
2. 树中所有节点到重心的距离之和最小，如果有两个重心，那么它们距离之和相等
3. 两个树通过一条边合并，新的重心在原树两个重心的路径上
4. 树删除或添加一个叶子节点，重心最多只移动一条边

const int SIZE=1e6+5;
int n,tot,head[SIZE],to[SIZE<<1],nxt[SIZE<<1];
inline void add_edge(int a,int b){to[++tot]=b;nxt[tot]=head[a];head[a]=tot;
}
int sz[SIZE],miID,mival;
void dfs(int x,int f){sz[x]=1;int res=0;for(int i=head[x];i;i=nxt[i]){int y=to[i];if(y==f)continue;dfs(y,x);sz[x]+=sz[y];MAX(res,sz[y]);}MAX(res,n-sz[x]);if(res<mival)mival=res,miID=x;
}
int main(){#ifndef ONLINE_JUDGEfreopen("jiedai.in","r",stdin);
#endifrd(n);for(int i=1;i<n;i++){int a,b;rd(a),rd(b);add_edge(a,b),add_edge(b,a);}mival=1e9,miID=-1;dfs(1,0);printf("%d %d\n",miID,mival);return (0-0);
}

树上倍增+LCA

const int SIZE=5e5+5;
int tot,head[SIZE],to[SIZE<<1],nxt[SIZE<<1];
int fa[SIZE][20],deep[SIZE];
void add_edge(int a,int b){to[++tot]=b;nxt[tot]=head[a];head[a]=tot;
}
void dfs(int x,int f){deep[x]=deep[f]+1;fa[x][0]=f;for(int i=1;i<20;i++)fa[x][i]=fa[fa[x][i-1]][i-1];for(int i=head[x];i;i=nxt[i]){int y=to[i];if(y==f)continue;dfs(y,x);}
}
int LCA(int x,int y){if(deep[x]<deep[y])swap(x,y);int step=deep[x]-deep[y];for(int i=0;i<20;i++)if(step&1<<i)x=fa[x][i];if(x==y)return x;for(int i=19;i>=0;i--)if(fa[x][i]!=fa[y][i])x=fa[x][i],y=fa[y][i];return fa[x][0];
}

树链剖分+LCA

const int SIZE=5e5+5;
int tot,head[SIZE],to[SIZE<<1],nxt[SIZE<<1];
int dfsid,fa[SIZE],deep[SIZE],sz[SIZE],son[SIZE],top[SIZE],dfn[SIZE],reid[SIZE];
void add_edge(int a,int b){to[++tot]=b;nxt[tot]=head[a];head[a]=tot;
}
void dfs(int x,int f){deep[x]=deep[f]+1;fa[x]=f;sz[x]=1;son[x]=0;for(int i=head[x];i;i=nxt[i]){int y=to[i];if(y==f)continue;dfs(y,x);sz[x]+=sz[y];if(sz[y]>sz[son[x]])son[x]=y;}
}
void chain(int x,int k){top[x]=k;reid[dfn[x]=++dfsid]=x;if(son[x])chain(son[x],k);for(int i=head[x];i;i=nxt[i]){int y=to[i];if(y==fa[x]||y==son[x])continue;chain(y,y);}
}
int LCA(int a,int b){while(top[a]!=top[b]){if(deep[top[a]]<deep[top[b]])swap(a,b);a=fa[top[a]];}return deep[a]<deep[b]?a:b;
}

RMQ+LCA

O(nlogn) 预处理 O(1) 求LCA

const int SIZE=5e5+5;
int n,tot,head[SIZE],to[SIZE<<1],nxt[SIZE<<1];
int dfsid,deep[SIZE],dfn[SIZE],st[SIZE<<1][20];
void add_edge(int a,int b){to[++tot]=b;nxt[tot]=head[a];head[a]=tot;
}
void dfs(int x,int f){deep[x]=deep[f]+1;dfn[x]=++dfsid;st[dfsid][0]=x;for(int i=head[x];i;i=nxt[i]){int y=to[i];if(y==f)continue;dfs(y,x);st[++dfsid][0]=x;}
}
void RMQ(){for(int i=dfsid;i>=1;i--)for(int j=1;i+(1<<j)-1<=dfsid;j++)st[i][j]=deep[st[i][j-1]]<deep[st[i+(1<<j-1)][j-1]]?st[i][j-1]:st[i+(1<<j-1)][j-1];
}
int LCA(int x,int y){int l=dfn[x],r=dfn[y];if(l>r)swap(l,r);int k=log2(r-l+1);return deep[st[l][k]]<deep[st[r-(1<<k)+1][k]]?st[l][k]:st[r-(1<<k)+1][k];
}

点分治

P3806

n个节点的边权树，m次询问，每次询问树上距离为k的点对是否存在

1≤n≤104,1≤m≤100,1≤k≤107,边权1≤w≤1041≤n≤10^4,1≤m≤100,1≤k≤10^7,边权1≤w≤10^41≤n≤104,1≤m≤100,1≤k≤107,边权1≤w≤104

const int M=10005;
const int K=1e7+5;
int tot,head[M],to[M<<1],co[M<<1],nxt[M<<1];
inline void add_edge(int a,int b,int c){to[++tot]=b;co[tot]=c;nxt[tot]=head[a];head[a]=tot;
}
int n,m,Q[M],ans[M],cnt[K];
int all,used[M],tmp[M],sz[M],mx[M];
void center(int x,int f){sz[x]=1;mx[x]=0;tmp[++all]=x;for(int i=head[x];i;i=nxt[i]){int y=to[i];if(y==f||used[x])continue;center(y,x);sz[x]+=sz[y];MAX(mx[x],sz[y]);}
}
void update(int x,int f,int d,int v){if(d<K)cnt[d]+=v;for(int i=head[x];i;i=nxt[i]){int y=to[i];if(y==f||used[y])continue;update(y,x,d+co[i],v);}
}
void query(int x,int f,int d){for(int i=1;i<=m;i++)if(Q[i]-d>=0&&cnt[Q[i]-d])ans[i]=1;for(int i=head[x];i;i=nxt[i]){int y=to[i];if(y==f||used[y])continue;query(y,x,d+co[i]);}
}
void divide(int x){all=0;center(x,0);for(int i=1;i<=all;i++){MAX(mx[tmp[i]],all-sz[tmp[i]]);if(mx[tmp[i]]<mx[x])x=tmp[i];}cnt[0]=1;for(int i=head[x];i;i=nxt[i]){int y=to[i];if(used[y])continue;query(y,x,co[i]);update(y,x,co[i],1);}update(x,0,0,-1);used[x]=1;for(int i=head[x];i;i=nxt[i]){int y=to[i];if(used[y])continue;divide(y);}
}
int main(){#ifndef ONLINE_JUDGEfreopen("jiedai.in","r",stdin);
#endifrd(n),rd(m);for(int i=1;i<n;i++){int a,b,c;rd(a),rd(b),rd(c);add_edge(a,b,c),add_edge(b,a,c);}for(int i=1;i<=m;i++)rd(Q[i]);divide(1);for(int i=1;i<=m;i++)puts(ans[i]?"AYE":"NAY");return (0-0);
}

DSU on tree

const int M=1e5+5;
int CAS,n,k,q,A[M],B[M],tot,head[M],to[M<<1],nxt[M<<1];
inline void add_edge(int a,int b){to[++tot]=b;nxt[tot]=head[a];head[a]=tot;
}
int sz[M],fa[M],son[M];
void dfs(int x,int f){fa[x]=f;sz[x]=1;son[x]=0;for(int i=head[x];i;i=nxt[i]){int y=to[i];if(y==f)continue;dfs(y,x);sz[x]+=sz[y];if(sz[y]>sz[son[x]])son[x]=y;}
}
int now,ans[M],cnt[M];
void insert(int val,int sign){if(cnt[val]==k)now--;cnt[val]+=sign;if(cnt[val]==k)now++;
}
void rdfs(int x,int sign){insert(A[x],sign);for(int i=head[x];i;i=nxt[i]){int y=to[i];if(y==fa[x])continue;rdfs(y,sign);}
}
void dsu(int x){for(int i=head[x];i;i=nxt[i]){int y=to[i];if(y==fa[x]||y==son[x])continue;dsu(y);rdfs(y,-1);}if(son[x])dsu(son[x]);insert(A[x],1);for(int i=head[x];i;i=nxt[i]){int y=to[i];if(y==fa[x]||y==son[x])continue;rdfs(y,1);}ans[x]=now;
}int main(){#ifndef ONLINE_JUDGEfreopen("jiedai.in","r",stdin);
#endifrd(CAS);for(int cas=1;cas<=CAS;cas++){memset(head,0,sizeof(head));memset(cnt,0,sizeof(cnt));now=tot=0;rd(n),rd(k);for(int i=1;i<=n;i++)rd(A[i]),B[i]=A[i];sort(B+1,B+1+n);int uni=unique(B+1,B+1+n)-B-1;for(int i=1;i<=n;i++)A[i]=lower_bound(B+1,B+1+uni,A[i])-B;for(int i=1;i<n;i++){int a,b;rd(a),rd(b);add_edge(a,b),add_edge(b,a);}dfs(1,0);dsu(1);printf("Case #%d:\n",cas);rd(q);while(q--){int x;rd(x);printf("%d\n",ans[x]);}if(cas<CAS)puts("");}return (0-0);
}

虚树*

图论

最短路

Dijkstra

单源最短路径 P4779

const int SIZE=2e5+5;
int n,m,s;
int tot,head[SIZE],to[SIZE<<1],co[SIZE<<1],nxt[SIZE<<1];
ll dis[SIZE];
inline void add_edge(int a,int b,int c){to[++tot]=b;co[tot]=c;nxt[tot]=head[a];head[a]=tot;
}
struct node{int x;ll d;bool operator <(const node &A)const{return d>A.d;}
};
void Dij(int s){memset(dis,-1,sizeof(dis));priority_queue<node>Q;dis[s]=0;Q.push((node){s,0});while(!Q.empty()){int x=Q.top().x;ll d=Q.top().d;Q.pop();if(d>dis[x])continue;for(int i=head[x];i;i=nxt[i]){int y=to[i];if(dis[y]==-1||d+co[i]<dis[y]){dis[y]=d+co[i];Q.push((node){y,d+co[i]});}}}
}

Floyd

输入的边为有向边

 for(int i=1;i<=n;i++)for(int j=1;j<=n;j++)dis[i][j]=1e18;for(int i=1;i<=n;i++)dis[i][i]=0;for(int i=1;i<=m;i++){int a,b,c;rd(a),rd(b),rd(c);MIN(dis[a][b],1ll*c);}for(int k=1;k<=n;k++)//第一个for枚举中转节点for(int i=1;i<=n;i++)for(int j=1;j<=n;j++)MIN(dis[i][j],dis[i][k]+dis[k][j]);

SPFA

const int SIZE=2e5+5;
int n,m,s;
int tot,head[SIZE],to[SIZE<<1],co[SIZE<<1],nxt[SIZE<<1];
inline void add_edge(int a,int b,int c){to[++tot]=b;co[tot]=c;nxt[tot]=head[a];head[a]=tot;
}
int used[SIZE],in[SIZE];
ll dis[SIZE];
bool SPFA(int s){memset(dis,-1,sizeof(dis));memset(used,0,sizeof(used));memset(in,0,sizeof(in));queue<int>Q;dis[s]=0;Q.push(s),used[s]=1,in[s]++;while(!Q.empty()){int x=Q.front();Q.pop();ll d=dis[x];used[x]=0;if(in[x]>n+1)return 0;//存在负环for(int i=head[x];i;i=nxt[i]){int y=to[i];if(dis[y]==-1||d+co[i]<dis[y]){dis[y]=d+co[i];if(!used[y])Q.push(y),used[y]=1,in[y]++;}}}return 1;
}

差分约束

n个未知数，m个不等式，求任意一组满足这个不等式组的解

{xa1−xb1≤c1xa2−xb2≤c2⋮xam−xbm≤cm\begin{cases} x_{a_1}-x_{b_1}≤c_1\\ x_{a_2}-x_{b_2}≤c_2\\ \quad \vdots\\ x_{a_m}-x_{b_m}≤c_m \end{cases}⎩⎪⎪⎪⎪⎨⎪⎪⎪⎪⎧​xa1​​−xb1​​≤c1​xa2​​−xb2​​≤c2​⋮xam​​−xbm​​≤cm​​

差分约束系统中的每个约束条件 xi−xj≤ckx_i-x_j≤c_kxi​−xj​≤ck​ 都可以变形成 xi≤xj+ckx_i≤x_j+c_kxi​≤xj​+ck​ ，这与单源最短路中的三角形不等式 dis[xi]≤dis[xj]+ckdis[x_i]≤dis[x_j]+c_kdis[xi​]≤dis[xj​]+ck​ 非常相似。因此，我们可以把每个变量 xix_ixi​ 看做图中的一个结点，对于每个约束条件 ，从结点 jjj 向结点 iii 连一条长度为 ckc_kck​ 的有向边。

如果 {a1,a2,...,an}\{a_1,a_2,...,a_n\}{a1​,a2​,...,an​} 是该差分约束系统的一组解，那么对于任意的常数 ddd ，{a1+d,a2+d,...,an+d}\{a_1+d,a_2+d,...,a_n+d\}{a1​+d,a2​+d,...,an​+d} 显然也是该差分约束系统的一组解。

设 dis[0]=0dis[0]=0dis[0]=0 并向每一个点连一条权重为 0 的边，用SPFA算法跑单源最短路，若图中存在负环，则给定的差分约束系统无解，否则， xi=dis[i]x_i=dis[i]xi​=dis[i] 为该差分约束系统的一组解。

const int SIZE=5005;
int n,m;
int tot,head[SIZE],to[SIZE<<1],co[SIZE<<1],nxt[SIZE<<1];
inline void add_edge(int a,int b,int c){to[++tot]=b;co[tot]=c;nxt[tot]=head[a];head[a]=tot;
}
int used[SIZE],in[SIZE];
ll dis[SIZE];
bool SPFA(int s){memset(dis,-1,sizeof(dis));memset(used,0,sizeof(used));memset(in,0,sizeof(in));queue<int>Q;dis[s]=0;Q.push(s),used[s]=1,in[s]++;while(!Q.empty()){int x=Q.front();Q.pop();ll d=dis[x];used[x]=0;if(in[x]>n+1)return 0;for(int i=head[x];i;i=nxt[i]){int y=to[i];if(dis[y]==-1||d+co[i]<dis[y]){dis[y]=d+co[i];if(!used[y])Q.push(y),used[y]=1,in[y]++;}}}return 1;
}
int main(){#ifndef ONLINE_JUDGEfreopen("jiedai.in","r",stdin);
#endifrd(n),rd(m);for(int i=1;i<=m;i++){int a,b,c;rd(a),rd(b),rd(c);add_edge(b,a,c);}for(int i=1;i<=n;i++)add_edge(0,i,0);if(SPFA(0)==0)puts("NO");else for(int i=1;i<=n;i++)printf("%lld%c",dis[i],i==n?'\n':' ');return (0-0);
}

二分图

二分图最大匹配

匈牙利算法 P3386

二分图左边n个点，右边m个点，e条边，求二分图最大匹配

const int SIZE=505;
int n,m,e;
int match[SIZE<<1],mark[SIZE<<1];
int tot,head[SIZE],to[SIZE*SIZE],nxt[SIZE*SIZE];
inline void add_edge(int a,int b){to[++tot]=b;nxt[tot]=head[a];head[a]=tot;
}
bool dfs(int x){for(int i=head[x];i;i=nxt[i]){int y=to[i];if(mark[y])continue;mark[y]=1;if(!match[y]||dfs(match[y])){match[x]=y;match[y]=x;return 1;}}return 0;
}
int main(){#ifndef ONLINE_JUDGEfreopen("jiedai.in","r",stdin);
#endifrd(n),rd(m),rd(e);for(int i=1;i<=e;i++){int a,b;rd(a),rd(b);add_edge(a,b+n);}int ans=0;for(int i=1;i<=n;i++){memset(mark,0,sizeof(mark));if(!match[i]&&dfs(i))ans++;}printf("%d\n",ans);return (0-0);
}

二分图最小点覆盖

选出最少的点，满足每条边至少有一个端点被选

二分图中，最小点覆盖 = 最大匹配

二分图最小边覆盖

选出最少的边，满足每个点至少有一条边被选

二分图中，最小边覆盖 = 总点数 - 最大匹配

二分图最大独立集

选出最多的点，满足两两之间没有边相连

二分图中，最大独立集 = 总点数 - 最大匹配

二分图最大团

选出最多的点，满足二分图两侧的点对之间都有边相连

二分图中，最大团 = 补图的最大独立集

DAG最小不相交路径覆盖

有向无环图中，找出最少的路径，使得这些路径经过了所有的点，且每一条路径经过的顶点各不相同

把原图中的点 V 拆成两个点 Vx 和 Vy ，对于原图中的边 A->B ，在新图中建立一条 Ax->By 的边

那么，最小不相交路径覆盖= 原图的点数 - 新图的最大匹配

DAG最小可相交路径覆盖

有向无环图中，找出最少的路径，使得这些路径经过了所有的点，路径可以相交

先用floyd求出原图的传递闭包：如果a到b有路， 那么就添加一条从a到b的边。然后就转化成了最小不相交路径覆盖问题

DAG最大独立集

选出最多的点，满足两两之间没有边相连

有向无环图中，最大独立集 = 最小不相交路径覆盖

DAG最长反链

选出最多的点，满足两两之间没有路径相连

有向无环图中，最长反链 = 最小可相交路径覆盖

网络流

最大流

n个点，m条有向边，源点为s，汇点为t，求从源点到汇点的最大流

const int SIZE=1e5+5;
const int inf=(1ll<<31)-1;
int n,m,s,t,tot,head[SIZE],to[SIZE<<1],nxt[SIZE<<1],flow[SIZE<<1];
inline void add_edge(int a,int b,int c){to[++tot]=b;flow[tot]=c;nxt[tot]=head[a];head[a]=tot;
}
int Q[SIZE],deep[SIZE],cur[SIZE];
bool BFS(){for(int i=1;i<=n;i++)deep[i]=0;for(int i=1;i<=n;i++)cur[i]=head[i];int L=0,R=0;deep[s]=1;Q[++R]=s;while(L<R){int x=Q[++L];for(int i=head[x];i;i=nxt[i]){int y=to[i];if(!deep[y]&&flow[i]>0){deep[y]=deep[x]+1;Q[++R]=y;}}}return deep[t]>0;
}
int dfs(int x,int limit){if(x==t||!limit)return limit;int res=0;for(int &i=cur[x];i;i=nxt[i]){int y=to[i];if(deep[y]==deep[x]+1){int f=dfs(y,min(limit,flow[i]));if(!f)continue;res+=f;limit-=f;flow[i]-=f;flow[i^1]+=f;if(!limit)break;}}return res;
}
int main(){#ifndef ONLINE_JUDGEfreopen("jiedai.in","r",stdin);
#endifmemset(head,0,sizeof(head));tot=1;rd(n),rd(m),rd(s),rd(t);for(int i=1;i<=m;i++){int a,b,c;rd(a),rd(b),rd(c);add_edge(a,b,c);add_edge(b,a,0);}ll mxflow=0;while(BFS())mxflow+=dfs(s,inf);printf("%lld\n",mxflow);return (0-0);
}

最大流最小割定理

最大流=最小割

例题：P2762

m个实验，n个仪器，每个实验有收益，每个仪器有费用，一个实验需要使用若干仪器，购买的仪器可以重复使用，求最大净收益

从源点向每个实验连一条边，容量为实验的收益；
从每个仪器向汇点连一条边，容量为仪器的费用；
从每个实验向每个仪器连一条边，容量为正无穷。
如果要进行一项实验，就割掉它所需要的仪器与汇点之间的边，
这样源点就无法通过这条实验边到达汇点。
如果不进行一项实验，就割掉源点与这个实验之间的边，
同样地，源点也无法通过这条实验边到达汇点。
可以发现，只要源点与汇点是不连通的，就对应存在一组合法的选择方案。
在这个图中，一个割的容量就等于不进行的实验的边权和加上需要的仪器的边权和。
即 割的容量=sum{不要的实验}+sum{要的仪器}割的容量=sum\{不要的实验\}+sum\{要的仪器\}割的容量=sum{不要的实验}+sum{要的仪器}
即 割的容量=sum{所有实验}−sum{要的实验}+sum{要的仪器}割的容量=sum\{所有实验\}-sum\{要的实验\}+sum\{要的仪器\}割的容量=sum{所有实验}−sum{要的实验}+sum{要的仪器}
即 割的容量=sum{所有实验}−净收益割的容量=sum\{所有实验\}-净收益割的容量=sum{所有实验}−净收益
即 净收益=sum{所有实验}−割的容量净收益=sum\{所有实验\}-割的容量净收益=sum{所有实验}−割的容量
那么 最大净收益=sum{所有实验}−最小割最大净收益=sum\{所有实验\}-最小割最大净收益=sum{所有实验}−最小割
根据最大流最小割定理，可得：
最大净收益=sum{所有实验}−最大流最大净收益=sum\{所有实验\}-最大流最大净收益=sum{所有实验}−最大流

费用流（MCMF）

最小费用最大流

const int SIZE=5e4+5;
const int inf=(1ll<<31)-1;
const ll INF=5e18;
int n,m,s,t,all,tot,head[SIZE],to[SIZE<<1],nxt[SIZE<<1],flow[SIZE<<1],cost[SIZE<<1];
void add_edge(int a,int b,int c,int d){to[++tot]=b;flow[tot]=c;cost[tot]=d;nxt[tot]=head[a];head[a]=tot;
}
int Q[SIZE],limit[SIZE],mark[SIZE],pre[SIZE],last[SIZE];
ll dis[SIZE],mxflow,micost;
bool SPFA(){for(int i=1;i<=all;i++)dis[i]=INF,limit[i]=inf,mark[i]=0;int L=0,R=0;Q[R=R%all+1]=s;dis[s]=0;pre[t]=0;while(L!=R){int x=Q[L=L%all+1];mark[x]=0;for(int i=head[x];i;i=nxt[i]){int y=to[i];if(flow[i]&&dis[x]+cost[i]<dis[y]){dis[y]=dis[x]+cost[i];pre[y]=x;last[y]=i;limit[y]=min(limit[x],flow[i]);if(!mark[y]){Q[R=R%all+1]=y;mark[y]=1;}}}}return pre[t]>0;
}
void MCMF(){while(SPFA()){  mxflow+=limit[t];micost+=limit[t]*dis[t];for(int i=t;i!=s;i=pre[i]){flow[last[i]]-=limit[t];flow[last[i]^1]+=limit[t];}}
}
int main(){#ifndef ONLINE_JUDGEfreopen("jiedai.in","r",stdin);
#endifmemset(head,0,sizeof(head));tot=1;rd(n),rd(m),rd(s),rd(t);all=n;for(int i=1;i<=m;i++){int a,b,c,d;rd(a),rd(b),rd(c),rd(d);add_edge(a,b,c,d);add_edge(b,a,0,-d);}MCMF();printf("%lld %lld\n",mxflow,micost);return (0-0);
}

Kruskal 重构树

在运行 Kruskal 最小生成树算法的过程中，当合并 x 和 y 所在的连通块时，若这两个连通块的根节点为 fx 和 fy ，我们新建一个节点 T ，将 fx 和 fy 的父亲指向 T ，同时把 (x,y) 之间的边权当作 T 的点权。这样就得到了 Kruskal 重构树。

Kruskal 重构树是一个大小为 2n-1 的二叉树，且只有非叶子节点有权值。

在最小 Kruskal 重构树中，每个非叶子节点的权值会小于其父亲节点的权值。

例题：n 个点 m 条无向边的图，q 次询问，每次询问从 u 到 v 的所有路径中，最大的边的最小值。BZOJ 3732

分析：建立最小 Kruskal 重构树，然后树上 u 和 v 的 lca 的权值就是答案

const int MAXN=1e5+5;
const int MAXM=2e5+5;
struct EDGE{int x,y,v;bool operator <(const EDGE &A)const{return v<A.v;}
}edge[MAXM];
int n,m,q,FA[MAXN<<1],fa[MAXN<<1][20];
int ls[MAXN<<1],rs[MAXN<<1],val[MAXN<<1],deep[MAXN<<1];
int getfa(int x){if(x==FA[x])return x;return FA[x]=getfa(FA[x]);
}
void init(int x){deep[x]=deep[fa[x][0]]+1;for(int i=1;i<20;i++)fa[x][i]=fa[fa[x][i-1]][i-1];if(x<=n)return;init(ls[x]),init(rs[x]);
}
int LCA(int x,int y){if(deep[x]<deep[y])swap(x,y);int step=deep[x]-deep[y];for(int i=0;i<20;i++)if(step&1<<i)x=fa[x][i];if(x==y)return x;for(int i=19;i>=0;i--)if(fa[x][i]!=fa[y][i])x=fa[x][i],y=fa[y][i];return fa[x][0];
}
int main(){#ifndef ONLINE_JUDGEfreopen("jiedai.in","r",stdin);
#endifrd(n),rd(m),rd(q);for(int i=1;i<=m;i++)rd(edge[i].x),rd(edge[i].y),rd(edge[i].v);sort(edge+1,edge+1+m);for(int i=1;i<n+n;i++)FA[i]=i;int all=n;for(int i=1;i<=m;i++){int x=getfa(edge[i].x),y=getfa(edge[i].y);if(x!=y){val[++all]=edge[i].v;ls[all]=x,rs[all]=y;FA[x]=FA[y]=fa[x][0]=fa[y][0]=all;}}init(all);while(q--){int x,y;rd(x),rd(y);printf("%d\n",val[LCA(x,y)]);}return (0-0);
}

Tarjan

有向图缩点

P3387 以下代码建完图之后会有重边

const int M=1e6+5;
int n,dfsid,all,top,stk[M],low[M],dfn[M],vis[M],bel[M];
vector<int>E[M],RE[M];
void tarjan(int x){dfn[x]=low[x]=++dfsid;stk[++top]=x;vis[x]=1;for(int i=0;i<E[x].size();i++){int y=E[x][i];if(!dfn[y]){tarjan(y);MIN(low[x],low[y]);}else if(vis[y])MIN(low[x],dfn[y]);}if(low[x]==dfn[x]){all++;while(top){int y=stk[top--];bel[y]=all;vis[y]=0;if(y==x)break;}}
}
void rebuild(){dfsid=all=0;for(int i=1;i<=n;i++)dfn[i]=0;for(int i=1;i<=n;i++)if(!dfn[i])tarjan(i);for(int x=1;x<=n;x++)for(int i=0;i<E[x].size();i++){int y=E[x][i];if(bel[x]!=bel[y])RE[bel[x]].push_back(bel[y]);}
}

无向图求割点

P3388

割点的定义：在无向连通图中，如果将其中一个点以及所有连接该点的边去掉，图就不再连通，那么这个点就叫割点

const int M=1e6+5;
int n,m,tot,head[M],to[M<<1],nxt[M<<1];
inline void add_edge(int a,int b){to[++tot]=b;nxt[tot]=head[a];head[a]=tot;
}
int dfsid,dfn[M],low[M],ans[M];
void tarjan(int x,int rt){dfn[x]=low[x]=++dfsid;int son=0;for(int i=head[x];i;i=nxt[i]){int y=to[i];if(!dfn[y]){tarjan(y,rt);MIN(low[x],low[y]);if(x!=rt&&low[y]>=dfn[x])ans[x]=1;if(x==rt)son++;}else MIN(low[x],dfn[y]);}if(son>=2)ans[rt]=1;
}
int main(){#ifndef ONLINE_JUDGEfreopen("jiedai.in","r",stdin);
#endifrd(n),rd(m);for(int i=1;i<=m;i++){int a,b;rd(a),rd(b);add_edge(a,b),add_edge(b,a);}for(int i=1;i<=n;i++)if(dfn[i]==0)tarjan(i,i);int res=0;for(int i=1;i<=n;i++)if(ans[i])res++;printf("%d\n",res);for(int i=1;i<=n;i++)if(ans[i])printf("%d%c",i,--res==0?'\n':' ');return (0-0);
}

无向图求割边

割边的定义：在无向连通图中，如果将一条边去掉，图就不再连通，那么这条边就叫割边

割边的求法与割点类似，若 low[y]>dfn[x]，那么 x 与 y 之间的边就是割边

且不需要考虑根节点的问题

const int M=1e6+5;
int n,m,tot,head[M],to[M<<1],nxt[M<<1];
inline void add_edge(int a,int b){to[++tot]=b;nxt[tot]=head[a];head[a]=tot;
}
int dfsid,dfn[M],low[M];
vector<pair<int,int> >ans;
void tarjan(int x,int f){dfn[x]=low[x]=++dfsid;for(int i=head[x];i;i=nxt[i]){int y=to[i];if(y==f)continue;if(!dfn[y]){tarjan(y,x);MIN(low[x],low[y]);if(low[y]>dfn[x])ans.push_back(make_pair(x,y));}else MIN(low[x],dfn[y]);}
}
int main(){#ifndef ONLINE_JUDGEfreopen("jiedai.in","r",stdin);
#endifrd(n),rd(m);for(int i=1;i<=m;i++){int a,b;rd(a),rd(b);add_edge(a,b),add_edge(b,a);}for(int i=1;i<=n;i++)if(dfn[i]==0)tarjan(i,0);printf("%d\n",ans.size());for(int i=0;i<ans.size();i++)printf("%d %d\n",ans[i].first,ans[i].second);return (0-0);
}

最大团 极大团

const int M=55;
int n,ans,edge[M][M],all[M][M],some[M][M],none[M][M];
void dfs(int d,int an,int sn,int nn){if(!sn&&!nn)MAX(ans,an);int key=some[d][1];for(int i=1;i<=sn;i++){int v=some[d][i],snn=0,nnn=0;if(edge[key][v])continue;for(int j=1;j<=an;j++)all[d+1][j]=all[d][j];all[d+1][an+1]=v;for(int j=1;j<=sn;j++)if(edge[some[d][j]][v])some[d+1][++snn]=some[d][j];for(int j=1;j<=nn;j++)if(edge[none[d][j]][v])none[d+1][++nnn]=none[d][j];dfs(d+1,an+1,snn,nnn);some[d][i]=0,none[d][++nn]=v;}
}
int main(){#ifndef ONLINE_JUDGEfreopen("jiedai.in","r",stdin);
#endifwhile(scanf("%d",&n),n){for(int i=1;i<=n;i++)for(int j=1;j<=n;j++)rd(edge[i][j]);for(int i=1;i<=n;i++)some[1][i]=i;ans=0;dfs(1,0,n,0);printf("%d\n",ans);}return 0;
}

仙人掌 圆方树

P5236 给定一个有 n 个点和 m 条双向边的仙人掌图，q 次询问，每次询问两个点之间的最短路

圆方树的空间要开两倍

const int M=2e5+5;
int n,m,q;
vector<pair<int,ll> >E1[M],E2[M];
int dfsid,all,dfn[M],from[M],fa[M];
int sz[M],deep[M],son[M],top[M];
ll sum[M],dis[M];
void tarjan(int x,int f){fa[x]=f;dfn[x]=++dfsid;for(int i=0;i<E1[x].size();i++){int y=E1[x][i].first,d=E1[x][i].second;if(y==f)continue;if(!dfn[y]){sum[y]=sum[x]+d;tarjan(y,x);if(!from[y]){E2[x].push_back(make_pair(y,d));E2[y].push_back(make_pair(x,d));}}else if(dfn[y]<dfn[x]){sum[++all]=sum[x]-sum[y]+d;for(int k=x;k!=fa[y];k=fa[k]){if(k!=y)from[k]=all;ll tmp=sum[k]-sum[y];E2[k].push_back(make_pair(all,min(tmp,sum[all]-tmp)));E2[all].push_back(make_pair(k,min(tmp,sum[all]-tmp)));}}}
}
void dfs(int x,int f){fa[x]=f;deep[x]=deep[f]+1;sz[x]=1;son[x]=0;for(int i=0;i<E2[x].size();i++){int y=E2[x][i].first;if(y==f)continue;dis[y]=dis[x]+E2[x][i].second;dfs(y,x);sz[x]+=sz[y];if(sz[y]>sz[son[x]])son[x]=y;}
}
void chain(int x,int k){top[x]=k;if(son[x])chain(son[x],k);for(int i=0;i<E2[x].size();i++){int y=E2[x][i].first;if(y==fa[x]||y==son[x])continue;chain(y,y);}
}
int LCA(int a,int b){while(top[a]!=top[b]){if(deep[top[a]]<deep[top[b]])swap(a,b);a=fa[top[a]];}return deep[a]<deep[b]?a:b;
}
int jump(int x,int f){int res;while(top[x]!=top[f]){res=top[x];x=fa[top[x]];}return x==f?res:son[f];
}
int main(){#ifndef ONLINE_JUDGEfreopen("jiedai.in","r",stdin);
#endifrd(n),rd(m),rd(q);while(m--){int a,b,c;rd(a),rd(b),rd(c);E1[a].push_back(make_pair(b,c));E1[b].push_back(make_pair(a,c));}all=n;tarjan(1,0);dfs(1,0);chain(1,1);while(q--){int x,y;rd(x),rd(y);int lca=LCA(x,y);if(lca<=n)printf("%lld\n",dis[x]+dis[y]-2*dis[lca]);else{int a=jump(x,lca),b=jump(y,lca);ll tmp=abs(sum[a]-sum[b]);printf("%lld\n",dis[x]+dis[y]-dis[a]-dis[b]+min(tmp,sum[lca]-tmp));}}return 0;
}

2-sat*

Prufer序列

一个 n 个节点的树可以 n-2 个 [1,n] 之间的数表示，这些数形成的序列叫作 Prufer 序列，树与 Prufer 序列一一对应

无根树转化为 Prufer 序列：

找到编号最小的度数为1的点，删除该节点并在序列中添加与该节点相连的节点的编号，重复这样的过程直到整个树只剩下两个节点

Prufer 序列转化为无根树：

取出 Prufer 序列中最前面的元素 u ，在点集中找到编号最小的没有在 Prufer 序列中出现的元素v，给 u , v 连边，然后在点集中删除 v ，重复这样的过程直到点集中只剩两个点，最后给这两个点连边

Prufer 序列的性质：

1. Cayley 公式：n个点的无向完全图的生成树有 nn−2n^{n-2}nn−2 种
2. Prufer序列中某个点出现的次数等于这个点的度数-1
3. n 个节点的度依次为 d1,d2,...,dnd_1,d_2,...,d_nd1​,d2​,...,dn​ 的树共有 (n−2)!∏i=1n(di−1)!\frac{(n−2)!}{\prod_{i=1}^{n}(d_i−1)!}∏i=1n​(di​−1)!(n−2)!​ 个
4. n 个点的有根树有 nn−1n^{n-1}nn−1 种
5. n 个点 m 条边的无向图，有k个连通块，加 k-1 条边使整个图连通，方案数为 nk−2∏i=1ksin^{k-2}\prod\limits_{i=1}^{k}{s_i}nk−2i=1∏k​si​ （其中 sis_isi​ 为每个连通块的大小）

矩阵树定理

问题：给定一张n个节点，m条双向边的图，求生成树的数量

解法：令A为邻接矩阵， A[i][j]A[i][j]A[i][j] 表示节点 iii 与节点 jjj 之间的边数

D为度数矩阵，D[i][i]D[i][i]D[i][i] 表示节点 iii 的度数，其他无值

则基尔霍夫(Kirchhoff)矩阵为 : K=D−AK=D-AK=D−A

基尔霍夫矩阵去掉任意一行一列的行列式的值，就是生成树的数量

例题：P6178

题意：给定一张n个结点m条边的带权图（图可能是无向图，可能是有向图，有向图的生成树是以1为根的外向树）

定义生成树的权值为生成树中所有边权的乘积，求其所有不同生成树的权值之和，对 1e9+7 取模

解法：将边权为v的边看作v条边权为1的边，问题转化为求生成树个数

对于以1为根的外向树个数，由于是外向树，度数矩阵D中的 D[i][i]D[i][i]D[i][i] 表示 iii 节点的入边权之和（内向树则为出边权之和）

由于以1为根，基尔霍夫矩阵去掉第1行第1列的行列式的值，就是以1为根的生成树的数量（以 x 为根就删掉 x 行 x 列）

const int SIZE=305;
const int P=1e9+7;
int n,m,op,A[SIZE][SIZE];
int fast(int a,int b=P-2){int res=1;while(b){if(b&1)res=1ll*res*a%P;a=1ll*a*a%P;b>>=1;}return res;
}
int det(int n,int A[SIZE][SIZE]){for(int i=1;i<=n;i++)for(int j=1;j<=n;j++)A[i][j]=(A[i][j]%P+P)%P;int res=1,sign=1;for(int i=1;i<=n;i++){int mr=i;for(int j=i+1;j<=n;j++)if(A[j][i]>A[mr][i])mr=j;if(i!=mr)swap(A[i],A[mr]),sign=-sign;res=1ll*res*A[i][i]%P;if(res==0)break;int inv=fast(A[i][i]);for(int j=i+1;j<=n;j++){int d=1ll*A[j][i]*inv%P;for(int k=i;k<=n;k++)A[j][k]=(A[j][k]-1ll*d*A[i][k]%P+P)%P;}}return (res*sign+P)%P;
}
int main(){#ifndef ONLINE_JUDGEfreopen("jiedai.in","r",stdin);
#endifrd(n),rd(m),rd(op);for(int i=1;i<=m;i++){int a,b,c;rd(a),rd(b),rd(c);if(op==0){//无向图A[a][a]=(A[a][a]+c)%P;A[b][b]=(A[b][b]+c)%P;A[a][b]=(A[a][b]-c+P)%P;A[b][a]=(A[b][a]-c+P)%P;}else{//有向图，求以1为根的外向生成树个数A[b][b]=(A[b][b]+c)%P;A[a][b]=(A[a][b]-c+P)%P;}}for(int i=1;i<=n;i++)A[1][i]=A[i][1]=0;A[1][1]=1;printf("%d\n",det(n,A));return (0-0);
}

LGV引理

在一个有向无向图中，给定 n 个起点 A1,A2,...,AnA_1,A_2,...,A_nA1​,A2​,...,An​ 和 n 个终点 B1,B2,...,BnB_1,B_2,...,B_nB1​,B2​,...,Bn​ ，则 n 条简单路径 A1→B1,A2→B2,...,An→BnA_1→B_1,A_2→B_2,...,A_n→B_nA1​→B1​,A2​→B2​,...,An​→Bn​ 互不相交的方案数为

M=∣e(A1,B1)…e(A1,Bn)⋮⋱⋮e(An,B1)…e(An,Bn)∣M=\begin{vmatrix} e(A_1,B_1) & \ldots & e(A_1,B_n) \\ \vdots & \ddots & \vdots \\ e(A_n,B_1) & \ldots & e(A_n,B_n) \\ \end{vmatrix}M=∣∣∣∣∣∣∣​e(A1​,B1​)⋮e(An​,B1​)​…⋱…​e(A1​,Bn​)⋮e(An​,Bn​)​∣∣∣∣∣∣∣​

其中 e(u,v)e(u,v)e(u,v) 表示 u→vu→vu→v 的路径方案数

全局最小割

P5632 Stoer-Wagner算法

题意：给定一个无向连通图，割掉若干条边，使图变成两个连通分量，且边权和最小

考虑图中任意两点 s、t ,若最终 s 与 t 不在同一连通块中，那么使 s 和 t 不连通的最小代价就是 s 和 t 之间的最小割。若最终 s 与 t 在同一连通块中，对于任意一点 i ，如果要割边 s-i 的话，t-i 也要割，如果不割 s-i 的话， t-i 也没必要割，因此，可以把 s 和 t 合并成一个点

于是，找到两个点 s、t，求出这两个点之间的最小割，然后将这两个点合并，重复这样的过程，就能求出全局最小割

构造两个集合 A 和 B ，一开始 A 为空，B 是所有点的集合。定义B 中一个点的权值为这个点到 A 中所有点的边权之和

找某两点的最小割的方法是：不断将 B 中权值最大的点移入 A 中，直到 B 中只剩一个点，令这个点为 t , A 中任意一点为 s ，s 到 t 的最小割为 t 的权值

const int M=605;
int n,m,s,t,del[M],vis[M],ord[M];
ll dis[M][M],val[M];
ll cut(int k){for(int i=1;i<=n;i++)vis[i]=val[i]=0;val[0]=-1;for(int i=1;i<=k;i++){int mx=0;for(int j=1;j<=n;j++)if(!del[j]&&!vis[j]&&val[j]>val[mx])mx=j;vis[mx]=1;ord[i]=mx;for(int j=1;j<=n;j++)if(!del[j]&&!vis[j])val[j]+=dis[mx][j];}s=ord[k-1];t=ord[k];return val[t];
}
ll SW(){ll res=1e18;for(int i=1;i<n;i++){MIN(res,cut(n+1-i));del[t]=1;for(int j=1;j<=n;j++){dis[s][j]+=dis[t][j];dis[j][s]+=dis[j][t];}}return res;
}
int main(){#ifndef ONLINE_JUDGEfreopen("jiedai.in","r",stdin);
#endifrd(n),rd(m);for(int i=1;i<=m;i++){int a,b,c;rd(a),rd(b),rd(c);dis[a][b]+=c;dis[b][a]+=c;}printf("%lld\n",SW());return 0;
}

数学

质数判断

代码一：朴素算法 O(n)O(\sqrt{n})O(n​)

template<class T>bool is_prime(T x){if(x<=3)return x>=2;if(x%2==0||x%3==0)return 0;for(T i=5;i*i<=x;i+=6)if(x%i==0||x%(i+2)==0)return 0;return 1;
}

代码二：Miller_Rabin质数判断

费马小定理： n是一个奇质数，a是任何整数$(1≤a≤n−1) $，则 an−1≡1(modn)a^{n-1}\equiv 1(mod \ n)an−1≡1(mod n)

推论：如果n是一个奇质数，则方程 a2≡1(modn)a^2≡1(mod \ n)a2≡1(mod n) 只有±1两个解，即1和n-1

根据费马小定理，判断n是质数，我们随机选取一个数a，如果不满足an−1≡1(modn)a^{n-1}\equiv 1(mod \ n)an−1≡1(mod n)，则n不是质数

我们将 n−1n-1n−1 拆分成 n−1=d∗2rn-1=d*2^rn−1=d∗2r 的形式，则 an−1=ad∗a2ra^{n-1}=a^d*a^{2^r}an−1=ad∗a2r

也就是说 an−1a^{n-1}an−1 是由 ada^dad 不断进行平方得到的

根据上面的推论，判断n是质数，如果 ada^dad 在不断平方的过程中出现 (ad)2≡1(a^d)^2≡1(ad)2≡1且ad≠±1a^d≠±1ad​=±1 ，则n不是质数

选取足够多个质数进行测试，如果没有出现不符合质数特征的情况，可以认为这个数就是质数

固定选择{2,3,5,7,13,29,37,89}这8个质数作为测试的底数，在long long范围内可以保证质数测试的正确性

inline ll mul(ll a,ll b,ll mod){return (__int128)a*b%mod;
}
inline ll pow(ll a,ll b,ll mod){ll res=1;if(a>=mod)a%=mod;while(b){if(b&1)res=mul(res,a,mod);a=mul(a,a,mod);b>>=1;}return res;
}
bool miller_rabin(ll a,ll n){ll d=n-1,r=0;while(!(d&1))d>>=1,r++;ll x=pow(a,d,n);if(x==1)return 1;for(int i=0;i<r;i++){if(x==n-1)return 1;x=mul(x,x,n);}return 0;
}
bool is_prime(ll n){if(n<=1)return 0;static int num[]={2,3,5,7,13,29,37,89};for(int i=0;i<8;i++)if(n==num[i])return 1;for(int i=0;i<8;i++)if(!miller_rabin(num[i],n))return 0;return 1;
}

质因数分解

代码一：朴素算法 O(n)O(\sqrt{n})O(n​)

template<class T>void prime_factor(T x,vector<T>&Q){Q.clear();while(x%2==0)Q.push_back(2),x/=2;while(x%3==0)Q.push_back(3),x/=3;for(T i=5;i*i<=x;i+=6)if(x%i==0||x%(i+2)==0){while(x%i==0)Q.push_back(i),x/=i;while(x%(i+2)==0)Q.push_back(i+2),x/=i+2;}if(x>1)Q.push_back(x);
}

代码二：Pollard_Rho 算法 O(n14)O(n^\frac{1}{4})O(n41​)

inline ll mul(ll a,ll b,ll mod){return (__int128)a*b%mod;
}
inline ll pow(ll a,ll b,ll mod){ll res=1;if(a>=mod)a%=mod;while(b){if(b&1)res=mul(res,a,mod);a=mul(a,a,mod);b>>=1;}return res;
}
bool miller_rabin(ll a,ll n){ll d=n-1,r=0;while(!(d&1))d>>=1,r++;ll x=pow(a,d,n);if(x==1)return 1;for(int i=0;i<r;i++){if(x==n-1)return 1;x=mul(x,x,n);}return 0;
}
bool is_prime(ll x){if(x<=1)return 0;static int num[]={2,3,5,7,13,29,37,89};for(int i=0;i<8;i++)if(x==num[i])return 1;for(int i=0;i<8;i++)if(!miller_rabin(num[i],x))return 0;return 1;
}
ll fun(ll x,ll c,ll mod){return (mul(x,x,mod)+c)%mod;
}
ll gcd(ll n,ll m){if(m==0)return n;return gcd(m,n%m);
}
ll pollard_rho(ll x){ll c=rand()%(x-1)+1;ll s=0,t=0;for(int goal=1;;goal<<=1,s=t){ll val=1;for(int step=1;step<=goal;step++){t=fun(t,c,x);val=mul(val,abs(s-t),x);if(step%127==0){ll d=gcd(val,x);if(d>1)return d;}}ll d=gcd(val,x);if(d>1)return d;}
}
void find_fac(ll x){if(x==1)return;if(is_prime(x)){//找到了质因子xreturn;}ll y=x;while(y==x)y=pollard_rho(x);find_fac(y),find_fac(x/y);
}

组合数

组合数常见变换：

1. (nm)=(nn−m)\displaystyle \binom{n}{m} = \binom{n}{n-m}(mn​)=(n−mn​)
2. (nk)=nk(n−1k−1)\displaystyle \binom{n}{k} = \frac{n}{k} \binom{n-1}{k-1}(kn​)=kn​(k−1n−1​)，递推式
3. (nm)=(n−1m)+(n−1m−1)\displaystyle \binom{n}{m} = \binom{n-1}{m} + \binom{n-1}{m-1}(mn​)=(mn−1​)+(m−1n−1​)，杨辉三角的公式表达，O(n2)O(n^2)O(n2) 推导组合数
4. ∑i=0n(ni)=2n\displaystyle \sum_{i=0}^{n}\binom{n}{i} = 2^ni=0∑n​(in​)=2n，二项式定理 (a+b)n(a+b)^n(a+b)n 中 a=1,b=1a=1,b=1a=1,b=1 的特殊情况
5. ∑i=0n(−1)i(nm)=0\displaystyle \sum_{i=0}^{n}(-1)^i \binom{n}{m} = 0i=0∑n​(−1)i(mn​)=0，二项式定理中 a=1,b=−1a=1,b=-1a=1,b=−1 的特殊情况
6. ∑i=0n(ni)(mi)=∑i=0n(ni)(mm−i)=(n+mm)(n≥m)\displaystyle \sum_{i=0}^{n} \binom{n}{i} \binom{m}{i} = \sum_{i=0}^{n} \binom{n}{i} \binom{m}{m-i} = \binom{n+m}{m} \ \ (n≥m)i=0∑n​(in​)(im​)=i=0∑n​(in​)(m−im​)=(mn+m​)  (n≥m)
7. ∑i=0n(ni)2=(2nn)\displaystyle \sum_{i=0}^{n}\binom{n}{i}^2 = \binom{2n}{n}i=0∑n​(in​)2=(n2n​)，上式中 n=mn=mn=m 的特殊情况
8. ∑i=0ni(ni)=n2n−1\displaystyle \sum_{i=0}^{n}i\binom{n}{i} = n2^{n-1}i=0∑n​i(in​)=n2n−1，带权和的一个式子
9. ∑i=0ni2(ni)=n(n+1)2n−2\displaystyle \sum_{i=0}^{n} i^2\binom{n}{i} = n(n+1)2^{n-2}i=0∑n​i2(in​)=n(n+1)2n−2，带权和的另一个式子
10. ∑i=0n(ik)=(n+1k+1)\displaystyle \sum_{i=0}^{n}\binom{i}{k} = \binom{n+1}{k+1}i=0∑n​(ki​)=(k+1n+1​)
11. (nr)(rk)=(nk)(n−kr−k)\displaystyle \binom{n}{r} \binom{r}{k} = \binom{n}{k} \binom{n-k}{r-k}(rn​)(kr​)=(kn​)(r−kn−k​)
12. ∑i=0n(ia)(n−ib−a)=(n+1b+1)\displaystyle \sum_{i=0}^{n}\binom{i}{a} \binom{n-i}{b-a} = \binom{n+1}{b+1}i=0∑n​(ai​)(b−an−i​)=(b+1n+1​)，与 aaa 无关，a∈[0,b]a \in [0,b]a∈[0,b]
13. ∑i=0n(n−ii)=Fn+1\displaystyle \sum_{i=0}{n}\binom{n-i}{i} = F_{n+1}i=0∑​n(in−i​)=Fn+1​，其中 FFF 是斐波拉契数列

数论分块

向下取整

求 ∑i=1n⌊ni⌋\sum\limits_{i=1}^{n}{⌊\frac{n}{i}⌋}i=1∑n​⌊in​⌋

ll calc(int n){ll res=0;for(int l=1,r;l<=n;l=r+1){r=n/(n/l);res+=(r-l+1)*(n/l);}return res;
}

向上取整

求 ∑i=1n⌈ni⌉\sum\limits_{i=1}^{n}{⌈\frac{n}{i}⌉}i=1∑n​⌈in​⌉

ll calc(int n){ll res=0;for(int l=1,r;l<=n;l=r+1){if(l!=n)r=(n-1)/((n-1)/l);else r=n;res+=(r-l+1)*((n+l-1)/l);//n+l-1可能会爆int}return res;
}

线性筛

SIZE>70000时，质数的个数<SIZE/10

筛质数

pre[i]表示i最小的质数，判断质数：pre[i]==i（1需要特判）

const int SIZE=1e7+5;
int ptot,prime[SIZE/10],pre[SIZE];
void init_prime(){for(int i=2;i<SIZE;i++){if(!pre[i])prime[++ptot]=pre[i]=i;for(int j=1;j<=ptot;j++){int t=i*prime[j];if(t>=SIZE)break;pre[t]=prime[j];if(i%prime[j]==0)break;}}
}

欧拉函数

定义：不大于n的正整数中与n互质的数的个数，用 ϕ(i)\phi(i)ϕ(i) 表示，ϕ(1)=1\phi{(1)}=1ϕ(1)=1

通式：ϕ(n)=n∗∏i=1xpi−1pi\phi(n)=n*\prod\limits_{i=1}^{x}{\frac{p_i-1}{p_i}}ϕ(n)=n∗i=1∏x​pi​pi​−1​ ，其中p1,p2...pxp_1,p_2...p_xp1​,p2​...px​为n的所有不同质因子

性质：欧拉函数是积性函数，若m，n互质，则ϕ(nm)=ϕ(n)∗ϕ(m)\phi(nm)=\phi(n)*\phi(m)ϕ(nm)=ϕ(n)∗ϕ(m)

当n为奇数时，ϕ(2n)=ϕ(n)\phi(2n)=\phi(n)ϕ(2n)=ϕ(n)

当n为质数时，ϕ(n)=n−1,Xϕ(n)≡1(modn)\phi(n)=n-1,X^{\phi(n)}\equiv 1(mod \ n)ϕ(n)=n−1,Xϕ(n)≡1(mod n) （n，x互质）

还有一个性质，∑d∣nϕ(d)=n\sum\limits_{d|n}{\phi(d)}=nd∣n∑​ϕ(d)=n

欧拉定理：

若 gcd(a,p)=1gcd(a,p)=1gcd(a,p)=1 ， aϕ(p)≡1(modp)a^{\phi(p)}\equiv 1(mod \ p)aϕ(p)≡1(mod p)

扩展欧拉定理：
ab≡{abmodϕ(p)gcd(a,p)=1abgcd(a,p)≠1,b<ϕ(p)a(bmodϕ(p))+ϕ(p)gcd(a,p)≠1,b≥ϕ(p)(modp)a^b \equiv \begin{cases} a^{b \bmod \phi(p)} & gcd(a,p)=1 \\ a^{b} & gcd(a,p)\neq 1,b<\phi(p) \\ a^{(b \bmod \phi(p))+\phi(p)} & gcd(a,p)\neq 1,b\geq \phi(p) \end{cases} \pmod p ab≡⎩⎪⎨⎪⎧​abmodϕ(p)aba(bmodϕ(p))+ϕ(p)​gcd(a,p)=1gcd(a,p)​=1,b<ϕ(p)gcd(a,p)​=1,b≥ϕ(p)​(modp)

const int SIZE=1e7+5;
int ptot,prime[SIZE/10],phi[SIZE];
bool vis[SIZE];
void init_phi(){phi[1]=1;for(int i=2;i<SIZE;i++){if(!vis[i])prime[++ptot]=i,phi[i]=i-1;for(int j=1;j<=ptot;j++){int t=i*prime[j];if(t>=SIZE)break;vis[t]=1;if(i%prime[j]==0){phi[t]=phi[i]*prime[j];break;}phi[t]=phi[i]*phi[prime[j]];}}
}

莫比乌斯函数

const int SIZE=1e7+5;
int ptot,prime[SIZE/10],mul[SIZE];
bool vis[SIZE];
void init_mul(){mul[1]=1;for(int i=2;i<SIZE;i++){if(!vis[i])prime[++ptot]=i,mul[i]=-1;for(int j=1;j<=ptot;j++){int t=i*prime[j];if(t>=SIZE)break;vis[t]=1;if(i%prime[j]==0){mul[t]=0;break;}mul[t]=-mul[i];}}
}

杜教筛*

Min_25筛*

扩展欧几里得定理

给定 a,ba,ba,b ，求 ax+by=gcd(a,b)ax+by=gcd(a,b)ax+by=gcd(a,b) 的整数解

设 ax1+by1=gcd(a,b)ax_1+by_1=gcd(a,b)ax1​+by1​=gcd(a,b)

bx2+(amodb)y2=gcd(b,amodb)bx_2+(a \ mod \ b)y_2=gcd(b,a \ mod \ b)bx2​+(a mod b)y2​=gcd(b,a mod b)

根据欧几里得定理：gcd(a,b)=gcd(b,amodb)gcd(a,b)=gcd(b,a \ mod \ b)gcd(a,b)=gcd(b,a mod b)

ax1+by1=bx2+(amodb)y2ax_1+by_1=bx_2+(a \ mod \ b)y_2ax1​+by1​=bx2​+(a mod b)y2​

ax1+by1=bx2+(a−⌊ab⌋∗b)y2ax_1+by_1=bx_2+(a-\lfloor \frac{a}{b} \rfloor *b)y_2ax1​+by1​=bx2​+(a−⌊ba​⌋∗b)y2​

ax1+by1=ay2+b(x2−⌊ab⌋y2)ax_1+by_1=ay_2+b(x_2-\lfloor \frac{a}{b} \rfloor y_2)ax1​+by1​=ay2​+b(x2​−⌊ba​⌋y2​)

因此：x1=y2,y1=x2−⌊ab⌋y2x_1=y_2,y_1=x_2-\lfloor \frac{a}{b} \rfloor y_2x1​=y2​,y1​=x2​−⌊ba​⌋y2​

将 x,yx,yx,y 不断代入递归求解，直至 b=0b=0b=0 时， x=1,y=0x=1,y=0x=1,y=0 是 ax+by=gcd(a,b)ax+by=gcd(a,b)ax+by=gcd(a,b) 的一组解

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,x,y);ll tmp=x;x=y,y=tmp-(a/b)*y;return d;
}

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

中国剩余定理

求解下列同余方程组（模数互质）

{x≡a1(modm1)x≡a2(modm2)⋮x≡an(modmn)\begin{cases} x \equiv a_1 \ (\ mod \ m_1\ )\\ x \equiv a_2 \ (\ mod \ m_2\ )\\ \quad \vdots\\ x \equiv a_n \ (\ mod \ m_n) \end{cases}⎩⎪⎪⎪⎪⎨⎪⎪⎪⎪⎧​x≡a1​ ( mod m1​ )x≡a2​ ( mod m2​ )⋮x≡an​ ( mod mn​)​

构造法
令：lcm=∏i=1nmi\displaystyle lcm=\prod\limits_{i=1}^{n}{m_i}lcm=i=1∏n​mi​，Mi=lcmmi\displaystyle M_i=\frac{lcm}{m_i}Mi​=mi​lcm​，Mi∗ti≡1(modmi)M_i*t_i \equiv 1 \ (mod \ m_i)Mi​∗ti​≡1 (mod mi​)

则：res=∑i=1nai∗Mi∗ti\displaystyle res=\sum\limits_{i=1}^{n}{a_i*M_i*t_i}res=i=1∑n​ai​∗Mi​∗ti​

struct CRT{vector<ll>A,B;//res%B[i]==A[i]void add(ll a,ll b){A.push_back(a);B.push_back(b);}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;}ll solve(){ll lcm=1,res=0;for(int i=0;i<A.size();i++)lcm*=B[i];for(int i=0;i<A.size();i++){ll tmp=lcm/B[i],x,y;exgcd(tmp,B[i],x,y);x=(x%B[i]+B[i])%B[i];res=(res+(__int128)A[i]*x%lcm*tmp)%lcm;}return res;}
};

拓展中国剩余定理

求解下列同余方程组（模数不一定互质）

{x≡a1(modm1)x≡a2(modm2)⋮x≡an(modmn)\begin{cases} x \equiv a_1 \ (\ mod \ m_1\ )\\ x \equiv a_2 \ (\ mod \ m_2\ )\\ \quad \vdots\\ x \equiv a_n \ (\ mod \ m_n) \end{cases}⎩⎪⎪⎪⎪⎨⎪⎪⎪⎪⎧​x≡a1​ ( mod m1​ )x≡a2​ ( mod m2​ )⋮x≡an​ ( mod mn​)​

假设已经求出前 k−1k-1k−1 个方程组成的同余方程组的一个解为 xxx ，且令M=LCMi=1k−1miM=LCM_{i=1}^{k-1}{m_i}M=LCMi=1k−1​mi​

则前 k−1k-1k−1 个方程的方程组通解为 x+i∗M(i∈Z)x+i*M(i∈Z)x+i∗M(i∈Z)

那么加入第 kkk 个方程后，我们需要求一个正整数 ttt，使得 x+t∗M≡ak(modmk)x+t*M \equiv a_k \ (mod \ m_k)x+t∗M≡ak​ (mod mk​)

即 t∗M≡ak−x(modmk)t*M \equiv a_k-x \ (mod \ m_k)t∗M≡ak​−x (mod mk​) ，可以用拓展欧几里得求出 ttt

若该同余式无解，则整个方程组无解

若有解，则前 kkk 个同余式组成的方程组的一个解为 xk=x+t∗Mx_k=x+t*Mxk​=x+t∗M

struct EXCRT{vector<ll>A,B;//res%B[i]==A[i]void add(ll a,ll b){A.push_back(a);B.push_back(b);}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;}ll solve(){ll res=0,lcm=1;for(int i=0;i<A.size();i++){ll x,y,g=exgcd(lcm,B[i],x,y),tmp=B[i]/g;if((A[i]-res)%g)return -1;exgcd(lcm/g,tmp,x,y);x=(((__int128)(A[i]-res)/g*x)%tmp+tmp)%tmp;lcm*=tmp;res=((res+(__int128)x*(lcm/tmp))%lcm+lcm)%lcm;}return res;}
};

卢卡斯定理

Cab≡C⌊aP⌋⌊bP⌋⋅CamodPbmodP(modP)C_{a}^{b} \equiv C_{\lfloor \frac{a}{P} \rfloor }^{\lfloor \frac{b}{P} \rfloor }·C_{a \ mod \ P}^{b \ mod \ P} \ (mod \ P)Cab​≡C⌊Pa​⌋⌊Pb​⌋​⋅Ca mod Pb mod P​ (mod P)

模数P为质数，建议在 10710^7107 以内

struct Lucas{int P;vector<int>fac,inv;Lucas(int p=2):P(p){fac.resize(P);inv.resize(P);fac[0]=fac[1]=inv[0]=inv[1]=1;for(int i=2;i<P;i++)fac[i]=1ll*fac[i-1]*i%P;for(int i=2;i<P;i++)inv[i]=1ll*(P-P/i)*inv[P%i]%P;for(int i=2;i<P;i++)inv[i]=1ll*inv[i]*inv[i-1]%P;}int comb(ll n,ll m){if(n<0||m<0||n<m)return 0;if(n<P)return 1ll*fac[n]*inv[m]%P*inv[n-m]%P;return 1ll*comb(n/P,m/P)*comb(n%P,m%P)%P;}
};

原根

原根的定义：m 是正整数，g 是整数，若 g 模 m 的阶等于 φ(m)\varphi(m)φ(m) ，则称 g 是 m 的原根

原根存在定理：一个数 m 存在原根当且仅当 m=2,4,pa,2pa（p是奇质数，a是正整数）m=2,4,p^a,2p^a（p是奇质数，a是正整数）m=2,4,pa,2pa（p是奇质数，a是正整数）

原根判定定理：若 g 是 m 的原根，则对于任意 k∣φ(m)，k≠1，k≠φ(m)k|\varphi(m)，k≠1，k≠\varphi(m)k∣φ(m)，k​=1，k​=φ(m) ，都有 gφ(m)/k≢1(modm)g^{\varphi(m)/k}\not\equiv 1 (mod \ m)gφ(m)/k​≡1(mod m)

求所有原根：若 g 是 m 的原根，对于所有 1≤s≤φ(m),gcd(s,(φ(m))=11≤s≤\varphi(m),gcd(s,(\varphi(m))=11≤s≤φ(m),gcd(s,(φ(m))=1，gsg^sgs 为 m 的所有原根

原根个数：若 m 有原根，则它的原根个数为 φ(φ(m))\varphi(\varphi(m))φ(φ(m))

template<class T>void prime_factor(T x,vector<T>&Q){Q.clear();while(x%2==0)Q.push_back(2),x/=2;while(x%3==0)Q.push_back(3),x/=3;for(T i=5;i*i<=x;i+=6)if(x%i==0||x%(i+2)==0){while(x%i==0)Q.push_back(i),x/=i;while(x%(i+2)==0)Q.push_back(i+2),x/=i+2;}if(x>1)Q.push_back(x);
}
int min_G(int n){vector<int>Q;prime_factor(n,Q);if(n==2||n==4);else if(n%2==1&&Q.size()>=1&&Q[0]==Q.back());else if(n%4==2&&Q.size()>=2&&Q[1]==Q.back());else return -1;prime_factor(n-1,Q);Q.erase(unique(Q.begin(),Q.end()),Q.end());for(int i=1;i<n;i++){bool flag=1;for(auto j:Q){int res=1,a=i,b=(n-1)/j;while(b){if(b&1)res=1ll*res*a%n;a=1ll*a*a%n;b>>=1;}if(res==1){flag=0;break;}}if(flag)return i;}
}

二次剩余

定义：如果关于x的方程 x2≡n(modp)x^2 \equiv n \ (mod \ p)x2≡n (mod p) 有解

则称是n模p意义下的二次剩余，否则称为模意p意义下的二次非剩余

(np)(\frac{n}{p})(pn​) 表示n为关于p的勒让德符号：

若 (np)=0(\frac{n}{p})=0(pn​)=0 ，则p是n的倍数

若 (np)=1(\frac{n}{p})=1(pn​)=1 ，则n是模p意义下的二次剩余

若 (np)=−1(\frac{n}{p})=-1(pn​)=−1 ，则n是模p意义下的二次非剩余

模数为奇质数

欧拉准则：(np)≡np−12(modp)(\frac{n}{p}) \equiv n^{\frac{p-1}{2}} \ (mod \ p)(pn​)≡n2p−1​ (mod p)

Cipolla算法：求解 x2≡n(modp)x^2 \equiv n \ (mod \ p)x2≡n (mod p)

随机找到一个数a，使得 (a2−np)=−1(\frac{a^2-n}{p})=-1(pa2−n​)=−1 ，即 （a2−n）p−12≡p−1(modp)（a^2-n）^{\frac{p-1}{2}} \equiv p-1 \ (mod \ p)（a2−n）2p−1​≡p−1 (mod p)

可以证明，随机找a的期望次数为2

定义 w2≡a2−nw^2 \equiv a^2-nw2≡a2−n ，由于 a2−na^2-na2−n 不能开方，那么 www 就类似于一个虚数单位

这样的复数可以用 a+bwa+bwa+bw 来表示，简写成 (a,b)(a,b)(a,b)

(a,b)∗(c,d)≡(ac+bdw2,ad+bc)(a,b)*(c,d) \equiv (ac+bdw^2,ad+bc)(a,b)∗(c,d)≡(ac+bdw2,ad+bc)

结论：上述方程的一个解为 x≡(a+w)p+12x \equiv (a+w)^{\frac{p+1}{2}}x≡(a+w)2p+1​ ，另外一个解为 p−xp-xp−x

int w,P;
struct CP{int x,y;CP(int x=0,int y=0):x(x),y(y){}CP operator *(const CP &A)const{return CP((1ll*x*A.x+1ll*y*A.y%P*w)%P,(1ll*x*A.y+1ll*y*A.x)%P);}
};
int fast(int a,int b=P-2){int res=1;while(b){if(b&1)res=1ll*res*a%P;a=1ll*a*a%P;b>>=1;}return res;
}
CP fast(CP a,int b){CP res=1;while(b){if(b&1)res=res*a;a=a*a;b>>=1;}return res;
}
int Cipolla(int n){//x and P-xif(n==0)return 0;if(P==2)return n;if(fast(n,(P-1)/2)==P-1)return -1;int a;while(1){a=rand()%P;w=(1ll*a*a-n+P)%P;if(fast(w,(P-1)/2)==P-1)break;}return fast(CP(a,1),(P+1)/2).x;
}

模数为奇质数的幂次

模数为2的幂次

模数任意

BSGS*

斐波那契数列

定义：F0=0,F1=1,Fn=Fn−1+Fn−2F_0=0,F_1=1,F_n=F_{n-1}+F_{n-2}F0​=0,F1​=1,Fn​=Fn−1​+Fn−2​

通项公式：Fn=15[(1+52)n−(1−52)n]F_n=\frac{1}{\sqrt5}[(\frac{1+\sqrt5}{2})^n-(\frac{1-\sqrt5}{2})^n]Fn​=5​1​[(21+5​​)n−(21−5​​)n]

tip：5≡383008016或616991993(mod1000000009)\sqrt{5} \equiv 383008016 \ 或 \ 616991993 \ (mod \ 1000000009)5​≡383008016 或 616991993 (mod 1000000009)

性质：

1. Fn+m+1=FnFm+Fn+1Fm+1F_{n+m+1}=F_nF_m+F_{n+1}F_{m+1}Fn+m+1​=Fn​Fm​+Fn+1​Fm+1​

   推论：F2n+1=Fn2+Fn+12F_{2n+1}=F_n^2+F_{n+1}^2F2n+1​=Fn2​+Fn+12​

   推论：F2n=Fn+12−Fn−12=2FnFn+1−Fn2F_{2n}=F_{n+1}^2-F_{n-1}^2=2F_nF_{n+1}-F_n^2F2n​=Fn+12​−Fn−12​=2Fn​Fn+1​−Fn2​

2. Fn−kFn+k=Fn2+(−1)n+kFk2F_{n-k}F_{n+k}=F_n^2+(-1)^{n+k}F_k^2Fn−k​Fn+k​=Fn2​+(−1)n+kFk2​

   推论：Fn−1Fn+1=Fn2+(−1)nF_{n-1}F_{n+1}=F_n^2+(-1)^nFn−1​Fn+1​=Fn2​+(−1)n

3. F1+F2+...+Fn=Fn+2−1F_1+F_2+...+F_n=F_{n+2}-1F1​+F2​+...+Fn​=Fn+2​−1

4. F12+F22+...+Fn2=FnFn+1F_1^2+F_2^2+...+F_n^2=F_n F_{n+1}F12​+F22​+...+Fn2​=Fn​Fn+1​

5. F1+F3+...+F2n−1=F2nF_1+F_3+...+F_{2n-1}=F_{2n}F1​+F3​+...+F2n−1​=F2n​

6. F2+F4+...+F2n=F2n+1−1F_2+F_4+...+F_{2n}=F_{2n+1}-1F2​+F4​+...+F2n​=F2n+1​−1

7. F1+2F2+3F3+...+nFn=nFn+2−Fn+3+2F_1+2F_2+3F_3+...+nF_n=nF_{n+2}-F_{n+3}+2F1​+2F2​+3F3​+...+nFn​=nFn+2​−Fn+3​+2

8. −F1+F2−F3+...+(−1)nFn=(−1)nFn−1−1-F_1+F_2-F_3+...+(-1)^nF_n=(-1)^nF_{n-1}-1−F1​+F2​−F3​+...+(−1)nFn​=(−1)nFn−1​−1

9. 3Fn=Fn−2+Fn+23F_n=F_{n-2}+F_{n+2}3Fn​=Fn−2​+Fn+2​

10. F2n−2m−2(F2n+F2n+2)=F2m+2+F4n−2mF_{2n-2m-2}(F_{2n}+F_{2n+2})=F_{2m+2}+F_{4n-2m}F2n−2m−2​(F2n​+F2n+2​)=F2m+2​+F4n−2m​ (n>m≥−1,n≥1)(n>m≥-1,n≥1)(n>m≥−1,n≥1)

11. Fn∣Fm⇔n∣mF_n|F_m⇔n|mFn​∣Fm​⇔n∣m

12. gcd(Fn,Fm)=Fgcd(n,m)gcd(F_n,F_m)=F_{gcd(n,m)}gcd(Fn​,Fm​)=Fgcd(n,m)​

循环节：

Fn%2F_n\%2Fn​%2 的循环节为3，Fn%3F_n\%3Fn​%3 的循环节为8，Fn%5F_n\%5Fn​%5 的循环节为20

Fn%10F_n\%10Fn​%10 的循环节为60， Fn%100F_n\%100Fn​%100 的循环节为300，Fn%1000F_n\%1000Fn​%1000 的循环节为1500

Fn%PF_n\%PFn​%P 的循环节不超过6P

利用性质1的推论求FnF_nFn​：

fib(n)返回 (Fn,Fn+1)(F_n,F_{n+1})(Fn​,Fn+1​) 的二元组，复杂度为 O(logn)O(logn)O(logn)

pair<int,int> fib(ll n){//fib(n).firstif(n==0)return make_pair(0,1);pair<int,int>tmp=fib(n>>1);int a=tmp.first,b=tmp.second;int c=(2ll*b-a+P)*a%P,d=(1ll*a*a+1ll*b*b)%P;if(n&1)return make_pair(d,(c+d)%P);else return make_pair(c,d);
}

广义斐波那契数列

定义一：f1=1,f2=a+b,fn=(a+b)fn−1−abfn−2f_1=1,f_2=a+b,f_n=(a+b)f_{n-1}-abf_{n-2}f1​=1,f2​=a+b,fn​=(a+b)fn−1​−abfn−2​

通项公式：fn=an−bna−bf_n=\frac{a^n-b^n}{a-b}fn​=a−ban−bn​

当 a+b=1,ab=−1a+b=1,ab=-1a+b=1,ab=−1 时即为斐波那契数列

定义二：f0=a,f1=b,fn=fn−1+fn−2f_0=a,f_1=b,f_n=f_{n-1}+f_{n-2}f0​=a,f1​=b,fn​=fn−1​+fn−2​

通项公式：fn=Fn−1a+Fnbf_n=F_{n-1}a+F_{n}bfn​=Fn−1​a+Fn​b

当 a=0,b=1a=0,b=1a=0,b=1 时即为斐波那契数列

矩阵

矩阵加法、矩阵乘法

const int SIZE=305;
const int P=1e9+7;
int n;
struct matrix{int a[SIZE][SIZE];matrix(){memset(a,0,sizeof(a));}matrix operator +(const matrix &A)const{matrix T;for(int i=1;i<=n;i++)for(int j=1;j<=n;j++)T.a[i][j]=(a[i][j]+A.a[i][j])%P;return T;}matrix operator *(const matrix &A)const{matrix T;for(int i=1;i<=n;i++)for(int k=1;k<=n;k++)for(int j=1;j<=n;j++)T.a[i][j]=(T.a[i][j]+1ll*a[i][k]*A.a[k][j])%P;return T;}
};

矩阵快速幂

matrix fast(matrix a,ll b){matrix res;for(int i=1;i<=n;i++)res.a[i][i]=1;while(b){if(b&1)res=res*a;a=a*a;b>>=1;}return res;
}

矩阵等比数列求和

计算 A1+A2+...+AkA^1+A^2+...+A^kA1+A2+...+Ak

matrix calc(matrix A,ll k){static int stk[105],top;matrix E,res=A,tmp=A;if(k==0)return E;for(int i=1;i<=n;i++)E.a[i][i]=1;while(k)stk[++top]=k%2,k/=2;for(top--;top;top--){res=(E+tmp)*res;tmp=tmp*tmp;if(stk[top]){tmp=tmp*A;res=res+tmp;}}return res;
}

矩阵求逆

int fast(int a,int b=P-2){int res=1;while(b){if(b&1)res=1ll*res*a%P;a=1ll*a*a%P;b>>=1;}return res;
}
bool inverse(matrix &mat,matrix &res){static int A[SIZE][SIZE<<1];for(int i=1;i<=n;i++)for(int j=1;j<=n;j++)A[i][j]=mat.a[i][j],A[i][j+n]=(i==j);for(int i=1;i<=n;i++){for(int j=i;j<=n;j++)if(A[j][i]){swap(A[i],A[j]);break;}if(A[i][i]==0)return 0;int inv=fast(A[i][i]);for(int k=i;k<=n+n;k++)A[i][k]=1ll*A[i][k]*inv%P;for(int j=1;j<=n;j++)if(i!=j)for(int k=n+n;k>=i;k--)A[j][k]=(A[j][k]-1ll*A[j][i]*A[i][k]%P+P)%P;}for(int i=1;i<=n;i++)for(int j=1;j<=n;j++)res.a[i][j]=A[i][j+n];return 1;
}

伴随矩阵

设矩阵 A=(aij)n✖nA=(a_{ij})_{n✖n}A=(aij​)n✖n​ ，将矩阵 AAA 中的元素 aija_{ij}aij​ 所在的第 i 行第 j 列元素去掉后，剩余的各元素组成的n-1阶矩阵所确定的行列式称为元素 aija_{ij}aij​ 的余子式，记为 MijM_{ij}Mij​ 。

称 Aij=(−1)i+jMijA_{ij}=(-1)^{i+j}M_{ij}Aij​=(−1)i+jMij​ 为元素 aija_{ij}aij​ 的代数余子式。

矩阵 AAA 的伴随矩阵 A∗A^{*}A∗ 的构造如下：

[A11A21...An1A12A22...An2⋮⋮⋮A1nA2n...Ann]\left[ \begin{matrix} A_{11} & A_{21} & ... & A_{n1} \\ A_{12} & A_{22} & ... & A_{n2} \\ \vdots & \vdots & & \vdots \\ A_{1n} & A_{2n} & ... & A_{nn} \end{matrix} \right]⎣⎢⎢⎢⎡​A11​A12​⋮A1n​​A21​A22​⋮A2n​​.........​An1​An2​⋮Ann​​⎦⎥⎥⎥⎤​

如果 AAA 可逆，则 A∗=∣A∣A−1A^{*}=|A|A^{-1}A∗=∣A∣A−1

高斯消元

代码一：判断解的情况，并求唯一解

返回值：0表示无解，1表示存在唯一解，2表示存在无穷多解

如果存在唯一解，系数矩阵会被化简为n*n的单位矩阵，则常数列就是对应的解

判断无解要在化简结束之后，要检查所有行是否出现一边为零一边不为零的情况

const int SIZE=305;
const double eps=1e-10;
int Gauss(int n,double A[SIZE][SIZE]){bool Inf=0;for(int i=1;i<=n;i++){for(int j=i;j<=n;j++)if(fabs(A[j][i])>=eps){swap(A[i],A[j]);break;}if(fabs(A[i][i])<eps){Inf=1;continue;}for(int k=n+1;k>=i;k--)A[i][k]/=A[i][i];for(int j=1;j<=n;j++)if(j!=i)for(int k=n+1;k>=i;k--)A[j][k]-=A[j][i]*A[i][k];}for(int i=1;i<=n;i++){int flag=1;for(int j=1;j<=n;j++)if(fabs(A[i][j])>=eps){flag=0;break;}if(flag&&fabs(A[i][n+1])>=eps)return 0;}if(Inf)return 2;else return 1;
}

行列式求值

行列式的性质：

• 矩阵转置，行列式不变
• 矩阵行（列）交换，行列式取反
• 矩阵行（列）相加或相减，行列式不变
• 矩阵行（列）所有元素同时乘以数 ，行列式等比例变大

对矩阵进行高斯消元，得到一个对角线矩阵，矩阵的行列式就是对角线元素之积。其符号由交换行的数量来确定（如果为奇数，则行列式的符号应颠倒）

代码一：模数为质数，复杂度为O(n3)O(n^3)O(n3)

const int SIZE=305;
const int P=1e9+7;
int fast(int a,int b=P-2){int res=1;while(b){if(b&1)res=1ll*res*a%P;a=1ll*a*a%P;b>>=1;}return res;
}
int det(int n,int A[SIZE][SIZE]){for(int i=1;i<=n;i++)for(int j=1;j<=n;j++)A[i][j]=(A[i][j]%P+P)%P;int res=1,sign=1;for(int i=1;i<=n;i++){int mr=i;for(int j=i+1;j<=n;j++)if(A[j][i]>A[mr][i])mr=j;if(i!=mr)swap(A[i],A[mr]),sign=-sign;res=1ll*res*A[i][i]%P;if(res==0)break;int inv=fast(A[i][i]);for(int j=i+1;j<=n;j++){int d=1ll*A[j][i]*inv%P;for(int k=i;k<=n;k++)A[j][k]=(A[j][k]-1ll*d*A[i][k]%P+P)%P;}}return (res*sign+P)%P;
}

代码二：任意模数，复杂度为O(n3logP)O(n^{3}logP)O(n3logP)

const int SIZE=305;
int det(int n,int A[SIZE][SIZE],int P){int res=1;for(int i=1;i<=n;i++){for(int j=i+1;j<=n;j++)while(A[j][i]){int t=A[i][i]/A[j][i];for(int k=i;k<=n;k++)A[i][k]=(A[i][k]-1ll*A[j][k]*t)%P;swap(A[i],A[j]);res=-res;}res=1ll*res*A[i][i]%P;if(!res)return 0;}return (res+P)%P;
}

线性基

线性基的定义

由原集合AAA得到线性基TTT

使得TTT中元素互相异或所形成的集合，等价于原序列AAA的元素互相异或形成的集合

可以理解为线性基将原序列进行了压缩

线性基的性质

1. 线性基能相互异或得到原集合的所有相互异或得到的值
2. 线性基是满足性质111的最小的集合
3. 线性基没有异或和为000的子集

线性基的构造

设一个数组ddd，表示序列aaa的线性基，下标从000开始

d[i]表示的是线性基中第i+1i+1i+1个位置上所存储的数字

如果这个数字不为000

那么这个数字转化为二进制数后的第i+1i+1i+1位是111

且这个二进制数的最高位是第i+1i+1i+1位

我们都将序列aaa中每一个数依次插入线性基中，得到了序列aaa的线性基

线性基的插入

从下面代码可知，在线性基中插入数valvalval时，从高位到低位依次扫描它为111的二进制位

当扫描到第i位时，如果d[i]d[i]d[i]不存在，就令d[i]=vald[i]=vald[i]=val，插入结束

如果d[i]存在，此时valvalval与d[i]d[i]d[i]的第i+1i+1i+1位都为111，就令valvalval=valvalval^d[i]d[i]d[i]

于是valvalval与d[i]d[i]d[i]异或后第i+1i+1i+1位变成了000，然后继续向下扫描

最终，valvalval会有插入成功和插入不成功这两种结局

如果插入成功，就说明当前线性基里的一些数异或起来不能等于valvalval，因此valvalval是不可替代的

如果插入不成功，是因为当前线性基里的一些数异或起来可以等于valvalval

并且valvalval在一系列操作之后变成了000，因此valvalval是多余的

bool insert(ll *d,ll val){//向线性基中插入一个数for(ll i=59;i>=0;i--)if(val&1ll<<i){if(d[i])val^=d[i];else{d[i]=val;break;}}return val>0;//判断val是否插入成功
}

查询最大异或和

准确地说，是求一个序列中的若干个数的异或和的最大值

从高位到低位扫描线性基，如果答案异或d[i]d[i]d[i]后变大了，就让答案异或d[i]d[i]d[i]

这其实是个贪心的过程，我们只需让答案的高位尽可能大

当扫到d[i]d[i]d[i]且d[i]d[i]d[i]不为000时

由于d[i]d[i]d[i]的第i+1i+1i+1位是111，且d[i]的i+1i+1i+1位以上都是000

所以如果答案的第i+1i+1i+1位是000，则答案异或上d[i]d[i]d[i]之后一定会变大

如果答案的第i+1i+1i+1位是111，则答案异或上d[i]d[i]d[i]之后一定会变小

ll query_max(ll *d){ll res=0;for(ll i=59;i>=0;i--)if((res^d[i])>res)res^=d[i];return res;
}

查询最小异或和

显然，如果这个线性基有无法插入的数，那么最小异或和就为000

否则，最小值是线性基中最小的d[i]d[i]d[i]

这是因为如果让最小的d[i]d[i]d[i]去异或其它的d[i]d[i]d[i]，那么它一定会变大，所以它自己就是最小的

查询第k小异或和

准确地说：从一个序列中取任意个元素进行异或，求能异或出的所有数字中第kkk小的值。

要求第kkk小值，首先将线性基进行改造，改造后每一个d[i]d[i]d[i]相互独立

对于每一个d[i]d[i]d[i]，枚举jjj=iii to 111，如果d[i]d[i]d[i]的第jjj位为1，那么让d[i]d[i]d[i]异或d[j−1]d[j−1]d[j−1]

这样改造后，如果d[i]d[i]d[i]不为000，那么所有d[j]d[j]d[j]的第i+1i+1i+1位上为111的只有d[i]d[i]d[i]

于是线性基中的元素，作用其实都是提供自己最高位上的111

那么只要使提供出来的111可以和kkk的每一位上的111对应

那么求出来的答案就是第kkk小的

void rebuild(ll *d){for(int i=59;i>=0;i--)for(int j=i-1;j>=0;j--)if(d[i]&1ll<<j)d[i]^=d[j];
}
ll query_kth(ll *d,ll k){rebuild(d);int cnt=0;for(int i=0;i<=59;i++)if(d[i])cnt++;if(cnt<n&&k==1)return 0;//最小异或和为0的情况（n为向线性基中插入的次数）if(cnt<n)k--;//最小异或和为0的情况（n为向线性基中插入的次数）if(k>=(1ll<<cnt))return -1;//不存在第k小异或和ll ans=0;for(int i=0;i<=59;i++)if(d[i]){if(k&1)ans^=d[i];k>>=1;}return ans;
}

线性基求交

要得到线性基A与B的交集，构造一个线性基C，首先，C的初始值为A。

遍历线性基B，当遍历到元素x时：

1. 若元素x不能被线性基C表示，说明x与线性基C中的元素线性无关（即这些元素们之间不能互相表示），那么就将元素x插入到线性基C中。

2. 若元素x能被线性基C表示，即线性基C中的若干个元素的异或和为x，那么这些元素中必定有线性基A中的元素，也可能有线性基B中的元素，然后将那些只属于线性基A的元素加入到答案线性基中。（或者只将线性基B中的元素插入也行）

void merge(ll *res,ll *A,ll *B){static ll C[60],D[60],tmp[60];for(int i=0;i<60;i++){C[i]=A[i];if(A[i])D[i]=1ll<<i;tmp[i]=0;}for(int i=0;i<60;i++)if(B[i]){ll x=B[i],k=0,flag=1;for(int j=59;j>=0;j--){if(x&1ll<<j){if(C[j])x^=C[j],k^=D[j];else{C[j]=x,D[j]=k;flag=0;break;}}}if(flag){ll x=0;for(int j=0;j<60;j++)if(k&1ll<<j)x^=C[j];insert(tmp,x);}}for(int i=0;i<60;i++)res[i]=tmp[i];
}

前缀线性基

HDOJ 6579：有n个数，m次操作，强制在线

操作0 l r：询问[l,r]的最大异或和

操作1 x：序列的最后添加一个数x

暴力的做法可以用数据结构维护区间线性基，但肯定过不了。贪心地维护序列的前缀线性基 (上三角形态)，对于每个线性基，将出现位置靠右的数 字尽可能地放在高位，也就是说在插入新数字的时候，要同时记录对应位置上数字的出现位 置，并且在找到可以插入的位置的时候，如果新数字比位置上原来的数字更靠右，就将该位 置上原来的数字向低位推。 在求最大值的时候，从高位向低位遍历，如果该位上的数字出现在询问中区间左端点的 右侧且可以使答案变大，就异或到答案里。 对于线性基的每一位，与它异或过的线性基更高位置上的数字肯定都出现在它右侧 (否 则它就会被插入在那个位置了)，因此做法的正确性显然。

const int SIZE=1e6+5;
int n,m,cas,A[SIZE];
int d[SIZE][30],id[SIZE][30];
void build(int val,int pos){int tmp=pos;for(int i=29;i>=0;i--){d[pos][i]=d[pos-1][i];id[pos][i]=id[pos-1][i];}for(int i=29;i>=0;i--)if(val&1<<i){if(d[pos][i]){if(tmp>id[pos][i]){swap(d[pos][i],val);swap(id[pos][i],tmp);}val^=d[pos][i];}else{d[pos][i]=val;id[pos][i]=tmp;break;}}
}
int main(){#ifndef ONLINE_JUDGEfreopen("jiedai.in","r",stdin);
#endifrd(cas);while(cas--){rd(n),rd(m);for(int i=1;i<=n;i++)rd(A[i]),build(A[i],i);int ans=0,op,x,l,r;while(m--){rd(op);if(op)rd(A[++n]),build(A[n]^=ans,n);else{rd(l),rd(r);l=(l^ans)%n+1,r=(r^ans)%n+1;if(l>r)swap(l,r);ans=0;for(int i=29;i>=0;i--)if(id[r][i]>=l)MAX(ans,ans^d[r][i]);printf("%d\n",ans);}}for(int i=1;i<=n;i++)for(int j=0;j<30;j++)d[i][j]=id[i][j]=0;}return (0-0);
}

拉格朗日插值*

FFT

多项式乘法，可跑 1e6​ 的数据

const int SIZE=(1<<21)+5;
const double PI=acos(-1);
struct CP{double x,y;CP(double x=0,double y=0):x(x),y(y){}CP operator +(const CP &A)const{return CP(x+A.x,y+A.y);}CP operator -(const CP &A)const{return CP(x-A.x,y-A.y);}CP operator *(const CP &A)const{return CP(x*A.x-y*A.y,x*A.y+y*A.x);}
};
int limit,rev[SIZE];
void DFT(CP *F,int op){for(int i=0;i<limit;i++)if(i<rev[i])swap(F[i],F[rev[i]]);for(int mid=1;mid<limit;mid<<=1){CP wn(cos(PI/mid),op*sin(PI/mid));for(int len=mid<<1,k=0;k<limit;k+=len){CP w(1,0);for(int i=k;i<k+mid;i++){CP tmp=w*F[i+mid];F[i+mid]=F[i]-tmp;F[i]=F[i]+tmp;w=w*wn;}}}if(op==-1)for(int i=0;i<limit;i++)F[i].x/=limit;
}
void FFT(int n,int m,CP *F,CP *G){for(limit=1;limit<=n+m;limit<<=1);for(int i=0;i<limit;i++)rev[i]=(rev[i>>1]>>1)|((i&1)?limit>>1:0);DFT(F,1),DFT(G,1);for(int i=0;i<limit;i++)F[i]=F[i]*G[i];DFT(F,-1);
}int n,m;
CP F[SIZE],G[SIZE];
int main(){#ifndef ONLINE_JUDGEfreopen("jiedai.in","r",stdin);
#endifscanf("%d%d",&n,&m);for(int i=0;i<=n;i++)scanf("%lf",&F[i].x);for(int i=0;i<=m;i++)scanf("%lf",&G[i].x);FFT(n,m,F,G);for(int i=0;i<=n+m;i++)printf("%d%c",int(F[i].x+0.5),i==n+m?'\n':' ');return (0-0);
}

NTT

jxl板子：

const int maxn = (1<<19)+7;
const int mod = 998244353, G = 3, sqr2 = 116195171;
int n, m, nm, up;
ll f[maxn], g[maxn];
int tr[maxn];ll qpow(ll a,ll k=mod-2,ll p=mod){ll s=1;while(k){if(k&1)s=s*a%p;a=a*a%p;k>>=1;}return s;}void NTT(ll *f, int op) {for(int i=0; i<up; ++i) if(i<tr[i]) swap(f[i],f[tr[i]]);for(int len=2; len<=up; len<<=1) {ll w1=qpow(qpow(G,(mod-1)/len),~op?1:mod-2);for(int l=0, hf=len>>1; l<up; l+=len) {ll w=1;for(int i=l; i<l+hf; ++i) {ll tp=w*f[i+hf]%mod;f[i+hf]=(f[i]-tp+mod)%mod;f[i]=(f[i]+tp)%mod;w=w*w1%mod;}}}if(op==-1) for(int i=0, inv=qpow(up,mod-2); i<up; ++i) f[i]=f[i]*inv%mod;
}int main() {n=read(), m=read();for(int i=0; i<=n; ++i) f[i]=read();for(int i=0; i<=m; ++i) g[i]=read();for(nm=n+m, up=1; up<=nm; up<<=1);for(int i=0; i<up; ++i) tr[i]=(tr[i>>1]>>1)|((i&1)?up>>1:0);NTT(f,1); NTT(g,1);for(int i=0; i<up; ++i) f[i]=f[i]*g[i]%mod;NTT(f,-1);ll ans=0;for(int i=0; i<=nm; ++i) ans=(ans+f[i])%mod;printf("%lld\n", ans);
}

FWT

FWT(f,3)，就是子集和。

jxl板子：

const int maxn = (1<<20)+7;
const int mod = 998244353, inv2 = 499122177;int n, up;
ll f[maxn], g[maxn], h[maxn];void FWT(ll *f, int op) {for(int len=2; len<=up; len<<=1) {for(int l=0, hf=len>>1; l<up; l+=len) {for(int i=l; i<l+hf; ++i) {ll x=f[i], y=f[i+hf];if(op>0) {if(op==1) f[i]=(x+y)%mod, f[i+hf]=(x-y+mod)%mod; //xorelse if(op==2) f[i]=(x+y)%mod; //andelse f[i+hf]=(x+y)%mod; //or}else {if(op==-1) f[i]=(x+y)*inv2%mod, f[i+hf]=(x-y+mod)*inv2%mod; //xorelse if(op==-2) f[i]=(x-y+mod)%mod; //andelse f[i+hf]=(y-x+mod)%mod; //or}}}}
}int main() {for(int i=0; i<up; ++i) f[i]=read();for(int i=0; i<up; ++i) g[i]=read();FWT(f,3); FWT(g,3);for(int i=0; i<up; ++i) h[i]=f[i]*g[i]%mod;FWT(f,-3); FWT(g,-3); FWT(h,-3);for(int i=0; i<up; ++i) printf("%lld%c", h[i], " \n"[i==up-1]);FWT(f,2); FWT(g,2);for(int i=0; i<up; ++i) h[i]=f[i]*g[i]%mod;FWT(f,-2); FWT(g,-2); FWT(h,-2);for(int i=0; i<up; ++i) printf("%lld%c", h[i], " \n"[i==up-1]);FWT(f,1); FWT(g,1);for(int i=0; i<up; ++i) h[i]=f[i]*g[i]%mod;FWT(f,-1); FWT(g,-1); FWT(h,-1);for(int i=0; i<up; ++i) printf("%lld%c", h[i], " \n"[i==up-1]);
}

生成函数*

自适应辛普森积分

利用二次函数的积分来近似需要的积分，表达式为：∫lrf(x)dx=(r−l)∗(f(l)+f(r)+4f(l+r2))6\int_l^r{f(x)dx}=\frac{(r-l)*(f(l)+f(r)+4f(\frac{l+r}{2}))}{6}∫lr​f(x)dx=6(r−l)∗(f(l)+f(r)+4f(2l+r​))​

利用递归计算积分，根据误差大小判断是否继续划分

调用时：double res=Simpson(l,r,simpson(l,r),1e-6);

inline double fun(double x){return (0.0);//需要积分的函数
}
inline double simpson(double l,double r){return (r-l)*(fun(l)+fun(r)+4*fun((l+r)/2))/6;
}
double Simpson(double l,double r,double res,double eps){//注意函数名的大小写double mid=(l+r)/2;double L=simpson(l,mid);double R=simpson(mid,r);if(fabs(L+R-res)<=15*eps)return L+R+(L+R-res)/15;return Simpson(l,mid,L,eps/2)+Simpson(mid,r,R,eps/2);
}

数学杂项

卡特兰数

1, 1, 2, 5, 14, 42, 132, 429, 1430, 4862, 16796, 58786, …
C1=1,C2=1Cn+1=C0Cn+C1Cn−1+...+CnC0(n−3)Cn=n2(C3Cn−1+C4Cn−2+...+Cn−2C4+Cn−1C3)hn=cn+1,h0=1,h1=1hn=hn−1(4n−2)n+1hn=C2nnn+1=C2nn−C2nn−1C_1=1, C_2=1 \\ C_{n+1}=C_0C_n+C_1C_{n-1}+...+C_nC_0 \\ (n-3)C_n=\frac{n}{2}(C_3C_{n-1}+C_4C_{n-2}+...+C_{n-2}C_4+C_{n-1}C_3) \\ h_n=c_{n+1}, h_0=1,h_1=1 \\ h_n=\frac{h_{n-1}(4n-2)}{n+1} \\ h_n=\frac{C_{2n}^n}{n+1}=C_{2n}^n-C_{2n}^{n-1} C1​=1,C2​=1Cn+1​=C0​Cn​+C1​Cn−1​+...+Cn​C0​(n−3)Cn​=2n​(C3​Cn−1​+C4​Cn−2​+...+Cn−2​C4​+Cn−1​C3​)hn​=cn+1​,h0​=1,h1​=1hn​=n+1hn−1​(4n−2)​hn​=n+1C2nn​​=C2nn​−C2nn−1​

约瑟夫环问题

n个人从1开始编号，站成一圈，从1号开始依次报数，报到k的倍数的人出局，问最后剩下的人的编号

代码一：线性算法 O(n)

int josephus(int n,int k){int res=0;for(int i=1;i<=n;i++)res=(res+k)%i;return res+1;
}

代码二：对数算法 O(k⋅logn)O(k·logn)O(k⋅logn)

int josephus(int n,int k){if(n==1)return 1;if(k==1)return n;if(k>n)return (josephus(n-1,k)-1+k)%n+1;int res=josephus(n-n/k,k)-n%k;if(res<1)res+=n;else res+=(res-1)/(k-1);return res;
}

已知年月日，求星期几

int day(int y,int m,int d){if(m<=2)m+=12,y--;return (d+2*m+3*(m+1)/5+y+y/4-y/100+y/400)%7+1;
}

[1,2n]内选a个奇数b个偶数互不相邻的方案数

ans=Cn−baCn−abans=C_{n-b}^{a}C_{n-a}^{b}ans=Cn−ba​Cn−ab​

求区间[l,r]内经过点x的长度小于等于k的线段个数

ll calc(ll l,ll r,ll x,ll k){//l<=i<=x<=j<=r&&j-i+1<=kll res=0;ll a=min(x+1,max(l,x+1-k));ll b=max(l-1,min(x,r+1-k));if(b>=a)res=(b+k-x)*(b+k-x+1)/2-(a+k-x-1)*(a+k-x)/2;a=min(x+1,max(l,r-k+2));b=x;if(b>=a)res+=(b-a+1)*(r-x+1);return res;
}

博弈论

SG函数*

Nim Game

尼姆博弈：n堆石子，每次可以从一堆中取出任意个石子，不能取的一方失败

一堆石子的SG值为其石子数，整个游戏的SG值是所有石子堆的SG值的异或和

必胜：SG不为0，必败：SG为0

Bash Game

巴什博弈：n堆石子，每次可以从一堆中取出最多k个石子，不能取的一方失败

一堆石子的SG值为其石子数 mod (k+1)，整个游戏的SG值是所有石子堆的SG值的异或和

必胜：SG不为0，必败：SG为0

Nim-k Game

Moore’s Nimk博弈：n堆石子，每次可以从最多k堆中每堆取出任意个石子，不能取的一方失败

一堆石子的SG值为其石子数，对于每个二进制位，求出有多少堆石子的SG值在该位上为1

对于每一个二进制位，若该位上为1的数量 mod (k+1)都为0，则为必败，否则必胜

Anti-Nim Game

反Nim博弈：n堆石子，每次可以从一堆中取出任意个石子，不能取的一方获胜

一堆石子的SG值为其石子数，整个游戏的SG值是所有石子堆的SG值的异或和

必胜：(1)SG为0且每一堆的石子数都不超过1

​ (2)SG不为0且至少有一堆石子数大于1

必败：其余为必败

Anti-SG Game

反SG博弈：SG游戏中不能行动的一方获胜

SJ定理前提：对于任意一个 Anti-SG 游戏，当局面中所有的单一游戏的SG值都为0时，游戏结束

必胜：(1)游戏的SG值为0且每一个单一游戏的SG值都不超过1

​ (2)游戏的SG值不为0且至少有一个单一游戏的SG值大于1

必败：其余为必败

Staircase Game

阶梯博弈：n层阶梯，每层阶梯上有一堆石子，每次可以从从一层阶梯上取出任意个石子放到下一层阶梯，不能取的一方失败

整个游戏的SG值为奇数阶梯上的石子数的异或和

如果移动偶数层的石子到奇数层，对手一定可以继续移动这些石子到偶数层，使得局面的SG不变

必胜：SG不为0，必败：SG为0

Multi-SG Game

Multi-SG游戏规定，在符合拓扑原则的前提下，一个单一游戏的后继可以为多个单一游戏

我们可以通过异或运算把多个单一游戏连接到一起，作为一个后继状态

Lasker’s Nim博弈：n堆石子，每次可以从一堆中取出任意个石子，或者将一堆石子分成两堆非空石子，不能行动的一方失败

一堆石子的SG值为：SG(0)=0，SG(4k+1)=4k+1，SG(4k+2)=4k+2，SG(4k+3)=4k+4，SG(4k+4)=4k+3

整个游戏的SG值是所有石子堆的SG值的异或和

必胜：SG不为0，必败：SG为0

拓展Lasker’s Nim博弈：n堆石子，每次可以从一堆中取出任意个石子，或者将一堆石子分成三堆非空石子，不能行动的一方失败

一堆石子的SG值为：SG(0)=0，SG(8k+x)=8k+x (1≤x≤6)，SG(8k+7)=8k+8，SG(8k+8)=8k+7

Wythoff Game

威佐夫博弈：两堆石子，每次可以从一堆中取出任意个石子，或者同时从两堆中取出相同数量的石子，不能行动的一方失败

必败态为(1,2),(3,5),(4,7),(6,10),(8,13)…

(1) 每对数的第一个数为前面没有出现过的最小正整数

(2) 每对数的第二个数为第一个数加上项数

Betty定理：α,β为正无理数且 1α+1β=1\frac{1}{α}+\frac{1}{β}=1α1​+β1​=1，则数列 {⌊αn⌋} 与 {⌊βn⌋} 无交集，且能恰好覆盖正整数集合

令 an=⌊αn⌋,bn=⌊βn⌋a_n=⌊αn⌋,b_n=⌊βn⌋an​=⌊αn⌋,bn​=⌊βn⌋ ，根据(2)，⌊βn⌋=⌊αn⌋+n=⌊(α+1)n⌋，所以 β=α+1，则 1α+1α+1=1\frac{1}{α}+\frac{1}{α+1}=1α1​+α+11​=1 ，解得 α=1+52α=\frac{1+\sqrt{5}}{2}α=21+5​​

必败态为 an=⌊1+52n⌋,bn=an+na_n=⌊\frac{1+\sqrt{5}}{2}n⌋,b_n=a_n+nan​=⌊21+5​​n⌋,bn​=an​+n

拓展威佐夫博弈：两堆石子，每次可以从一堆中取出任意个石子，或者同时从两堆中取出数量差小于等于k的石子，不能行动的一方失败

必败态为 an=⌊1−k+k2+2k+52n⌋,bn=an+(k+1)na_n=⌊\frac{1-k+\sqrt{k^2+2k+5}}{2}n⌋,b_n=a_n+(k+1)nan​=⌊21−k+k2+2k+5​​n⌋,bn​=an​+(k+1)n

ll calc_an(int n,int k){return (sqrt((long double)k*k+k+k+5)+1-k)*0.5*n;
}
ll calc_bn(int n,int k){return calc_an(n,k)+1ll*(k+1)*n;
}

Every-SG Game

问题模型：n个游戏同时进行，对于还没有结束的所有单一游戏，游戏者必须对其进行决策，没有单一游戏可决策的一方失败

对于先手方来说，必胜的单一游戏，玩得越久越好；必败的单一游戏，越早结束越好

对于必胜的单一游戏，需要知道最慢几步可以结束游戏；对于必败的单一游戏，需要知道最快几步可以结束游戏

令step表示结束游戏的步数

若状态S是结束状态，那么 step(s)=0

若状态S是先手必胜态，那么 step(s)=max（step(t)）+1，其中 t 是先手必败的后继状态

若状态S是先手必败态，那么 step(s)=min（step(t)）+1

先手必胜当且仅当所有单一游戏中最大的 step 为奇数

Fibonacci-Nim Game

斐波拉契博弈：一堆石子，大小为n，先手第一次取不能取完，之后每次取的数量不能超过对手刚取的数量的两倍，不能取的一方失败

(1) 如果 n 是斐波拉契数，设为 F(k) ，则 F(k)=F(k-1)+F(k-2) ，若先手取的石子数大于等于 F(k-2) ,那么后手可以把剩下的石子全部取完；那么先手只能在 F(k-2) 中取，这是一个子问题，用归纳法可以证明这是必败的，剩下的 F(k-1) 也是个子问题，此时先手必败

(2) 如果 n 不是斐波拉契数，根据Zeckendorf定理：正整数都可以表示成若干个不连续的斐波那契数之和，先手可以取完分解后最小的一堆石子，后手肯定取不完剩下的最小的一堆石子，参考(1)，最小的一堆会被先手取完，之后的每一堆，先手都能取完，此时先手必胜

结论：先手必败当且仅当 n 是斐波拉契数，先手必胜的策略是取分解之后最小的斐波拉契数

k倍动态减法博弈

一堆石子，大小为n，先手第一次取不能取完，之后每次取的数量不能超过对手刚取的数量的k倍，不能取的一方失败

这是斐波拉契博弈的拓展问题。与斐波拉契博弈相似，我们要求一个数列，将 n 分解成数列中一些项的和

我们用 a 数组表示要求的数列，b[i] 表示 a[0…i] 能够构造的最大数字，满足这个数字分解之后的项在 a[0…i] 中

(1) a[i]=b[i-1]+1

(2) b[i]=a[i]+b[t] 其中 t 是最大的整数满足 a[t]*k<=a[i]

结论：先手必败当且仅当 n 是 a 数列中的数，先手必胜的策略是取分解之后最小的项

int calc(int n,int k){//n<=1e8 k<=1e5static const int M=750005;static int A[M],B[M];int res,i,t=0;A[1]=B[1]=1;for(i=2;;i++){A[i]=B[i-1]+1;if(A[i]>=n)break;while(A[t+1]*k<A[i])t++;B[i]=A[i]+B[t];}while(A[i]>n)i--;if(A[i]==n)return -1;while(n){while(A[i]>n)i--;n-=A[i];res=A[i];}return res;
}

树上删边游戏

在n个点的有根树中，每次可以删掉一条边以及删边后不再与边相连的部分，不能行动的一方失败

叶子节点的SG值为0，非叶子节点的SG值为其所子节点 (SG值+1)的异或和

必胜：根节点SG不为0，必败：根节点SG为0

无向图删边游戏

在n个点的无向联通图中，每次可以删掉一条边以及删边后不再与边相连的部分，不能行动的一方失败

将任意一个偶环缩成一个新点

将任意一个奇环缩成一个新点，并在该点上加一条新边

所有原来连到环上的边改为与新点相连，问题转化为树上删边游戏

翻硬币游戏

问题模型：n个硬币排成一排，双方根据某些约束轮流翻硬币，满足翻转的硬币中，最右边的必须是从正面翻到反面，不能翻的一方失败

结论：局面的SG值等于局面中每个正面朝上的硬币单一存在时的SG值的异或和。

每次只能翻转一个硬币：SG(0)=0,SG(n)=1SG(0)=0,SG(n)=1SG(0)=0,SG(n)=1

每次只能翻转一个或两个硬币：SG(n)=nSG(n)=nSG(n)=n

每次只能翻转一个或两个或三个硬币（Mock Turtles Game）：SG(n)=2n−1−popcount(n−1)%2SG(n)=2n-1-popcount(n-1)\%2SG(n)=2n−1−popcount(n−1)%2

每次只能翻转连续的k个硬币：SG(n)=[n%k=0]SG(n)=[n\%k=0]SG(n)=[n%k=0]

每次只能翻转连续的任意个硬币（Ruler Game）：SG(n)=lowbit(n)SG(n)=lowbit(n)SG(n)=lowbit(n)

每次只能翻转两个硬币，翻转的硬币的距离不超过k：SG(n)=(n−1)%(k+1)SG(n)=(n-1)\%(k+1)SG(n)=(n−1)%(k+1)

高维组合游戏

Nim积：Nim积与Nim和的关系就类似于乘法与加法

Tartan定理：对于一个高维组合游戏，玩家的操作是一维操作的笛卡尔积的形式，那么这个游戏的SG值为每一维单独的SG值的Nim积

例如HDU3404：在一个n+1行m+1列的矩阵中，坐标从(0,0)到(n,m)，每个格子上有一个硬币，正面朝上或反面朝上。双方轮流操作，每次操作必须翻转一个矩形四个角上的硬币，而且其中行列坐标都较大的硬币是从正面翻到反面。不能操作的一方失败。

1. 局面的SG值等于局面中每个正面朝上的硬币单一存在时的SG值的Nim和
2. 每个硬币单独存在的SG值等于其行列坐标的Nim积

计算Nim积的复杂度为 O(log2x)O(log^{2}{x})O(log2x)

int Nim_multi_power(int x,int y){//x=1<<(1<<i)if(x<=1||y<=1)return x*y;int m;for(int i=0;1<<(1<<i)<=x;i++)m=1<<(1<<i);int p=x/m,s=y/m,t=y%m;int d1=Nim_multi_power(p,s);int d2=Nim_multi_power(p,t);return (m*(d1^d2))^Nim_multi_power(m/2,d1);
}
int Nim_multi(int x,int y){if(x<=1||y<=1)return x*y;if(x<y)swap(x,y);int m;for(int i=0;1<<(1<<i)<=x;i++)m=1<<(1<<i);int p=x/m,q=x%m,s=y/m,t=y%m;int c1=Nim_multi(p,s);int c2=Nim_multi(p,t)^Nim_multi(q,s);int c3=Nim_multi(q,t);return (m*(c1^c2))^c3^Nim_multi_power(m/2,c1);
}

不平等博弈

超现实数（Surreal Number）：在数学中，超现实数是一个包含无穷大数、无穷小数还有实数的域。一个超现实数由左集合与右集合组成，记作 x={L|R} , 其中 x 大于 L 中任意一个元素且小于 R 中任意一个元素，且满足 x=A2Bx=\frac{A}{2^B}x=2BA​ ， A 是整数，B 是非负整数，在 B 尽可能小的情况下 A 的绝对值尽可能小

在一个组合游戏中，玩家被称为左玩家和右玩家。一个单一游戏的局面可以用一个超现实数SN表示，超现实数的左集合为左玩家操作一次后所有局面的集合，右集合为右玩家操作一次后所有局面的集合。SN值可以理解为左玩家比右玩家多操作的次数。整个游戏局面的SN值是所有单一游戏SN值之和

结论：若 SN>0，则左玩家必胜；若 SN<0，则右玩家必胜；若 SN=0，则后手方必胜

Blue-Red Hackenbush

若干个 BW 串，玩家 W 只能拿 W，玩家 B 只能拿 B，每次拿走一个字符后其后缀也会消失，不能操作的一方输

例如串WWWBWWB，其SN值为 3−2−1+2−2+2−3−2−4=45163-2^{-1}+2^{-2}+2^{-3}-2^{-4}=\frac{45}{16}3−2−1+2−2+2−3−2−4=1645​

计算每个串的SN值并求和，若 SN>0 则 W 必胜；若 SN<0 则 B 必胜；若 SN=0 则后手必胜

以下代码计算超现实数，用于打表找规律

double surreal(double l,double r){if(l<0&&r>0)return 0;else if(floor(l+1)>l&&floor(l+1)<r){if(l>=0)return floor(l+1);else return ceil(r-1);}else{for(int i=1;i<1e9;i*=2)for(int j=floor(l*i+1);j<r*i;j++)if(j>l*i&&j<r*i)return 1.0*j/i;}
}

其他博弈

擦数游戏

n个数1~n，每次选一个数划去它的所有约数，不能行动的一方失败

先手必胜

字符串*

动态规划

区间覆盖计数问题

n个位置，m个区间，求选出的区间能够覆盖[1,n]的方案数

分析：设 dp[i]dp[i]dp[i] 表示恰好覆盖了区间 [1,i][1,i][1,i] 的方案数，设 sum[i]=∑k=0idp[k]sum[i]=\sum\limits_{k=0}^{i}{dp[k]}sum[i]=k=0∑i​dp[k]

将区间按照右端点排序，从小到大扫描每个右端点

若只有一个右端点为 iii 的区间 [l,i][l,i][l,i] ，那么 dp[i]=∑k=l−1i−1dp[k]=sum[i−1]−sum[l−2]dp[i]=\sum\limits_{k=l-1}^{i-1}{dp[k]}=sum[i-1]-sum[l-2]dp[i]=k=l−1∑i−1​dp[k]=sum[i−1]−sum[l−2]

若有多个右端点相同的区间，则将这些区间再按照左端点从小到大排序

对于这些右端点相同的区间，枚举第一个选取的区间，由于按左端点排序后排在前面的区间会包含排在后面的区间，那么当一个区间选取后，排在后面的区间无论是否选取都可以

于是 dp[i]=∑j=ST(sum[i−1]−sum[lj−2])⋅2T−jdp[i]=\sum\limits_{j=S}^{T}{(sum[i-1]-sum[l_j-2])·2^{T-j}}dp[i]=j=S∑T​(sum[i−1]−sum[lj​−2])⋅2T−j

复杂度为 O(n+mlogm)O(n+mlogm)O(n+mlogm)

const int SIZE=1e6+5;
const int P=1e9+7;
int solve(int n,int m,pair<int,int>*Q){static int dp[SIZE],sum[SIZE],two[SIZE];sort(Q+1,Q+1+m,[&](pair<int,int>a,pair<int,int>b){if(a.second!=b.second)return a.second<b.second;else return a.first<b.first;});dp[0]=sum[0]=two[0]=1;for(int i=1;i<SIZE;i++)two[i]=2ll*two[i-1]%P;int S=1,T=0;for(int i=1;i<=n;i++){dp[i]=0;if(S<=m&&Q[S].second==i){while(T<m&&Q[T+1].second==i)T++;for(int j=S;j<=T;j++){int tmp=(Q[j].first>=2?sum[Q[j].first-2]:0);dp[i]=(dp[i]+1ll*(sum[i-1]-tmp+P)*two[T-j])%P;}S=T+1;}sum[i]=(sum[i-1]+dp[i])%P;}return dp[n];
}

仙人掌DP

BZOJ 4316

题意：求仙人掌的最大独立集，即选出最多的点，这些点之间互相没有边连接

dp[i][j][k]dp[i][j][k]dp[i][j][k] 表示在 i 的子树中，j 为 i 点是否选择，k 为 i 到 i 的父亲这条边所在的环的最低点是否选择，的最大值

根据当前点是否是一个环的最低点，来初始化 DP 数组

转移时，分三种情况讨论：桥边、别的环的边、当前环的边

const int M=1e5+5;
int n,m,tot,head[M],to[M<<1],nxt[M<<1];
inline void add_edge(int a,int b){to[++tot]=b;nxt[tot]=head[a];head[a]=tot;
}
int dfsid,all,fa[M],dfn[M],from[M],bot[M];
int dp[M][2][2],tmp[2][2];
void dfs(int x,int f){fa[x]=f;dfn[x]=++dfsid;for(int i=head[x];i;i=nxt[i]){int y=to[i];if(y==f)continue;if(!dfn[y])dfs(y,x);else if(dfn[y]<dfn[x]){all++;bot[x]=1;for(int k=x;k!=y;k=fa[k])from[k]=all;}}//以下为DP部分if(bot[x])dp[x][1][1]=1;else dp[x][1][0]=1;for(int i=head[x];i;i=nxt[i]){int y=to[i];if(fa[y]!=x)continue;memset(tmp,0,sizeof(tmp));for(int j=0;j<2;j++)for(int k=0;k<2;k++)for(int a=0;a<2;a++)for(int b=0;b<2;b++){if(j&&a)continue;if(from[y]==0){//桥边MAX(tmp[j][k],dp[x][j][k]+dp[y][a][b]);}else if(from[y]!=from[x]){//别的环的边if(j&&b)continue;MAX(tmp[j][k],dp[x][j][k]+dp[y][a][b]);}else if(from[y]==from[x]){//当前环的边MAX(tmp[j][b],dp[x][j][k]+dp[y][a][b]);}}for(int j=0;j<2;j++)for(int k=0;k<2;k++)dp[x][j][k]=tmp[j][k];}
}
int main(){#ifndef ONLINE_JUDGEfreopen("jiedai.in","r",stdin);
#endifrd(n),rd(m);for(int i=1;i<=m;i++){int a,b;rd(a),rd(b);add_edge(a,b);add_edge(b,a);}memset(dp,0,sizeof(dp));int ans=0;for(int i=1;i<=n;i++)if(!dfn[i]){dfs(i,0);int res=0;for(int j=0;j<2;j++)for(int k=0;k<2;k++)MAX(res,dp[i][j][k]);ans+=res;}printf("%d\n",ans);return 0;
}

DP套DP

给定一个仅由26个大写字母构成的字符串S，统计有多少个仅由大写字母构成的字符串与S的编辑距离恰好为d

|S|≤10，d≤10

const int M=15;
const int P=998244353;
int n,d;
char str[M];
map<vector<int>,int>dp[2];
void add(int &x,int y){x+=y;if(x>=P)x-=P;
}
vector<int>calc(vector<int>lst,char k){vector<int>now(n+1,1e9);for(int i=0;i<=n;i++){if(i)MIN(now[i],lst[i-1]+1);if(i)MIN(now[i],now[i-1]+1);MIN(now[i],lst[i]+1);if(i&&str[i]==k)MIN(now[i],lst[i-1]);}return now;
}
int main(){#ifndef ONLINE_JUDGEfreopen("jiedai.in","r",stdin);
#endifscanf("%s%d",str+1,&d);n=strlen(str+1);vector<int>vi(n+1);for(int i=0;i<=n;i++)vi[i]=i;dp[0][vi]=1;int res=0;for(int i=1;;i++){int cur=i&1;if(dp[!cur].size()==0)break;dp[cur].clear();for(auto it:dp[!cur]){if((it.first)[n]==d)add(res,it.second);for(char k='A';k<='Z';k++){vector<int>vi=calc(it.first,k);int flag=0;for(int j=0;j<=n;j++)if(vi[j]<=d){flag=1;break;}if(flag)add(dp[cur][vi],it.second);}}}printf("%d\n",res);return 0;
}

WQS二分 DP凸优化

给定一个选物品的限制条件，要求恰好选 m 个，让你最大化（最小化）权值，其特点是权值关于选的物品数的函数满足凸包

例题：P2619 给定一个无向带权连通图，每条边是黑色或白色，求恰好有 need 条白边的最小生成树

可以发现，最小生成树的边权和关于白边数量的函数满足下凸包，于是可以套用 WQS 二分

用二分找到一个值，使得所有白边加上这个值后，最小生成树中可以有 need 条白边

const int M=1e5+5;
struct node{int a,b,c,type;bool operator <(const node &A)const{if(c!=A.c)return c<A.c;return type<A.type;}
}E[M];
int n,m,need,fa[M];
int getfa(int x){return x==fa[x]?x:fa[x]=getfa(fa[x]);
}
pair<int,int>solve(int mid){int cnt=0,res=0;for(int i=1;i<=n;i++)fa[i]=i;for(int i=1;i<=m;i++)if(E[i].type==0)E[i].c+=mid;sort(E+1,E+1+m);for(int i=1;i<=m;i++){int x=getfa(E[i].a),y=getfa(E[i].b);if(x!=y){if(E[i].type==0)cnt++;res+=E[i].c;fa[x]=y;}}for(int i=1;i<=m;i++)if(E[i].type==0)E[i].c-=mid;return make_pair(cnt,res);
}
int main(){#ifndef ONLINE_JUDGEfreopen("jiedai.in","r",stdin);
#endifrd(n),rd(m),rd(need);for(int i=1;i<=m;i++)rd(E[i].a),rd(E[i].b),rd(E[i].c),rd(E[i].type),E[i].a++,E[i].b++;int l=-100,r=100,res=0;while(l<=r){int mid=l+r>>1;int cnt=solve(mid).first;if(cnt>=need)l=mid+1,res=mid;else r=mid-1;}printf("%d\n",solve(res).second-res*need);return (0-0);
}

计算几何

基本函数与常量

const double eps=1e-10;
const double PI=acos(-1);
inline int sgn(double x){return fabs(x)<eps?0:(x<0?-1:1);}

二维点类

struct point{double x,y;point(double a=0,double b=0):x(a),y(b){}point operator +(const point &A)const{return point(x+A.x,y+A.y);}point operator -(const point &A)const{return point(x-A.x,y-A.y);}point operator *(const double v)const{return point(x*v,y*v);}point operator /(const double v)const{return point(x/v,y/v);}bool operator ==(const point &A)const{return sgn(x-A.x)==0&&sgn(y-A.y)==0;}double norm(){return sqrt(x*x+y*y);}double norm2(){return x*x+y*y;}
};

点积

double dot(point a,point b){return a.x*b.x+a.y*b.y;
}

叉积

double det(point a,point b){return a.x*b.y-a.y*b.x;
}

两点之间距离

double dist(point a,point b){return (a-b).norm();
}

极角

注意：atan2函数比较慢

double angle(point p){return atan2(p.y,p.x);
}

点绕原点逆时针旋转（弧度）

点p绕原点O逆时针旋转弧度A：rotate_point(p,A)

点a绕点b逆时针旋转弧度c：rotate_point(a-b,c)+b

point rotate_point(point p,double A){return point(p.x*cos(A)-p.y*sin(A),p.x*sin(A)+p.y*cos(A));
}

二维线段直线类

struct line{point a,b;line(){}line(point x,point y):a(x),b(y){}
};

点到直线的距离

double point_to_line(point p,point s,point t){return fabs(det(s-p,t-p)/dist(s,t));
}

点到线段的距离

double point_to_segment(point p,point s,point t){if(sgn(dot(p-s,t-s))==-1)return dist(p,s);if(sgn(dot(p-t,s-t))==-1)return dist(p,t);return fabs(det(s-p,t-p)/dist(s,t));
}

点在直线上的垂足

point point_proj_line(point p,point s,point t){double r=dot(t-s,p-s)/dot(t-s,t-s);return s+(t-s)*r;
}

判断点在直线上

bool point_on_line(point p,point s,point t){return sgn(det(p-s,p-t))==0;
}

判断点在线段上

包含端点，不含端点则将<=改成<

bool point_on_segment(point p,point s,point t){return sgn(det(p-s,p-t))==0&&sgn(dot(p-s,p-t))<=0;
}

判断直线平行

bool parallel(line a,line b){return sgn(det(a.a-a.b,b.a-b.b))==0;
}

两条直线的交点

如果有交点则返回true，且交点保存在res中

需要注意的是，两直线重合也返回false

bool line_make_point(line a,line b,point &res){if(parallel(a,b))return 0;double s1=det(a.a-b.a,b.b-b.a);double s2=det(a.b-b.a,b.b-b.a);res=(a.b*s1-a.a*s2)/(s1-s2);return 1;
}

若两直线不平行，以下代码直接返回交点

point line_make_point(line a,line b){double s1=det(a.a-b.a,b.b-b.a);double s2=det(a.b-b.a,b.b-b.a);return (a.b*s1-a.a*s2)/(s1-s2);
}

三角形的有向面积（需*0.5）

double cross(point o,point a,point b){//*0.5return (a.x-o.x)*(b.y-o.y)-(a.y-o.y)*(b.x-o.x);
}

判断直线和线段是否相交

判断直线a和线段b是否相交

1表示规范相交（只有一个交点且线段端点不在直线上）

2表示不规范相交（线段的一个或两个端点在直线上）

0表示不相交

int line_and_segment(line a,line b){int d1=sgn(cross(a.a,a.b,b.a));int d2=sgn(cross(a.a,a.b,b.b));if((d1^d2)==-2)return 1;if(d1==0||d2==0)return 2;return 0;
}

判断线段和线段是否相交

bool segment_and_segment(line a,line b){if(point_on_line(a.a,b.a,b.b)&&point_on_line(a.b,b.a,b.b))//四点共线return point_on_segment(a.a,b.a,b.b)||point_on_segment(a.b,b.a,b.b)||point_on_segment(b.a,a.a,a.b)||point_on_segment(b.b,a.a,a.b);return sgn(det(a.b-a.a,b.a-a.a)*det(a.b-a.a,b.b-a.a))<=0&&sgn(det(b.b-b.a,a.a-b.a)*det(b.b-b.a,a.b-b.a))<=0;
}

将直线沿法向量方向平移

沿a.a->a.b的左侧平移距离len

line move_d(line a,double len){point e=a.b-a.a;e=rotate_point(e/e.norm(),PI/2);return line(a.a+e*len,a.b+e*len);
}

多边形类

多边形顶点坐标按逆时针排序

struct polygon{static const int SIZE=1005;int n;point A[SIZE];polygon(){}
};

计算多边形面积

 double area(){double res=0;A[n+1]=A[1];for(int i=1;i<=n;i++)res+=det(A[i],A[i+1]);return res/2.0;}

计算多边形周长

 double perimeter(){double res=0;A[n+1]=A[1];for(int i=1;i<=n;i++)res+=(A[i]-A[i+1]).norm();return res;}

判断点是否在多边形内部

 int point_in(point p){int res=0;A[n+1]=A[1];for(int i=1;i<=n;i++){if(point_on_segment(p,A[i],A[i+1]))return 2;int k=sgn(det(A[i+1]-A[i],p-A[i]));int d1=sgn(A[i].y-p.y);int d2=sgn(A[i+1].y-p.y);if(k>0&&d1>0&&d2<=0)res++;if(k<0&&d2>0&&d1<=0)res--;}return res!=0;}

计算多边形的重心

 point mass_center(){point res=point(0,0);if(sgn(area())==0)return res;A[n+1]=A[1];for(int i=1;i<=n;i++)res=res+(A[i]+A[i+1])*det(A[i],A[i+1]);return res/area()/6.0;}

多边形边界上的格点个数

 int gcd(int a,int b){if(b==0)return a;return gcd(b,a%b);}int border_int_point(){int res=0;A[n+1]=A[1];for(int i=1;i<=n;i++)res+=gcd(abs(int(A[i].x-A[i+1].x)),abs(int(A[i].y-A[i+1].y)));return res;}

多边形内的格点个数

Pick定理：对于顶点均为整点的简单多边形，其内部（不含边界）的格点个数为a，其边界上的格点个数为b，其面积为S

它们满足关系式：S=a+b2−1S=a+\frac{b}{2}-1S=a+2b​−1

变形可得：a=S−b2+1a=S-\frac{b}{2}+1a=S−2b​+1

皮克定理与欧拉公式（V−E+F=2V-E+F=2V−E+F=2）等价

 int inside_int_point(){return int(area())-border_int_point()/2+1;}

凸包类

struct polygon_convex{vector<point>p;polygon_convex(int sz=0){p.resize(sz);}
};

计算凸包面积

 double area(){double res=0;for(int i=0;i<p.size();i++)res+=det(p[i],p[(i+1)%p.size()]);return res/2.0;}

计算凸包周长

 double perimeter(){double res=0;for(int i=0;i<p.size();i++)res+=(p[i]-p[(i+1)%p.size()]).norm();return res;}

求点集的凸包（逆时针顺序）

void convex_hull(vector<point>vi,polygon_convex &res){res.p.resize(2*vi.size()+5);sort(vi.begin(),vi.end(),[&](point &a,point &b){if(sgn(a.x-b.x)!=0)return a.x<b.x;return a.y<b.y;});vi.erase(unique(vi.begin(),vi.end()),vi.end());int m=0;for(int i=0;i<vi.size();i++){while(m>1&&sgn(det(res.p[m-1]-res.p[m-2],vi[i]-res.p[m-2]))<=0)m--;res.p[m++]=vi[i];}int k=m;for(int i=int(vi.size())-2;i>=0;i--){while(m>k&&sgn(det(res.p[m-1]-res.p[m-2],vi[i]-res.p[m-2]))<=0)m--;res.p[m++]=vi[i];}res.p.resize(m);if(vi.size()>1)res.p.resize(m-1);
}

判断点是否在凸包内部

需要凸包逆时针

写法一：O(n)O(n)O(n)，true表示点在凸包内部或边界上

bool contain_On(const polygon_convex &a,const point &b){int n=a.p.size(),flag=0;for(int i=0;i<n;i++){int x=sgn(det(a.p[i]-b,a.p[(i+1)%n]-b));if(x){if(flag==0)flag=x;else if(flag!=x)return 0;}}return 1;
}

写法二：O(log2n)O(log_{2}{n})O(log2​n)

半平面交

P4196 凸包求交可以转化为求半平面交

以下代码能求出半平面交所形成的闭合凸包，向量的左侧为半平面。如果闭合凸包不存在，或者半平面无交集，需要另行判断。

const double eps=1e-10;
inline int sgn(double x){return fabs(x)<eps?0:(x<0?-1:1);}
struct point{double x,y;point(double a=0,double b=0):x(a),y(b){}point operator -(const point &A)const{return point(x-A.x,y-A.y);}point operator *(const double v)const{return point(x*v,y*v);}point operator /(const double v)const{return point(x/v,y/v);}
};
struct line{point a,b;double ang;line(){}line(point x,point y):a(x),b(y){}
};
double det(point a,point b){return a.x*b.y-a.y*b.x;
}
point line_make_point(line a,line b){double s1=det(a.a-b.a,b.b-b.a);double s2=det(a.b-b.a,b.b-b.a);return (a.b*s1-a.a*s2)/(s1-s2);
}
struct polygon_convex{vector<point>p;polygon_convex(int sz=0){p.resize(sz);}
};
bool on_right(line l1,line l2,line l3){point o=line_make_point(l2,l3);return sgn(det(l1.b-l1.a,o-l1.a))<0;
}
double angle(line l){return atan2(l.b.y-l.a.y,l.b.x-l.a.x);
}
void halfplane_intersection(vector<line>vi,polygon_convex &res){res.p.clear();for(int i=0;i<vi.size();i++)vi[i].ang=angle(vi[i]);sort(vi.begin(),vi.end(),[&](line l1,line l2){if(sgn(l1.ang-l2.ang)==0)return det(l1.b-l1.a,l2.b-l1.a)<0;return l1.ang<l2.ang;});static const int SIZE=1e6+5;static line Q[SIZE];int L=0,R=0;for(int i=0;i<vi.size();i++){if(i>0&&sgn(angle(vi[i])-angle(vi[i-1]))==0)continue;while(L+1<R&&on_right(vi[i],Q[R-1],Q[R]))R--;while(L+1<R&&on_right(vi[i],Q[L+1],Q[L+2]))L++;Q[++R]=vi[i];}while(L+1<R&&on_right(Q[L+1],Q[R-1],Q[R]))R--;while(L+1<R&&on_right(Q[R],Q[L+1],Q[L+2]))L++;Q[R+1]=Q[L+1];for(int i=L+1;i<=R;i++)res.p.push_back(line_make_point(Q[i],Q[i+1]));
}

平面最远点对（凸包直径）

P1452 平面上有n个点 求出直线距离最近的点对距离

先求出点集的凸包，再在凸包上旋转卡壳

以下代码中，vi为点集，点集的凸包保存在p中，返回值为最远点对距离的平方

double longest_pair(vector<point>vi,vector<point>&p){p.resize(2*vi.size()+5);sort(vi.begin(),vi.end(),[&](point &a,point &b){if(sgn(a.x-b.x)!=0)return a.x<b.x;return a.y<b.y;});vi.erase(unique(vi.begin(),vi.end()),vi.end());int m=0;for(int i=0;i<vi.size();i++){while(m>1&&sgn(det(p[m-1]-p[m-2],vi[i]-p[m-2]))<=0)m--;p[m++]=vi[i];}int k=m;for(int i=int(vi.size())-2;i>=0;i--){while(m>k&&sgn(det(p[m-1]-p[m-2],vi[i]-p[m-2]))<=0)m--;p[m++]=vi[i];}p.resize(m);if(vi.size()>1)p.resize(m-1);int n=p.size();if(n<=1)return 0.0;double res=0;for(int i=0,j=0;i<n;i++){while((j+1)%n!=i&&sgn(det(p[(i+1)%n]-p[i],p[j]-p[i])-det(p[(i+1)%n]-p[i],p[(j+1)%n]-p[i]))<=0)j=(j+1)%n;MAX(res,(p[i]-p[j]).norm2());MAX(res,(p[(i+1)%n]-p[j]).norm2());}return res;
}

三维点类

struct point3{double x,y,z;point3(double a=0,double b=0,double c=0):x(a),y(b),z(c){}point3 operator +(const point3 &A)const{return point3(x+A.x,y+A.y,z+A.z);}point3 operator -(const point3 &A)const{return point3(x-A.x,y-A.y,z-A.z);}point3 operator *(const double v)const{return point3(x*v,y*v,z*v);}point3 operator /(const double v)const{return point3(x/v,y/v,z/v);}bool operator ==(const point3 &A)const{return sgn(x-A.x)==0&&sgn(y-A.y)==0&&sgn(z-A.z)==0;}
};

向量长度

 double length(){return sqrt(x*x+y*y+z*z);}

单位向量

 point3 unit(){return *this/length();}

点积

double dot(point3 a,point3 b){return a.x*b.x+a.y*b.y+a.z*b.z;
}

叉积

point3 det(point3 a,point3 b){return point3(a.y*b.z-a.z*b.y,a.z*b.x-a.x*b.z,a.x*b.y-a.y*b.x);
}

混合积

除以6就是a,b,c三个向量所构成的四面体的体积

double mix(point3 a,point3 b,point3 c){return dot(a,det(b,c));
}

两点之间距离

double dist(point3 a,point3 b){return (a-b).length();
}

三维直线与平面类

struct line3{point3 a,b;line3(){}line3(point3 x,point3 y):a(x),b(y){}double length(){return (a-b).length();}
};

struct plane3{point3 a,b,c;plane3(){}plane3(point3 x,point3 y,point3 z):a(x),b(y),c(z){}
};

平面法向量

point3 normal(point3 a,point3 b,point3 c){return det(a-b,b-c);
}

point3 normal(plane3 s){return det(s.a-s.b,s.b-s.c);
}

判断三点一线

bool points_one_line(point3 a,point3 b,point3 c){return sgn((det(a-b,b-c)).length())==0;
}

判断四点共面

bool points_one_plane(point3 a,point3 b,point3 c,point3 d){return sgn(dot(normal(a,b,c),d-a))==0;
}

判断直线平行

bool parallel(line3 l1,line3 l2){return sgn(det(l1.a-l1.b,l2.a-l2.b).length())==0;
}

判断直线垂直

bool perpendicular(line3 l1,line3 l2){return sgn(dot(l1.a-l1.b,l2.a-l2.b))==0;
}

判断平面内两点在直线同侧

bool same_side(point3 a,point3 b,line3 l){return dot(det(l.a-l.b,a-l.b),det(l.a-l.b,b-l.b))>eps;
}

判断平面内两点在直线异侧

bool opposite_side(point3 a,point3 b,line3 l){return dot(det(l.a-l.b,a-l.b),det(l.a-l.b,b-l.b))<-eps;
}

判断点在线段上（包含端点）

bool point_on_segment_in(point3 p,line3 l){return sgn(det(p-l.a,p-l.b).length())==0&&(l.a.x-p.x)*(l.b.x-p.x)<eps&&(l.a.y-p.y)*(l.b.y-p.y)<eps&&(l.a.z-p.z)*(l.b.z-p.z)<eps;
}

判断点在线段上（不含端点）

bool point_on_segment_ex(point3 p,line3 l){return point_on_segment_in(p,l)&&(sgn(p.x-l.a.x)||sgn(p.y-l.a.y)||sgn(p.z-l.a.z))&&(sgn(p.x-l.b.x)||sgn(p.y-l.b.y)||sgn(p.z-l.b.z));
}

判断线段是否相交（包含端点）

bool segment_and_segment_in(line3 l1,line3 l2){if(points_one_plane(l1.a,l1.b,l2.a,l2.b)==0)return 0;if(points_one_line(l1.a,l2.a,l2.b)&&points_one_line(l1.b,l2.a,l2.b))return point_on_segment_in(l1.a,l2)||point_on_segment_in(l1.b,l2)||point_on_segment_in(l2.a,l1)||point_on_segment_in(l2.b,l1);return same_side(l1.a,l1.b,l2)==0&&same_side(l2.a,l2.b,l1)==0;
}

判断线段是否相交（不含端点）

注意：当线段共线时输出0

bool segment_and_segment_ex(line3 l1,line3 l2){return points_one_plane(l1.a,l1.b,l2.a,l2.b)&&opposite_side(l1.a,l1.b,l2)&&opposite_side(l2.a,l2.b,l1);
}

两条直线的交点（需保证相交）

point3 line_make_point(line3 l1,line3 l2){double t=((l1.a.x-l2.a.x)*(l2.a.y-l2.b.y)-(l1.a.y-l2.a.y)*(l2.a.x-l2.b.x))/((l1.a.x-l1.b.x)*(l2.a.y-l2.b.y)-(l1.a.y-l1.b.y)*(l2.a.x-l2.b.x));return l1.a+(l1.b-l1.a)*t;
}

点到直线的距离

double point_to_line(point3 p,line3 l){return det(p-l.a,l.b-l.a).length()/l.length();
}

直线到直线的距离（不平行）

double line_to_line(line3 l1,line3 l2){point3 n=det(l1.a-l1.b,l2.a-l2.b);return fabs(dot(l1.a-l2.a,n))/n.length();
}

两条直线的夹角的cos值

double angle_cos(line3 l1,line3 l2){return dot(l1.a-l1.b,l2.a-l2.b)/(l1.a-l1.b).length()/(l2.a-l2.b).length();
}

点绕向量逆时针旋转（弧度）

点p绕向量Ol逆时针旋转弧度A：rotate_point(p,l,A)

点p绕直线l逆时针旋转弧度A：rotate_point(p-l.a,l.b-l.a,A)+l.a

point3 rotate_point(point3 p,point3 l,double A){l=l.unit();return p*cos(A)+det(l,p)*sin(A)+l*dot(l,p)*(1-cos(A));
}

判断两点在平面同侧

bool same_side(point3 a,point3 b,plane3 s){return dot(normal(s),a-s.a)*dot(normal(s),b-s.a)>eps;
}

判断两点在平面异侧

bool opposite_side(point3 a,point3 b,plane3 s){return dot(normal(s),a-s.a)*dot(normal(s),b-s.a)<-eps;
}

判断直线与平面平行

bool parallel(line3 l,plane3 s){return sgn(dot(l.a-l.b,normal(s)))==0;
}

判断两平面平行

bool parallel(plane3 s1,plane3 s2){return sgn(det(normal(s1),normal(s2)).length())==0;
}

判断直线与平面垂直

bool perpendicular(line3 l,plane3 s){return sgn(det(l.a-l.b,normal(s)).length())==0;
}

判断两平面垂直

bool perpendicular(plane3 s1,plane3 s2){return sgn(dot(normal(s1),normal(s2)))==0;
}

判断点在三角形内（包含边界）

bool point_in_plane_in(point3 p,plane3 s){return sgn(det(s.a-s.b,s.a-s.c).length()-det(p-s.a,p-s.b).length()-det(p-s.b,p-s.c).length()-det(p-s.c,p-s.a).length())==0;
}

判断点在三角形内（不含边界）

bool point_in_plane_ex(point3 p,plane3 s){return point_in_plane_in(p,s)&&det(p-s.a,p-s.b).length()>eps&&det(p-s.b,p-s.c).length()>eps&&det(p-s.c,p-s.a).length()>eps;
}

判断线段与三角形相交（包含边界）

bool segment_plane_in(line3 l,plane3 s){return same_side(l.a,l.b,s)==0&&same_side(s.a,s.b,plane3(l.a,l.b,s.c))==0&&same_side(s.b,s.c,plane3(l.a,l.b,s.a))==0&&same_side(s.c,s.a,plane3(l.a,l.b,s.b))==0;
}

判断线段与三角形相交（不含边界）

bool segment_plane_ex(line3 l,plane3 s){return opposite_side(l.a,l.b,s)&&opposite_side(s.a,s.b,plane3(l.a,l.b,s.c))&&opposite_side(s.b,s.c,plane3(l.a,l.b,s.a))&&opposite_side(s.c,s.a,plane3(l.a,l.b,s.b));
}

直线与平面的交点

point3 line_plane_make_point(line3 l,plane3 s){point3 n=normal(s);double t=(n.x*(s.a.x-l.a.x)+n.y*(s.a.y-l.a.y)+n.z*(s.a.z-l.a.z))/(n.x*(l.b.x-l.a.x)+n.y*(l.b.y-l.a.y)+n.z*(l.b.z-l.a.z));return l.a+(l.b-l.a)*t;
}

两平面的交线

bool plane_make_line(plane3 s1,plane3 s2,line3 &res){if(parallel(s1,s2))return 0;if(parallel(line3(s2.a,s2.b),s1))res.a=line_plane_make_point(line3(s2.b,s2.c),s1);else res.a=line_plane_make_point(line3(s2.a,s2.b),s1);res.b=res.a+det(normal(s1),normal(s2));return 1;
}

line3 plane_make_line(plane3 s1,plane3 s2){line3 res;if(parallel(line3(s2.a,s2.b),s1))res.a=line_plane_make_point(line3(s2.b,s2.c),s1);else res.a=line_plane_make_point(line3(s2.a,s2.b),s1);res.b=res.a+det(normal(s1),normal(s2));return res;
}

点到面的距离

double point_to_plane(point3 p,plane3 s){return fabs(dot(normal(s),p-s.a))/normal(s).length();
}

两平面的夹角的cos值

需要注意，区分正负

double angle_cos(plane3 s1,plane3 s2){return dot(normal(s1),normal(s2))/normal(s1).length()/normal(s2).length();
}

直线与平面的夹角的sin值

需要注意，区分正负

double angle_sin(line3 l,plane3 s){return dot(l.a-l.b,normal(s))/(l.a-l.b).length()/normal(s).length();
}

圆类

struct circle{point p;double r;circle(){}circle(point p,double r):p(p),r(r){}
};

inline double mysqrt(double x){return sqrt(max(0.0,x));}

圆与直线的交点

inline double mysqrt(double x){return sqrt(max(0.0,x));}
vector<point>circle_line_make_point(point a,point b,point o,double r){double dx=b.x-a.x,dy=b.y-a.y;double A=dx*dx+dy*dy;double B=2*dx*(a.x-o.x)+2*dy*(a.y-o.y);double C=(a.x-o.x)*(a.x-o.x)+(a.y-o.y)*(a.y-o.y)-r*r;double delta=B*B-4*A*C;vector<point>vi;if(sgn(delta)>=0){double t1=(-B+mysqrt(delta))/(2*A);double t2=(-B-mysqrt(delta))/(2*A);vi.push_back(point(a.x+t1*dx,a.y+t1*dy));if(sgn(delta)>0)vi.push_back(point(a.x+t2*dx,a.y+t2*dy));}return vi;
}

最小圆覆盖

平面上n个点，求一个半径最小的圆，能覆盖所有的点

时间复杂度期望为 O(n)O(n)O(n)

struct point{double x,y;point(double a=0,double b=0):x(a),y(b){}point operator +(const point &A)const{return point(x+A.x,y+A.y);}point operator -(const point &A)const{return point(x-A.x,y-A.y);}point operator /(const double v)const{return point(x/v,y/v);}double norm2(){return x*x+y*y;}
};
point circle_center(point a,point b,point c){double a1=b.x-a.x,b1=b.y-a.y,c1=(a1*a1+b1*b1)/2;double a2=c.x-a.x,b2=c.y-a.y,c2=(a2*a2+b2*b2)/2;double d=a1*b2-a2*b1;return point(a.x+(c1*b2-c2*b1)/d,a.y+(a1*c2-a2*c1)/d);
}
void min_circle_cover(int n,point A[],point &o,double &r){srand(time(NULL));random_shuffle(A+1,A+n+1);o=point(0,0),r=0.0;for(int i=1;i<=n;i++)if((A[i]-o).norm2()>r){o=A[i],r=0;for(int j=1;j<i;j++)if((A[j]-o).norm2()>r){o=(A[i]+A[j])/2,r=(A[j]-o).norm2();for(int k=1;k<j;k++)if((A[k]-o).norm2()>r)o=circle_center(A[i],A[j],A[k]),r=(A[k]-o).norm2();}}r=sqrt(r);
}

其他知识或技巧

平面最近点对

P1429 平面上有n个点 求出直线距离最近的点对距离

const double eps=1e-10;
inline int sgn(double x){return fabs(x)<eps?0:(x<0?-1:1);}
struct point{double x,y;point(double a=0,double b=0):x(a),y(b){}point operator -(const point &A)const{return point(x-A.x,y-A.y);}double norm(){return sqrt(x*x+y*y);}
};const int SIZE=2e5+5;
int n;
point p[SIZE];
double res;
void shortest_pair(int l,int r){if(l==r)return;int mid=l+r>>1;shortest_pair(l,mid),shortest_pair(mid+1,r);sort(p+l,p+mid+1,[](point a,point b){return a.y<b.y;});sort(p+mid+1,p+r+1,[](point a,point b){return a.y<b.y;});int now=mid+1;for(int i=l;i<=mid;i++){while(now<=r&&p[i].y-p[now].y>=res)now++;for(int j=now;j<=r&&p[j].y-p[i].y<res;j++)MIN(res,(p[i]-p[j]).norm());}
}int main(){#ifndef ONLINE_JUDGEfreopen("jiedai.in","r",stdin);
#endifcin>>n;for(int i=1;i<=n;i++)cin>>p[i].x>>p[i].y;sort(p+1,p+1+n,[](point a,point b){if(sgn(a.x-b.x)!=0)return a.x<b.x;return a.y<b.y;});res=1e18;shortest_pair(1,n);printf("%.4lf\n",res);return (0-0);
}

三角形重心、外心、垂心、内心

point triangle_mass_center(point a,point b,point c){//重心return (a+b+c)/3;
}
point triangle_circle_center(point a,point b,point c){//外心double a1=b.x-a.x,b1=b.y-a.y,c1=(a1*a1+b1*b1)/2;double a2=c.x-a.x,b2=c.y-a.y,c2=(a2*a2+b2*b2)/2;double d=a1*b2-a2*b1;return point(a.x+(c1*b2-c2*b1)/d,a.y+(a1*c2-a2*c1)/d);
}
point triangle_orthocenter(point a,point b,point c){//垂心return triangle_mass_center(a,b,c)*3-triangle_circle_center(a,b,c)*2;
}
point triangle_innercenter(point a,point b,point c){//内心double la=(b-c).norm();double lb=(c-a).norm();double lc=(a-b).norm();return point((la*a.x+lb*b.x+lc*c.x)/(la+lb+lc),(la*a.y+lb*b.y+lc*c.y)/(la+lb+lc));
}

四面体体积

给定四面体六条棱的长度，求四面体的体积

已知四面体棱 ab,ac,ad,bc,bd,cdab,ac,ad,bc,bd,cdab,ac,ad,bc,bd,cd 的长度分别为 l,n,a,m,b,cl,n,a,m,b,cl,n,a,m,b,c

double tetrahedral_volume(double l,double n,double a,double m,double b,double c){double x,y;x=4*a*a*b*b*c*c-a*a*(b*b+c*c-m*m)*(b*b+c*c-m*m)-b*b*(c*c+a*a-n*n)*(c*c+a*a-n*n);y=c*c*(a*a+b*b-l*l)*(a*a+b*b-l*l)-(a*a+b*b-l*l)*(b*b+c*c-m*m)*(c*c+a*a-n*n);return sqrt(x-y)/12;
}

jiedai算法模板合集（正在肝2021.8.15）相关推荐

1. 最全19000+国外AE模板合集包

   最全19000+国外AE模板合集包 此套模板 包含整理好的 2015-2021年6月 国外AE模板 19000多个 大小:3500+G 模板已编号,含预览视频,方便预览和查找,(少数模板无预览,少数模 ...

2. 研究百度下拉360下拉搜狗下拉神马下拉头条下拉抖音下拉的优化算法（合集帖）

   研究百度下拉360下拉搜狗下拉神马下拉头条下拉抖音下拉的优化算法(合集) 这个标题有点长哈,哈哈哈哈,看的是不是有点懵呢?缩减一点,就是研究百度.360.搜狗.神马.头条.抖音下拉词框的优化算法. 开 ...

3. 【推荐】700套高端简历模板合集

   简历(英语:resume),顾名思义,就是对个人学历.经历.特长.爱好及其它有关情况所作的简明扼要的书面介绍.简历是有针对性的自我介绍的一种规范化.逻辑化的书面表达.对应聘者来说,简历是求职的&quo ...

4. 模板合集TOP50！覆盖8大行业领域，20+业务场景，打包好了直接送

   2022年这么快就进入倒计时啦,工作也陆续进入到冲刺阶段. 打开设计器准备开干的你,有些疑惑是不是像冒泡一样出现: 别着急!大师兄从模板库里,精选出了模板合集TOP50!覆盖8大行业领域,20+业务场 ...

5. 打工人必备，850多份合同范文模板合集免费下载

   打工人必备,850多份合同范文模板合集免费下载  大家好呀,在日常工作中,大家应该很烦找各种合同模板.遇到很多流氓的文库下载还要不少的钱.小编特意整理了常用的各种合同范文模板,包含各行各业,希望能帮助 ...

6. goland 注释模板_有关新学期学习计划模板合集六篇

   有关新学期学习计划模板合集六篇 时间就如同白驹过隙般的流逝,我们又将迎来新一轮的学习,迎来新的喜悦.新的收获,现在这个时候,你会有怎样的计划呢?估计许多人是想得很多,但不会写,下面是小编为大家收集的新 ...

7. 思维导图模板合集及绘制思维导图方法介绍

   思维导图的绘制除了手绘思维导图对框架进行搭建使用之外,还可以选择套用模板进行编辑使用,这两种方法使用起来都很方便并且操作简单,下面是分享的思维导图模板合集以及利用迅捷画图在线网站编辑思维导图操作方法介 ...

8. 课外培训计算机小学课程计划,小学教学计划模板合集四篇

   小学教学计划模板合集四篇 日子如同白驹过隙,不经意间,又将迎来新的工作,新的挑战,现在就让我们好好地规划一下吧.相信大家又在为写计划犯愁了?以下是小编帮大家整理的小学教学计划5篇,供大家参考借鉴,希望 ...

9. 加密货币工具和算法大合集

   2019独角兽企业重金招聘Python工程师标准>>> 维基百科- 加密货币加密货币是一种利用密码学原理来保证账户之间交易的安全性并且控制货币发行量的数字资产. 现在市面上已经存在各 ...

最新文章

1. 说的都是概念——有关编程范式
2. python 代码-20个Python代码段，你需要立刻学会，好用到哭！
3. python银行系统-python银行系统实现源码
4. 第二轮“双一流”名单：这 44 所高校有调整！
5. 屏幕实战效果解析：IPS/TFT/AMOLED/SLCD
6. Linux系统编程8-18总结项目：完成一个简单的自己的shell
7. 三维坐标转经纬度_Three.js 地理坐标和三维空间坐标的转换
8. Net 4.0并行库实用性演练
9. dubbo+rabbitmq+hystrix实现服务的高可用
10. 计算机上相同的用户,php会话为不同计算机上的不同用户发出相同的session_id
11. 【优化求解】基于matlab蚁群算法求解函数极值问题【含Matlab源码 1201期】
12. 广东工业大学通信原理复习笔记第六章数字信号的基带传输（思维导图）含链接和习题
13. Java实现中文汉字转换拼音，解决多音字问题
14. 天涯怎么引流到qq？天涯社区如何靠发帖子引流？
15. DBMS （数据库管理系统） 是什么
16. 浏览器关闭垃圾的百度热点推荐------简洁就是美
17. Gather-Excite:Exploiting Feature Context in Convolutional Neural Networks
18. 自考大专计算机专业英语翻译,大专英语自考(上册)课文翻译及习题答案（138页）-原创力文档...
19. 硬核观察 #612 谷歌正式推出“切换到安卓”应用
20. 武田宣布多项细胞疗法合作，以推进公司的新型免疫肿瘤学阵容

热门文章

1. 学人工智能需要哪些准备？AI基础技能
2. 包装类，正则表示式，Arrays类
3. scala 判断手机号码_手机上的Scala
4. 小i机器人要求苹果公司停售iPhone产品
5. python 使用 numpy 判断两个向量是否平行
6. overcommit_memory=1是否可以改善服务器内存不足的情况
7. Python Basic - GIL：global interpreter lock(全局解释器锁)
8. 快手私信大卡片逆向分析
9. 三、使用Arduino控制WS2812B灯
10. 写一个静态HTML页面，直接写HTML代码和用JS动态生成代码，哪种方式要好