数据结构
树链剖分(bzoj1036)
1 #include<cstring> 2 #include<algorithm> 3 #include<cstdio> 4 #include<vector> 5 using namespace std; 6 const int maxn=30000,inf=1e9; 7 struct wjmzbmr 8 { 9 int tid,top,size,son,dep,father,v; 10 }tree[maxn+50]; 11 struct fotile96 12 { 13 int maxnum,sum; 14 }f[4*maxn+50]; 15 vector<int> g[maxn+50]; 16 int pos[maxn+50]; 17 int n,label=0,q; 18 void dfs(int x,int fa,int dep) 19 { 20 tree[x]=(wjmzbmr){0,0,1,0,dep,fa,tree[x].v}; 21 int m=0; 22 for(int i=0;i<g[x].size();++i) 23 if(g[x][i]!=fa) 24 { 25 dfs(g[x][i],x,dep+1); 26 tree[x].size+=tree[g[x][i]].size; 27 if(tree[g[x][i]].size>m) m=tree[g[x][i]].size,tree[x].son=g[x][i]; 28 } 29 } 30 void connect(int x,int top) 31 { 32 tree[x].tid=++label,pos[tree[x].tid]=x; 33 tree[x].top=top; 34 if(tree[x].son!=0) connect(tree[x].son,top); 35 for(int i=0;i<g[x].size();++i) 36 if(g[x][i]!=tree[x].father&&g[x][i]!=tree[x].son) connect(g[x][i],g[x][i]); 37 } 38 void make(int k,int l,int r) 39 { 40 if(l>r) return; 41 if(l==r) f[k].maxnum=f[k].sum=tree[pos[l]].v; 42 else 43 { 44 int mid=(l+r)>>1; 45 make(k*2,l,mid); 46 make(k*2+1,mid+1,r); 47 f[k].sum=f[k*2].sum+f[k*2+1].sum; 48 f[k].maxnum=max(f[k*2].maxnum,f[k*2+1].maxnum); 49 } 50 } 51 void change(int k,int l,int r,int x,int y) 52 { 53 if(l>r||r<x||l>x) return; 54 if(l==r) 55 { 56 if(l==x) f[k].maxnum=y,f[k].sum=y; 57 return; 58 } 59 int mid=(l+r)>>1; 60 change(k*2,l,mid,x,y); 61 change(k*2+1,mid+1,r,x,y); 62 f[k].sum=f[k*2].sum+f[k*2+1].sum; 63 f[k].maxnum=max(f[k*2].maxnum,f[k*2+1].maxnum); 64 } 65 int qqmax(int k,int l,int r,int ll,int rr) 66 { 67 if(l>r||l>rr||r<ll) return -inf; 68 if(l>=ll&&r<=rr) return f[k].maxnum; 69 int mid=(l+r)>>1; 70 return max(qqmax(k*2,l,mid,ll,rr),qqmax(k*2+1,mid+1,r,ll,rr)); 71 } 72 int qmax(int x,int y) 73 { 74 int ans=-inf; 75 while(tree[x].top!=tree[y].top) 76 { 77 if(tree[tree[x].top].dep<tree[tree[y].top].dep) swap(x,y); 78 ans=max(ans,qqmax(1,1,n,tree[tree[x].top].tid,tree[x].tid)); 79 x=tree[tree[x].top].father; 80 } 81 if(tree[x].dep>tree[y].dep) swap(x,y); 82 ans=max(ans,qqmax(1,1,n,tree[x].tid,tree[y].tid)); 83 return ans; 84 } 85 int qqsum(int k,int l,int r,int ll,int rr) 86 { 87 if(l>r||l>rr||r<ll) return 0; 88 if(l>=ll&&r<=rr) return f[k].sum; 89 int mid=(l+r)>>1; 90 return qqsum(k*2,l,mid,ll,rr)+qqsum(k*2+1,mid+1,r,ll,rr); 91 } 92 int qsum(int x,int y) 93 { 94 int ans=0; 95 while(tree[x].top!=tree[y].top) 96 { 97 if(tree[tree[x].top].dep<tree[tree[y].top].dep) swap(x,y); 98 ans+=qqsum(1,1,n,tree[tree[x].top].tid,tree[x].tid); 99 x=tree[tree[x].top].father; 100 } 101 if(tree[x].dep>tree[y].dep) swap(x,y); 102 ans+=qqsum(1,1,n,tree[x].tid,tree[y].tid); 103 return ans; 104 } 105 int main() 106 { 107 scanf("%d",&n); 108 for(int i=1;i<=n;++i) g[i].clear(); 109 for(int i=1;i<n;++i) 110 { 111 int x,y; 112 scanf("%d %d",&x,&y); 113 g[x].push_back(y); 114 g[y].push_back(x); 115 } 116 for(int i=1;i<=n;++i) scanf("%d",&tree[i].v); 117 dfs(1,0,0); 118 connect(1,1); 119 make(1,1,n); 120 scanf("%d ",&q); 121 for(int i=1;i<=q;++i) 122 { 123 int x,y; 124 char s[10]; 125 scanf("%s %d %d ",s,&x,&y); 126 if(s[0]=='C') change(1,1,n,tree[x].tid,y); 127 if(s[1]=='M') printf("%d ",qmax(x,y)); 128 if(s[1]=='S') printf("%d ",qsum(x,y)); 129 } 130 return 0; 131 }
主席树求区间K小(带修改)
1 //时间是O(nlog^2n) 2 //空间也是O(nlog^2n)，不过可以采用地址回收 3 #include<cstring> 4 #include<algorithm> 5 #include<cstdio> 6 #include<vector> 7 using namespace std; 8 const int maxn=5e4; 9 vector<int> L,R; 10 int sorta[maxn*2+50],a[maxn*2+50],rk[maxn+50]; 11 int ch[maxn*20*20+50][2],sum[maxn*20*20+50],root[maxn+50]; 12 int sz=0,len=0; 13 int n,m,T; 14 int lowbit(int x) 15 { 16 return x&(-x); 17 } 18 int change(int last,int l,int r,int x,int type)//type=1表示加上x，type=-1表示减去x 19 { 20 int k=++len; 21 ch[k][0]=ch[last][0],ch[k][1]=ch[last][1],sum[k]=sum[last]+type; 22 if(l==r) return k; 23 int mid=(l+r)>>1; 24 if(x<=mid) ch[k][0]=change(ch[last][0],l,mid,x,type);else ch[k][1]=change(ch[last][1],mid+1,r,x,type); 25 return k; 26 } 27 int query(int l,int r,int k)//询问[l,R]之间的第k小的值 28 { 29 if(l==r) return l; 30 int suml=0,sumr=0; 31 for(int i=0;i<L.size();++i) suml+=sum[ch[L[i]][0]]; 32 for(int i=0;i<R.size();++i) sumr+=sum[ch[R[i]][0]]; 33 int mid=(l+r)>>1; 34 if(k<=sumr-suml) 35 { 36 for(int i=0;i<L.size();++i) L[i]=ch[L[i]][0]; 37 for(int i=0;i<R.size();++i) R[i]=ch[R[i]][0]; 38 return query(l,mid,k); 39 } 40 else 41 { 42 for(int i=0;i<L.size();++i) L[i]=ch[L[i]][1]; 43 for(int i=0;i<R.size();++i) R[i]=ch[R[i]][1]; 44 return query(mid+1,r,k-(sumr-suml)); 45 } 46 } 47 int main() 48 { 49 scanf("%d",&T); 50 while(T--) 51 { 52 scanf("%d%d",&n,&m); 53 for(int i=1;i<=n;++i) scanf("%d",&a[i]),sorta[i]=a[i]; 54 sort(sorta+1,sorta+n+1); 55 len=0,sz=1; 56 rk[1]=sorta[1]; 57 for(int i=2;i<=n;++i) if(sorta[i]!=sorta[i-1]) sorta[++sz]=sorta[i],rk[sz]=sorta[i];//离散 58 memset(root,0,sizeof(root)); 59 for(int i=1;i<=n;++i) 60 { 61 int id=lower_bound(sorta+1,sorta+sz+1,a[i])-sorta; 62 for(int j=i;j<=n;j+=lowbit(j)) 63 root[j]=change(root[j],1,sz,id,1);//每个点自己建 64 } 65 while(m--) 66 { 67 char c=' '; 68 int x,y,k; 69 while(c!='Q'&&c!='C') c=getchar(); 70 if(c=='C') 71 { 72 scanf("%d%d",&x,&y); 73 int id=lower_bound(sorta+1,sorta+sz+1,a[x])-sorta; 74 for(int i=x;i<=n;i+=lowbit(i)) root[i]=change(root[i],1,sz,id,-1); 75 a[x]=y; 76 id=lower_bound(sorta+1,sorta+sz+1,a[x])-sorta; 77 for(int i=x;i<=n;i+=lowbit(i)) root[i]=change(root[i],1,sz,id,1); 78 } 79 else 80 { 81 scanf("%d%d%d",&x,&y,&k); 82 L.clear(),R.clear(); 83 for(int i=x-1;i;i-=lowbit(i)) L.push_back(root[i]); 84 for(int i=y;i;i-=lowbit(i)) R.push_back(root[i]); 85 printf("%d ",rk[query(1,sz,k)]); 86 } 87 } 88 } 89 return 0; 90 }
splay维护数列问题
1 #include<bits/stdc++.h> 2 using namespace std; 3 const int maxn=2e5,inf=1e9; 4 int pre[maxn+50],ch[maxn+50][2],flip[maxn+50],key[maxn+50],sz[maxn+50],mi[maxn+50],add[maxn+50]; 5 int a[maxn+50]; 6 int n,root=0,len=0,m; 7 struct wjmzbmr 8 { 9 int num,id; 10 bool operator < (const wjmzbmr& x) const 11 { 12 return num<x.num||(num==x.num&&id<x.id); 13 } 14 }b[maxn+50]; 15 void update(int k) 16 { 17 int l=ch[k][0],r=ch[k][1]; 18 sz[k]=1+sz[l]+sz[r]; 19 mi[k]=min(key[k],min(mi[l],mi[r])); 20 } 21 void pushdown(int k) 22 { 23 int& l=ch[k][0]; 24 int& r=ch[k][1]; 25 if(add[k]) 26 { 27 if(l) add[l]+=add[k],mi[l]+=add[k],key[l]+=add[k]; 28 if(r) add[r]+=add[k],mi[r]+=add[k],key[r]+=add[k]; 29 add[k]=0; 30 } 31 if(flip[k]) 32 { 33 flip[k]=0; 34 swap(l,r); 35 if(l) flip[l]^=1; 36 if(r) flip[r]^=1; 37 } 38 } 39 void rorate(int k,int d) 40 { 41 int f=pre[pre[k]],p=pre[k]; 42 pushdown(p),pushdown(k); 43 ch[p][d^1]=ch[k][d]; 44 pre[ch[k][d]]=p; 45 ch[k][d]=p; 46 pre[k]=f; 47 pre[p]=k; 48 if(f) ch[f][ch[f][1]==p]=k;else root=k; 49 update(p),update(k); 50 } 51 void splay(int x,int goal) 52 { 53 pushdown(x); 54 while(pre[x]!=goal) 55 { 56 pushdown(x); 57 if(pre[pre[x]]==goal) rorate(x,ch[pre[x]][0]==x); 58 else 59 { 60 int d1=ch[pre[x]][0]==x,d2=ch[pre[pre[x]]][0]==pre[x]; 61 if(d1==d2) rorate(pre[x],d2),rorate(x,d1); 62 else rorate(x,d1),rorate(x,d2); 63 } 64 } 65 update(x); 66 if(goal==0) root=x; 67 } 68 int build(int l,int r,int fa) 69 { 70 if(l>r) return 0; 71 if(l==r) 72 { 73 pre[l]=fa,ch[l][0]=ch[l][1]=0,add[l]=flip[l]=0,key[l]=mi[l]=a[l],sz[l]=1; 74 return l; 75 } 76 int mid=(l+r)>>1; 77 pre[mid]=fa,ch[mid][0]=ch[mid][1]=0,add[mid]=flip[mid]=0,key[mid]=mi[mid]=a[mid],sz[mid]=1; 78 ch[mid][0]=build(l,mid-1,mid),ch[mid][1]=build(mid+1,r,mid); 79 update(mid); 80 return mid; 81 } 82 void find(int k,int goal) 83 { 84 int t=root; 85 while(true) 86 { 87 pushdown(t); 88 if(sz[ch[t][0]]+1==k) break; 89 if(k<=sz[ch[t][0]]) t=ch[t][0];else k-=sz[ch[t][0]]+1,t=ch[t][1]; 90 } 91 splay(t,goal); 92 } 93 void make(int l,int r,int c)//将区间[l,r]剪掉，接在新序列的第c位后面 94 { 95 find(l-1,0); 96 find(r+1,root); 97 int u=ch[ch[root][1]][0]; 98 find(sz[ch[root][0]]+1+sz[u],ch[root][1]); 99 u=ch[ch[root][1]][0]; 100 ch[ch[root][1]][0]=0; 101 update(ch[root][1]),update(root); 102 find(c,0); 103 ch[u][1]=ch[root][1],ch[root][1]=u; 104 pre[ch[u][1]]=u,pre[u]=root; 105 update(u),update(root); 106 } 107 int main() 108 { 109 scanf("%d",&n); 110 for(int i=2;i<=n+1;++i) scanf("%d",&a[i]); 111 mi[0]=key[0]=inf; 112 root=build(1,n+2,0); 113 len=n+2; 114 scanf("%d",&m); 115 while(m--) 116 { 117 char s[10]; 118 int l,r,x; 119 scanf("%s",s); 120 if(s[0]=='A')//[l,r]加上x 121 { 122 scanf("%d%d%d",&l,&r,&x); 123 ++l,++r; 124 find(l-1,0); 125 find(r+1,root); 126 int u=ch[ch[root][1]][0]; 127 add[u]+=x,key[u]+=x,mi[u]+=x; 128 update(ch[root][1]),update(root); 129 } 130 if(s[0]=='R'&&s[3]=='E')//[l,r]反转 131 { 132 scanf("%d%d",&l,&r); 133 ++l,++r; 134 find(l-1,0); 135 find(r+1,root); 136 int u=ch[ch[root][1]][0]; 137 flip[u]^=1; 138 } 139 if(s[0]=='R'&&s[3]=='O')//[l,r]循环右移x次 140 { 141 scanf("%d%d%d",&l,&r,&x); 142 ++l,++r; 143 x%=r-l+1; 144 make(l,r-x,x+l-1); 145 } 146 if(s[0]=='I')//在l后插入一个位置，初值为x 147 { 148 scanf("%d%d",&l,&x); 149 ++l; 150 find(l,0); 151 int u=ch[root][1]; 152 ++len; 153 pre[len]=root,ch[len][0]=ch[len][1]=0,add[len]=flip[len]=0,key[len]=mi[len]=x,sz[len]=1; 154 ch[root][1]=len; 155 ch[len][1]=u; 156 pre[u]=len; 157 update(len),update(root); 158 } 159 if(s[0]=='D')//删除第i个位置 160 { 161 scanf("%d",&l); 162 ++l; 163 find(l-1,0); 164 find(l+1,root); 165 ch[ch[root][1]][0]=0; 166 update(ch[root][1]),update(root); 167 } 168 if(s[0]=='M')//求[l,r]最小值 169 { 170 scanf("%d%d",&l,&r); 171 ++l,++r; 172 find(l-1,0); 173 find(r+1,root); 174 printf("%d ",mi[ch[ch[root][1]][0]]); 175 } 176 } 177 return 0; 178 }
Link-Cut Tree
1 struct LCT 2 { 3 int ch[maxn+5][2],fa[maxn+5],flip[maxn+5]; 4 int top; 5 int q[maxn+5]; 6 int mx[maxn+5]; 7 int s1[maxn+5],s3[maxn+5],sz[maxn+5];//sz 实际子树大小 s1 指向某个点的非偏爱子树的大小和 s3 splay内子树中每个点的s1的和 8 void init() 9 { 10 for(int i=1;i<=n;++i) mx[i]=i,sz[i]=1,s3[i]=s1[i]=flip[i]=fa[i]=ch[i][0]=ch[i][1]=0; 11 } 12 bool isroot(int x) 13 { 14 return ch[fa[x]][0]!=x&&ch[fa[x]][1]!=x; 15 } 16 void update(int x) 17 { 18 int l=ch[x][0],r=ch[x][1]; 19 mx[x]=x; 20 if(val[mx[l]]>val[mx[x]])mx[x]=mx[l]; 21 if(val[mx[r]]>val[mx[x]])mx[x]=mx[r]; 22 sz[x]=s3[r]+s1[x]+1; 23 s3[x]=s3[l]+s3[r]+s1[x]+1; 24 } 25 void pushdown(int x) 26 { 27 int l=ch[x][0],r=ch[x][1]; 28 if(flip[x]) 29 { 30 flip[x]^=1;flip[l]^=1;flip[r]^=1; 31 swap(ch[x][0],ch[x][1]); 32 update(x); 33 } 34 } 35 void rotate(int &x) 36 { 37 int y=fa[x],z=fa[y],l,r; 38 if(ch[y][0]==x) l=0; 39 else l=1; 40 r=l^1; 41 if(!isroot(y)) 42 { 43 if(ch[z][0]==y) ch[z][0]=x; 44 else ch[z][1]=x; 45 } 46 fa[x]=z;fa[y]=x;fa[ch[x][r]]=y; 47 ch[y][l]=ch[x][r];ch[x][r]=y; 48 update(y);update(x); 49 } 50 void splay(int &x) 51 { 52 top=0; 53 q[++top]=x; 54 for(int i=x;!isroot(i);i=fa[i]) q[++top]=fa[i]; 55 while(top) pushdown(q[top--]); 56 while(!isroot(x)) 57 { 58 int y=fa[x],z=fa[y]; 59 if(!isroot(y)) 60 { 61 if(ch[y][0]==x^ch[z][0]==y) rotate(x); 62 else rotate(y); 63 } 64 rotate(x); 65 } 66 } 67 void access(int x) 68 { 69 for(int t=0;x;t=x,x=fa[x]) 70 { 71 splay(x); 72 s1[x]+=s3[ch[x][1]]; 73 ch[x][1]=t; 74 s1[x]-=s3[t]; 75 update(x); 76 } 77 } 78 void makeroot(int x) 79 { 80 access(x);splay(x);flip[x]^=1; 81 pushdown(x); 82 } 83 bool linked(int x,int y) 84 { 85 //判断点x和点y是否在一个连通块中 86 makeroot(x);access(y);splay(y); 87 return x==y||fa[x]; 88 } 89 void link(int x,int y) 90 { 91 makeroot(x);fa[x]=y; 92 access(y),splay(y),s1[y]+=s3[x]; 93 update(y); 94 } 95 void cut(int x,int y) 96 { 97 makeroot(x);access(y);splay(y); 98 ch[y][0]=fa[x]=0; 99 sz[y]-=s3[x]; 100 s3[y]-=s3[x]; 101 update(y); 102 } 103 int findfather(int root,int x) 104 { 105 //返回实际树上以root为根，x的父亲 106 makeroot(root); 107 access(x),splay(x); 108 pushdown(x); 109 int fa=ch[x][0]; 110 while(true) 111 { 112 pushdown(fa); 113 if(!ch[fa][1]) break; 114 fa=ch[fa][1]; 115 } 116 return fa; 117 118 } 119 int querymax(int x,int y) 120 { 121 makeroot(x);access(y);splay(y); 122 return mx[y]; 123 } 124 }lct;
Kdtree
1 #include<bits/stdc++.h> 2 using namespace std; 3 const int maxn=500000,inf=1e9; 4 int n,m,cur,ans,root; 5 int x[maxn+50],y[maxn+50]; 6 struct P 7 { 8 int mn[2],mx[2],d[2],lch,rch; 9 int& operator[](int x) {return d[x];} 10 friend bool operator<(P x,P y) {return x[cur]<y[cur];} 11 friend int dis(P x,P y) {return abs(x[0]-y[0])+abs(x[1]-y[1]);} 12 }p[maxn+50]; 13 struct kdtree 14 { 15 P t[2*maxn+50],T; 16 int ans; 17 void update(int k) 18 { 19 int l=t[k].lch,r=t[k].rch; 20 for (int i=0;i<2;i++) 21 { 22 t[k].mn[i]=t[k].mx[i]=t[k][i]; 23 if(l) t[k].mn[i]=min(t[k].mn[i],t[l].mn[i]); 24 if(r) t[k].mn[i]=min(t[k].mn[i],t[r].mn[i]); 25 if(l) t[k].mx[i]=max(t[k].mx[i],t[l].mx[i]); 26 if(r) t[k].mx[i]=max(t[k].mx[i],t[r].mx[i]); 27 } 28 } 29 int build(int l,int r,int now) 30 { 31 cur=now; 32 int mid=(l+r)/2; 33 nth_element(p+l,p+mid,p+r+1); 34 t[mid]=p[mid]; 35 for(int i=0;i<2;i++) t[mid].mx[i]=t[mid].mn[i]=t[mid][i]; 36 if(l<mid) t[mid].lch=build(l,mid-1,now^1); 37 if(r>mid) t[mid].rch=build(mid+1,r,now^1); 38 update(mid); 39 return mid; 40 } 41 void insert(int k,int now) 42 { 43 if(T[now]<t[k][now]) 44 { 45 if(t[k].lch) insert(t[k].lch,now^1); 46 else 47 { 48 t[k].lch=++n; 49 t[n]=T; 50 for(int i=0;i<2;++i) t[n].mx[i]=t[n].mn[i]=t[n][i]; 51 } 52 } 53 else 54 { 55 if(t[k].rch) insert(t[k].rch,now^1); 56 else 57 { 58 t[k].rch=++n; 59 t[n]=T; 60 for(int i=0;i<2;++i) t[n].mx[i]=t[n].mn[i]=t[n][i]; 61 } 62 } 63 update(k); 64 } 65 void ins(int x,int y) 66 { 67 T[0]=x,T[1]=y; 68 T.lch=T.rch=0; 69 insert(root,0); 70 } 71 int getmn(P x) 72 { 73 int ans=0; 74 for(int i=0;i<2;i++) 75 { 76 ans+=max(T[i]-x.mx[i],0); 77 ans+=max(x.mn[i]-T[i],0); 78 } 79 return ans; 80 }//估价函数，注意如果是欧几里得距离辣么估价函数要修改 81 int getmx(P x) 82 { 83 int ans=0; 84 for(int i=0;i<2;i++) ans+=max(abs(T[i]-x.mn[i]),abs(T[i]-x.mx[i])); 85 return ans; 86 } 87 void querymx(int k) 88 { 89 ans=max(ans,dis(t[k],T)); 90 int l=t[k].lch,r=t[k].rch,dl=-inf,dr=-inf; 91 if (l) dl=getmx(t[l]); 92 if (r) dr=getmx(t[r]); 93 if (dl>dr) 94 { 95 if (dl>ans) querymx(l); 96 if (dr>ans) querymx(r); 97 } 98 else 99 { 100 if (dr>ans) querymx(r); 101 if (dl>ans) querymx(l); 102 } 103 } 104 void querymn(int k) 105 { 106 //if(dis(t[k],T)) ans=min(ans,dis(t[k],T)); 107 ans=min(ans,dis(t[k],T)); 108 int l=t[k].lch,r=t[k].rch,dl=inf,dr=inf; 109 if(l) dl=getmn(t[l]); 110 if(r) dr=getmn(t[r]); 111 if (dl<dr) 112 { 113 if (dl<ans) querymn(l); 114 if (dr<ans) querymn(r); 115 } 116 else 117 { 118 if (dr<ans) querymn(r); 119 if (dl<ans) querymn(l); 120 } 121 } 122 int query(int f,int x,int y) 123 { 124 T[0]=x,T[1]=y; 125 T.lch=T.rch=0; 126 if (f==0) ans=-inf,querymx(root); 127 else 128 ans=inf,querymn(root); 129 return ans; 130 } 131 }kdtree; 132 int main() 133 { 134 135 scanf("%d %d",&n,&m); 136 for(int i=1;i<=n;++i) 137 { 138 scanf("%d%d",&x[i],&y[i]); 139 p[i][0]=x[i],p[i][1]=y[i]; 140 } 141 root=kdtree.build(1,n,0); 142 for(int i=1;i<=m;++i) 143 { 144 int t,x,y; 145 scanf("%d %d %d",&t,&x,&y); 146 if(t==1) kdtree.ins(x,y); 147 if(t==2) printf("%d ",kdtree.query(1,x,y)); 148 } 149 return 0; 150 }
Kdtree套替罪羊
1 #include<bits/stdc++.h> 2 using namespace std; 3 const int maxn=500000,inf=1e9; 4 int n,m,cur,ans,root; 5 int x[maxn+50],y[maxn+50]; 6 struct P 7 { 8 int mn[2],mx[2],d[2],lch,rch,sz; 9 int& operator[](int x) {return d[x];} 10 friend bool operator<(P x,P y) {return x[cur]<y[cur];} 11 friend int dis(P x,P y) {return abs(x[0]-y[0])+abs(x[1]-y[1]);} 12 }p[maxn*2+50]; 13 namespace kdtree 14 { 15 const double alpha=0.75; 16 P t[2*maxn+50],T; 17 int id[maxn+5]; 18 int ans,len,tot; 19 bool cmp(const int x,const int y) 20 { 21 return p[x][cur]<p[y][cur]; 22 } 23 void update(int k) 24 { 25 int l=t[k].lch,r=t[k].rch; 26 t[k].sz=t[t[k].lch].sz+t[t[k].rch].sz+1; 27 for (int i=0;i<2;i++) 28 { 29 t[k].mn[i]=t[k].mx[i]=t[k][i]; 30 if(l) t[k].mn[i]=min(t[k].mn[i],t[l].mn[i]); 31 if(r) t[k].mn[i]=min(t[k].mn[i],t[r].mn[i]); 32 if(l) t[k].mx[i]=max(t[k].mx[i],t[l].mx[i]); 33 if(r) t[k].mx[i]=max(t[k].mx[i],t[r].mx[i]); 34 } 35 } 36 int build(int l,int r,int now) 37 { 38 if(l>r) return 0; 39 cur=now; 40 int mid=(l+r)/2; 41 nth_element(id+l,id+mid,id+r+1,cmp); 42 int k=id[mid]; 43 t[k]=p[k]; 44 t[k].sz=1,t[k].lch=t[k].rch=0; 45 for(int i=0;i<2;i++) t[k].mx[i]=t[k].mn[i]=t[k][i]; 46 if(l<mid) t[k].lch=build(l,mid-1,now^1); 47 if(r>mid) t[k].rch=build(mid+1,r,now^1); 48 update(k); 49 return k; 50 } 51 void getid(int k) 52 { 53 id[++tot]=k; 54 p[k]=t[k]; 55 if(t[k].lch) getid(t[k].lch); 56 if(t[k].rch) getid(t[k].rch); 57 } 58 int rebuild(int k,int now) 59 { 60 tot=0; 61 getid(k); 62 return build(1,tot,now); 63 } 64 int insert(int k,int now) 65 { 66 if(T[now]<t[k][now]) 67 { 68 if(t[k].lch) t[k].lch=insert(t[k].lch,now^1); 69 else 70 { 71 t[k].lch=++len; 72 t[len]=T; 73 for(int i=0;i<2;++i) t[len].mx[i]=t[len].mn[i]=t[len][i]; 74 } 75 } 76 else 77 { 78 if(t[k].rch) t[k].rch=insert(t[k].rch,now^1); 79 else 80 { 81 t[k].rch=++len; 82 t[len]=T; 83 for(int i=0;i<2;++i) t[len].mx[i]=t[len].mn[i]=t[len][i]; 84 } 85 } 86 update(k); 87 if(t[k].sz*alpha+3<max(t[t[k].lch].sz,t[t[k].rch].sz)) k=rebuild(k,now); 88 return k; 89 } 90 void ins(int x,int y) 91 { 92 T[0]=x,T[1]=y; 93 T.sz=1; 94 T.lch=T.rch=0; 95 if(root==0) 96 { 97 len=1; 98 t[len]=T; 99 update(len); 100 root=len; 101 return; 102 } 103 root=insert(root,0); 104 } 105 int getmn(P x) 106 { 107 int ans=0; 108 for(int i=0;i<2;i++) 109 { 110 ans+=max(T[i]-x.mx[i],0); 111 ans+=max(x.mn[i]-T[i],0); 112 } 113 return ans; 114 }//估价函数，注意如果是欧几里得距离辣么估价函数要修改 115 int getmx(P x) 116 { 117 int ans=0; 118 for(int i=0;i<2;i++) ans+=max(abs(T[i]-x.mn[i]),abs(T[i]-x.mx[i])); 119 return ans; 120 } 121 void querymx(int k) 122 { 123 ans=max(ans,dis(t[k],T)); 124 int l=t[k].lch,r=t[k].rch,dl=-inf,dr=-inf; 125 if (l) dl=getmx(t[l]); 126 if (r) dr=getmx(t[r]); 127 if (dl>dr) 128 { 129 if (dl>ans) querymx(l); 130 if (dr>ans) querymx(r); 131 } 132 else 133 { 134 if (dr>ans) querymx(r); 135 if (dl>ans) querymx(l); 136 } 137 } 138 void querymn(int k) 139 { 140 //if(dis(t[k],T)) ans=min(ans,dis(t[k],T)); 141 ans=min(ans,dis(t[k],T)); 142 int l=t[k].lch,r=t[k].rch,dl=inf,dr=inf; 143 if(l) dl=getmn(t[l]); 144 if(r) dr=getmn(t[r]); 145 if (dl<dr) 146 { 147 if (dl<ans) querymn(l); 148 if (dr<ans) querymn(r); 149 } 150 else 151 { 152 if (dr<ans) querymn(r); 153 if (dl<ans) querymn(l); 154 } 155 } 156 int query(int f,int x,int y) 157 { 158 T[0]=x,T[1]=y; 159 T.lch=T.rch=0; 160 if (f==0) ans=-inf,querymx(root); 161 else 162 ans=inf,querymn(root); 163 return ans; 164 } 165 } 166 int main() 167 { 168 read(n); 169 read(m); 170 for(int i=1;i<=n;++i) 171 { 172 read(x[i]); 173 read(y[i]); 174 p[i][0]=x[i],p[i][1]=y[i]; 175 } 176 kdtree::len=n; 177 for(int i=1;i<=n;++i) kdtree::id[i]=i; 178 root=kdtree::build(1,n,0); 179 for(int i=1;i<=m;++i) 180 { 181 int t,x,y; 182 read(t); 183 read(x); 184 read(y); 185 if(t==1) kdtree::ins(x,y); 186 if(t==2) printf("%d ",kdtree::query(1,x,y)); 187 } 188 return 0; 189 }
左偏树
1 //APIO2012 dispatching 2 #include<bits/stdc++.h> 3 using namespace std; 4 const int maxn=1e5; 5 vector<int> g[maxn+50]; 6 long long cost[maxn+50],lead[maxn+50],sum[maxn+50],num[maxn+50]; 7 int l[maxn+50],r[maxn+50],dist[maxn+50]; 8 int n,m,root; 9 long long ans=0; 10 int merge(int u,int v)//左偏树核心操作——merge：合并两个左偏树 11 { 12 if(!u) return v; 13 if(!v) return u; 14 if(cost[u]<cost[v]) swap(u,v);//此处是大根堆 15 r[u]=merge(r[u],v); 16 if(dist[r[u]]>dist[l[u]]) swap(l[u],r[u]);//时刻保证dist(l)>=dist(r) 17 if(r[u]) dist[u]=dist[r[u]]+1;else dist[u]=0;//更新dist数组 18 num[u]=num[l[u]]+num[r[u]]+1; 19 sum[u]=sum[l[u]]+sum[r[u]]+cost[u];//维护节点信息 20 return u; 21 } 22 int del(int u) 23 { 24 return merge(l[u],r[u]);//删除操作就是去掉根节点,merge左右儿子 25 } 26 int dfs(int k) 27 { 28 int u=k; 29 num[u]=1,sum[u]=cost[u]; 30 while(sum[u]>m) u=del(u); 31 for(int i=0;i<g[k].size();++i) 32 { 33 int v=dfs(g[k][i]); 34 u=merge(u,v); 35 while(sum[u]>m) u=del(u); 36 } 37 ans=max(ans,lead[k]*num[u]); 38 return u; 39 } 40 int main() 41 { 42 scanf("%d%d",&n,&m); 43 for(int i=1;i<=n;++i) g[i].clear(); 44 for(int i=1;i<=n;++i) 45 { 46 int x; 47 scanf("%d%d%d",&x,&cost[i],&lead[i]); 48 if(lead[i]==0) root=i; 49 g[x].push_back(i); 50 } 51 memset(l,0,sizeof(l)); 52 memset(r,0,sizeof(r)); 53 memset(dist,0,sizeof(dist)); 54 memset(sum,0,sizeof(sum)); 55 memset(num,0,sizeof(num)); 56 dfs(1); 57 printf("%lld",ans); 58 return 0; 59 }
虚树
1 void buildtree() 2 { 3 len=0; 4 sort(a+1,a+m+1,cmp);//关键点 5 id.clear();//记录出现的关键点 6 id.push_back(a[1]); 7 s[++len]=a[1]; 8 for(int i=2;i<=m;i++) 9 { 10 int x=a[i],f=lca(s[len],x); 11 while(deep[f]<deep[s[len]]) 12 { 13 if(deep[f]>=deep[s[len-1]]) 14 { 15 addedge(f,s[len],abs(deep[f]-deep[s[len]]));//虚树addedge(u,v,w) 16 len--; 17 if(f!=s[len]) s[++len]=f; 18 break; 19 } 20 addedge(s[len-1],s[len],abs(deep[s[len-1]]-deep[s[len]])); 21 len--; 22 } 23 if(s[len]!= x) s[++len]=x; 24 } 25 while(--len) addedge(s[len],s[len+1],abs(deep[s[len+1]]-deep[s[len]])); 26 }
Split-Merge Treap
1 #include<bits/stdc++.h> 2 using namespace std; 3 const int maxn=1e5,inf=1e9,mod=1e6; 4 int pre[maxn+5],ch[maxn+5][2],val[maxn+5],key[maxn+5]; 5 int n,root,len; 6 int Rand() 7 { 8 int ans=0; 9 for(int i=1;i<=8;++i) ans=ans*10+rand()%10; 10 return ans; 11 } 12 void split(int k,int x,int &u,int &v) 13 { 14 if(!k) u=v=0; 15 else 16 { 17 pushdown(k); 18 if(x<=sz[ch[k][0]]) v=k,split(ch[k][0],x,u,ch[k][0]); 19 else u=k,split(ch[k][1],x-sz[ch[k][0]]-1,ch[k][1],v); 20 /* if(val[k]<=x) u=k,split(ch[k][1],x,ch[k][1],v); 21 else v=k,split(ch[k][0],x,u,ch[k][0]); */ 22 update(k); 23 } 24 } 25 int merge(int u,int v) 26 { 27 /* 28 将两个子树u,v合并，其中u子树的权值都小于v子树 29 */ 30 if(!u) return v; 31 if(!v) return u; 32 pushdown(u),pushdown(v); 33 if(key[u]<key[v]) 34 { 35 ch[u][1]=merge(ch[u][1],v); 36 update(u); 37 return u; 38 } 39 else 40 { 41 ch[v][0]=merge(u,ch[v][0]); 42 update(v); 43 return v; 44 } 45 } 46 int newnode(int x) 47 { 48 /* 49 新建一个权值为x的结点 50 */ 51 ++len; 52 val[len]=x; 53 ch[len][0]=ch[len][1]=0; 54 key[len]=Rand()%1000000000; 55 return len; 56 } 57 void insert(int x) 58 { 59 /* 60 插入一个权值为x的结点 61 */ 62 int u,v; 63 split(root,x,u,v); 64 root=merge(merge(u,newnode(x)),v); 65 } 66 void del(int x) 67 { 68 /* 69 删除一个权值为x的结点 70 */ 71 int u,v,c,d; 72 split(root,x,u,v); 73 split(u,x-1,c,d); 74 d=merge(ch[d][0],ch[d][1]); 75 u=merge(c,d); 76 root=merge(u,v); 77 } 78 int getpre(int k,int x) 79 { 80 /* 81 从以k为根的子树中找到x的前驱 82 */ 83 if(!k) return -inf; 84 if(x==val[k]) return x; 85 if(x<val[k]) return getpre(ch[k][0],x); 86 return max(val[k],getpre(ch[k][1],x)); 87 } 88 int getsuf(int k,int x) 89 { 90 /* 91 从以k为根的子树中找到x的后继 92 */ 93 if(!k) return inf; 94 if(x==val[k]) return x; 95 if(x>val[k]) return getsuf(ch[k][1],x); 96 return min(val[k],getsuf(ch[k][0],x)); 97 } 98 int main() 99 { 100 srand(time(0)); 101 scanf("%d",&n); 102 root=0; 103 len=0; 104 insert(-inf); 105 insert(inf); 106 int ans=0; 107 int a=0,b=0; 108 for(int i=1;i<=n;++i) 109 { 110 int op,x; 111 scanf("%d%d",&op,&x); 112 if(op==0) 113 { 114 if(b) 115 { 116 int ans1=getpre(root,x); 117 int ans2=getsuf(root,x); 118 if(ans2-x<x-ans1) ans=(ans+abs(ans2-x))%mod,del(ans2);else ans=(ans+abs(ans1-x))%mod,del(ans1); 119 --b; 120 } 121 else ++a,insert(x); 122 } 123 else 124 { 125 if(a) 126 { 127 int ans1=getpre(root,x); 128 int ans2=getsuf(root,x); 129 if(ans2-x<x-ans1) ans=(ans+abs(ans2-x))%mod,del(ans2);else ans=(ans+abs(ans1-x))%mod,del(ans1); 130 --a; 131 } 132 else ++b,insert(x); 133 } 134 } 135 printf("%d ",ans); 136 return 0; 137 }
点分治
1 int getroot(int u) 2 { 3 int t=1; 4 q[0]=u; 5 f[u]=-1; 6 for(int i=0;i<t;++i) 7 { 8 u=q[i]; 9 for(auto v:g[u]) 10 if (!vis[v]&&v!=f[u]) f[q[t++]=v]=u; 11 mx[q[i]]=0; 12 sz[q[i]]=1; 13 } 14 for (int i=t-1;i>=0;i--) 15 { 16 mx[q[i]]=max(mx[q[i]],t-sz[q[i]]); 17 if(mx[q[i]]*2<=t) return q[i]; 18 sz[f[q[i]]]+=sz[q[i]]; 19 mx[f[q[i]]]=max(mx[f[q[i]]],sz[q[i]]); 20 } 21 return 0; 22 } 23 void getdeep(int k,int fa,int f) 24 { 25 father[k]=fa; 26 dep[k]=dep[fa]+1; 27 maxd[k]=dep[k]; 28 par[k]=f; 29 L[k]=++t; 30 p[t]=k; 31 for(auto u:g[k]) 32 { 33 if(u==fa||vis[u]) continue; 34 getdeep(u,k,f); 35 maxd[k]=max(maxd[k],maxd[u]); 36 } 37 R[k]=t; 38 } 39 void solve(int k) 40 { 41 t=0; 42 maxd[k]=dep[k]=0; 43 L[k]=++t; 44 p[1]=k; 45 for(auto u:g[k]) 46 { 47 if(vis[u]) continue; 48 getdeep(u,k,u); 49 maxd[k]=max(maxd[k],maxd[u]); 50 } 51 R[k]=t; 52 53 for(int i=0;i<=maxd[k];++i) s[i]=0; 54 for(int i=1;i<=t;++i) 55 { 56 int u=p[i]; 57 if(u>n) continue; 58 ++s[dep[u]]; 59 if(dep[u]==D) res[k][par[u]]++; 60 } 61 ss[0]=s[0]; 62 for(int i=1;i<=maxd[k];++i) ss[i]=ss[i-1]+s[i]; 63 ans[k]+=ss[min(D-1,maxd[k])]; 64 for(int i=2;i<=t;++i) 65 { 66 int u=p[i]; 67 if(D-1-dep[u]>=0) ans[u]+=ss[min(D-1-dep[u],maxd[k])]; 68 if(D>=dep[u]&&D-dep[u]<=maxd[k]) 69 res[u][father[u]]+=s[D-dep[u]]; 70 } 71 72 for(auto u:g[k]) 73 { 74 if(vis[u]) continue; 75 for(int i=0;i<=maxd[u];++i) s[i]=ss[i]=0; 76 for(int i=L[u];i<=R[u];++i) 77 { 78 int v=p[i]; 79 if(v>n) continue; 80 ++s[dep[v]]; 81 } 82 ss[0]=s[0]; 83 for(int i=1;i<=maxd[u];++i) ss[i]=ss[i-1]+s[i]; 84 for(int i=L[u];i<=R[u];++i) 85 { 86 int v=p[i]; 87 if(D-1-dep[v]>=0) ans[v]-=ss[min(D-1-dep[v],maxd[u])]; 88 if(D>=dep[v]&&D-dep[v]<=maxd[u]) 89 res[v][father[v]]-=s[D-dep[v]]; 90 } 91 } 92 vis[k]=1; 93 for(auto u:g[k]) 94 { 95 if(vis[u]) continue; 96 sum=sz[u]; 97 root=getroot(u); 98 solve(root); 99 } 100 }
莫队算法
1 #include<bits/stdc++.h> 2 using namespace std; 3 const int maxn=40000; 4 int pos[maxn+5],a[maxn+5]; 5 struct wjmzbmr 6 { 7 int l,r,k1,k2,id,ans; 8 }q[maxn+5]; 9 int block,n,m; 10 int color[maxn+5],cnt[200+5],ans[maxn+5]; 11 int num[maxn+5],sum[maxn+5][200+5]; 12 bool cmp(const wjmzbmr a,const wjmzbmr b) 13 { 14 if(a.l/block!=b.l/block) return a.l<b.l; 15 if((a.l/block)&1) return a.r>b.r;else return a.r<b.r; 16 } 17 void update(int x,int type) 18 { 19 --color[num[x]]; 20 if(color[num[x]]==0) --cnt[num[x]/block]; 21 --sum[num[x]][x/block]; 22 num[x]+=type; 23 ++sum[num[x]][x/block]; 24 if(color[num[x]]==0) ++cnt[num[x]/block]; 25 ++color[num[x]]; 26 } 27 int query(int k1,int k2) 28 { 29 int i=0; 30 for(i=0;i<=n/block;++i) 31 if(cnt[i]<k1) k1-=cnt[i]; 32 else break; 33 int turn; 34 for(turn=block*i;turn<block*(i+1);++turn) 35 { 36 if(color[turn]>0) --k1; 37 if(k1==0) break; 38 } 39 for(i=0;i<=200;++i) 40 if(sum[turn][i]<k2) k2-=sum[turn][i]; 41 else break; 42 int ans; 43 for(ans=block*i;ans<block*(i+1);++ans) 44 { 45 if(num[ans]!=turn) continue; 46 if(num[ans]>0) --k2; 47 if(k2==0) break; 48 } 49 return ans; 50 } 51 int main() 52 { 53 scanf("%d",&n); 54 for(int i=1;i<=n;++i) scanf("%d",&a[i]); 55 scanf("%d",&m); 56 for(int i=1;i<=m;++i) scanf("%d%d%d%d",&q[i].l,&q[i].r,&q[i].k1,&q[i].k2),q[i].id=i; 57 block=sqrt(n); 58 sort(q+1,q+m+1,cmp); 59 block=200; 60 int l=1,r=0; 61 for(int i=1;i<=m;++i) 62 { 63 while(l>q[i].l) update(a[--l],1); 64 while(r<q[i].r) update(a[++r],1); 65 while(r>q[i].r) update(a[r--],-1); 66 while(l<q[i].l) update(a[l++],-1); 67 ans[q[i].id]=query(q[i].k1,q[i].k2); 68 } 69 for(int i=1;i<=m;++i) printf("%d ",ans[i]); 70 return 0; 71 }
带修改莫队
1 #include<bits/stdc++.h> 2 using namespace std; 3 const int maxn=1e6; 4 int block,n,m,num,s,l,r; 5 struct Q 6 { 7 int l,r,t,id; 8 bool operator < (const Q& x) const 9 { 10 if(l/block!=x.l/block) return l<x.l; 11 if(r/block!=x.r/block) return r<x.r; 12 return t<x.t; 13 } 14 }q[maxn+5]; 15 struct wjmzbmr 16 { 17 int pos,from,to,t; 18 }c[maxn+5]; 19 int a[maxn+5],b[maxn+5],color[maxn+5],ans[maxn+5]; 20 void add(int x,int type) 21 { 22 if(color[x]) --s; 23 color[x]+=type; 24 if(color[x]) ++s; 25 } 26 void ins(wjmzbmr q) 27 { 28 if(q.pos>=l&&q.pos<=r) add(q.from,-1),add(q.to,1); 29 a[q.pos]=q.to; 30 } 31 void del(wjmzbmr q) 32 { 33 if(q.pos>=l&&q.pos<=r) add(q.to,-1),add(q.from,1); 34 a[q.pos]=q.from; 35 } 36 int main() 37 { 38 int T; 39 scanf("%d%d",&n,&T); 40 for(int i=1;i<=n;++i) scanf("%d",&b[i]),a[i]=b[i]; 41 int time=0; 42 while(T--) 43 { 44 char s[3]; 45 int x,y; 46 scanf("%s%d%d",s,&x,&y); 47 if(s[0]=='Q') ++m,q[m]={x,y,++time,m}; 48 else c[++num]={x,0,y,++time}; 49 } 50 for(int i=1;i<=num;++i) c[i].from=b[c[i].pos],b[c[i].pos]=c[i].to; 51 c[num+1].t=10000000; 52 block=(int)(ceil(pow(n,2.0/3))); 53 sort(q+1,q+m+1); 54 l=1,r=0; 55 int now=0; 56 s=0; 57 for(int i=1;i<=m;++i) 58 { 59 while(r<q[i].r) add(a[++r],1); 60 while(l>q[i].l) add(a[--l],1); 61 while(r>q[i].r) add(a[r--],-1); 62 while(l<q[i].l) add(a[l++],-1); 63 while(c[now+1].t<q[i].t) ins(c[++now]); 64 while(c[now].t>q[i].t) del(c[now--]); 65 ans[q[i].id]=s; 66 } 67 for(int i=1;i<=m;++i) printf("%d ",ans[i]); 68 return 0; 69 }
可持久化堆优化K短路
1 /* 2 可持久化堆优化k短路 3 求最短路径树：O(nlogn) 4 求k短路：O(mlogm+klogk) 5 空间复杂度：O(mlogm)，但具体开数组要注意常数，要看执行merge操作的次数 6 */ 7 #include<bits/stdc++.h> 8 using namespace std; 9 const int maxn=100000,maxm=200000,maxsize=3000000;//maxsize是堆的最大个数 10 const long long inf=1000000000000000LL; 11 int a[maxn+5]; 12 long long s[maxn+5]; 13 int n,k,m,len,tot; 14 int S,T; 15 vector<int> g1[maxn+5]; 16 int mi[maxn+5],mx[maxn+5]; 17 int head[maxn+5],nx[maxm+5]; 18 long long dis[maxn+5]; 19 int pos[maxn+5],rt[maxn+5]; 20 struct Edge 21 { 22 int to; 23 long long w; 24 }e[maxm+5]; 25 struct node 26 { 27 /* 28 u是当前堆顶的节点编号 29 key是当前堆顶对应边的权值,边的权值定义为：走这条边要多绕多少路 30 l,r分别是堆左右儿子的地址 31 */ 32 int u; 33 long long key; 34 int l,r; 35 }H[maxsize+5]; 36 struct heapnode 37 { 38 /* 39 求k短路时候用到的数据结构 40 len表示1~倒数第二条边的边权和 41 root表示倒数第二条边的tail的H[tail]，其中堆顶就是最后一条边 42 */ 43 long long len; 44 int root; 45 bool operator < (const heapnode& x) const 46 { 47 return len+H[root].key>x.len+H[x.root].key; 48 } 49 }; 50 priority_queue<heapnode> q; 51 void addedge(int u,int v,long long w) 52 { 53 54 //printf("%d %d %lld ",u,v,w); 55 e[++len]={v,w}; 56 nx[len]=head[u]; 57 head[u]=len; 58 } 59 void dfs(int k,int fa) 60 { 61 s[k]=a[k]; 62 bool flag=0; 63 mi[k]=n+1; 64 mx[k]=0; 65 for(int i=0;i<g1[k].size();++i) 66 if(g1[k][i]!=fa) 67 { 68 flag=1; 69 dfs(g1[k][i],k); 70 s[k]+=s[g1[k][i]]; 71 mi[k]=min(mi[k],mi[g1[k][i]]); 72 mx[k]=max(mx[k],mx[g1[k][i]]); 73 } 74 if(!flag) mi[k]=mx[k]=++m; 75 } 76 int newnode(int u,long long key) 77 { 78 ++tot; 79 H[tot]={u,key,0,0}; 80 return tot; 81 } 82 int merge(int u,int v) 83 { 84 /* 85 merge两个堆u和v 86 这里是采用随机堆，方便合并，方便持久化 87 */ 88 if(!u) return v; 89 if(!v) return u; 90 if(H[v].key<H[u].key) swap(u,v); 91 int k=++tot; 92 H[k]=H[u]; 93 if(rand()%2) H[k].l=merge(H[k].l,v); 94 else H[k].r=merge(H[k].r,v); 95 return k; 96 } 97 void Kshort() 98 { 99 /* 100 求k短路 101 */ 102 dis[T]=0; 103 for(int i=0;i<T;++i) dis[i]=inf; 104 tot=0; 105 for(int i=m-1;i>=0;--i) 106 { 107 /* 108 DAG图求最短路径树 109 */ 110 int fa=0; 111 for(int j=head[i];j!=-1;j=nx[j]) 112 if(dis[i]>e[j].w+dis[e[j].to]) 113 { 114 dis[i]=e[j].w+dis[e[j].to]; 115 pos[i]=j; 116 fa=e[j].to; 117 } 118 rt[i]=rt[fa]; 119 for(int j=head[i];j!=-1;j=nx[j]) 120 if(j!=pos[i]) 121 { 122 //printf("ce : %d %d ",i,e[j].to); 123 rt[i]=merge(rt[i],newnode(e[j].to,e[j].w+dis[e[j].to]-dis[i])); 124 } 125 } 126 //printf("tot : %d ",tot); 127 //printf("len : %d ",len); 128 //printf("m : %d ",m); 129 //for(int i=0;i<=T;++i) printf("%d : %lld %d ",i,dis[i],pos[i]); 130 printf("%lld ",dis[S]+s[1]); 131 heapnode now={0LL,rt[S]}; 132 if(now.root) q.push(now); 133 while(--k&&!q.empty()) 134 { 135 /* 136 每次从优先队列队首取出最小的边集 137 每个边集对对应一种合法的k短路走法 138 有两种扩展方法 139 第一种：将最后一条边换成所在堆的次小元素（相当于从堆里把堆顶删除） 140 第二种：新加一条边，即从最后一条边继续往后走 141 */ 142 now=q.top(); 143 q.pop(); 144 printf("%lld ",now.len+H[now.root].key+dis[S]+s[1]); 145 int id=merge(H[now.root].l,H[now.root].r); 146 //printf("%d : %d %lld ",id,H[id].u,H[id].key); 147 if(id) 148 q.push({now.len,id}); 149 now.len+=H[now.root].key; 150 if(rt[H[now.root].u]) 151 q.push({now.len,rt[H[now.root].u]}); 152 } 153 } 154 int main() 155 { 156 srand(time(0)); 157 scanf("%d%d",&n,&k); 158 for(int i=1;i<=n;++i) scanf("%d",&a[i]); 159 for(int i=1;i<n;++i) 160 { 161 int u,v; 162 scanf("%d%d",&u,&v); 163 g1[u].push_back(v); 164 g1[v].push_back(u); 165 } 166 dfs(1,0); 167 for(int i=0;i<=n+1;++i) head[i]=-1; 168 for(int i=2;i<=n;++i) addedge(mi[i]-1,mx[i],-s[i]); 169 for(int i=0;i<m;++i) addedge(i,i+1,0); 170 S=0,T=m; 171 Kshort(); 172 return 0; 173 }
pbds里的rbtree
1 #include <bits/stdc++.h> 2 #include<ext/pb_ds/assoc_container.hpp> 3 //#include <bits/extc++.h> 4 using namespace std; 5 using namespace __gnu_pbds; 6 typedef tree<int,null_type,less<int>,rb_tree_tag,tree_order_statistics_node_update> rbtree; 7 /* 8 tree中不能有重复元素 9 若要能支持重复元素，可以将int改成pair<int,int> 10 */ 11 int main() 12 { 13 rbtree s; 14 s.insert(8); 15 s.insert(3); 16 s.insert(5); 17 printf("%d ",s.order_of_key(5));//s.order_of_key(x) s里面有多少数<x 18 printf("%d ",*s.find_by_order(1));//s.find_by_order(k) s里面第k+1小的数是谁 19 return 0; 20 }
笛卡尔树
1 pair<int,int> a[maxn+5]; 2 int lson[maxn+5],rson[maxn+5],fa[maxn+5],root; 3 int s[maxn+5]; 4 void buildtreap(pair<int,int> *a) 5 { 6 /* 7 大根堆 8 */ 9 int top=0; 10 for(int i=0;i<=n;++i) rson[i]=lson[i]=fa[i]=0; 11 for(int i=1;i<=n;++i) 12 { 13 int last=0; 14 while(top>=1&&a[i]>a[s[top]]) last=s[top--]; 15 if(top) rson[s[top]]=i,fa[i]=s[top]; 16 lson[i]=last; 17 if(!last) fa[last]=i; 18 s[++top]=i; 19 } 20 root=s[1]; 21 }
图论
匈牙利算法
1 #include<bits/stdc++.h> 2 using namespace std; 3 const int maxn=2000; 4 vector<int> g[maxn+5]; 5 bool v[maxn+5]; 6 int f[maxn+5]; 7 int n,m; 8 void addedge(int u,int v) 9 { 10 g[u].push_back(v); 11 g[v].push_back(u); 12 } 13 bool dfs(int k) 14 { 15 for(int i=0;i<g[k].size();++i) 16 if(v[g[k][i]]==0) 17 { 18 v[g[k][i]]=1; 19 if(f[g[k][i]]==0||dfs(f[g[k][i]])) 20 { 21 f[g[k][i]]=k; 22 f[k]=g[k][i]; 23 return true; 24 } 25 } 26 return false; 27 } 28 int main() 29 { 30 scanf("%d %d",&n,&m); 31 for(int i=0;i<=n+m;++i) g[i].clear(); 32 for(int i=n+1;i<=m+n;++i) 33 { 34 int x,y; 35 scanf("%d %d",&x,&y); 36 ++x,++y; 37 addedge(i,x); 38 addedge(i,y); 39 } 40 memset(f,0,sizeof(f)); 41 int ans=0; 42 for(int i=1;i<=n+m;++i) 43 { 44 if(f[i]!=0) continue; 45 memset(v,0,sizeof(v)); 46 if(dfs(i)) ++ans; 47 } 48 printf("%d",ans); 49 return 0; 50 }
Dinic
1 #include<bits/stdc++.h> 2 using namespace std; 3 const int maxn=2000,inf=1e9,maxm=90000; 4 struct Edge 5 { 6 int from,to,cap,flow; 7 }edge[maxm*2+5];; 8 vector <int> g[maxn+5]; 9 int step[maxn];//从源点到点x的距离 10 int iter[maxn];//定点x的第几条边开始有用 11 int n,m,S,T,len; 12 void addedge(int from,int to,int cap) 13 { 14 ++len; 15 edge[len]={from,to,cap,0}; 16 g[from].push_back(len); 17 ++len; 18 edge[len]={to,from,0,0}; 19 g[to].push_back(len); 20 } 21 void bfs(int S) 22 { 23 memset(step,-1,sizeof(step)); 24 step[S]=0; 25 queue<int> q; 26 q.push(S); 27 while(!q.empty()) 28 { 29 int v=q.front(); 30 q.pop(); 31 for(int i=0;i<g[v].size();++i) 32 { 33 Edge &e=edge[g[v][i]]; 34 if(e.cap>e.flow&&step[e.to]<0) 35 { 36 step[e.to]=step[v]+1; 37 q.push(e.to); 38 } 39 } 40 } 41 } 42 int dfs(int v,int t,int f) 43 { 44 if(v==t) return f; 45 for(int &i=iter[v];i<g[v].size();++i)//这里是引用，i++的同时iter也++，其实相当于上个的used，不过不用判断了 46 { 47 Edge &e=edge[g[v][i]]; 48 if(e.cap>e.flow&&step[e.to]>step[v]) 49 { 50 int d=dfs(e.to,t,min(e.cap-e.flow,f)); 51 if(d>0) 52 { 53 e.flow+=d; 54 edge[g[v][i]^1].flow-=d; 55 return d; 56 } 57 } 58 } 59 return 0; 60 } 61 int maxflow(int S,int T) 62 { 63 int flow=0; 64 for(;;) 65 { 66 bfs(S); 67 if(step[T]<0) return flow; 68 memset(iter,0,sizeof(iter)); 69 int f; 70 while((f=dfs(S,T,inf))>0) 71 flow+=f; 72 } 73 } 74 int main() 75 { 76 scanf("%d%d",&n,&m); 77 S=0,T=n+1; 78 for(int i=0;i<=T;++i) g[i].clear(); 79 len=-1; 80 for(int i=1;i<=n;++i) 81 { 82 int x; 83 scanf("%d",&x); 84 if(x==0) addedge(S,i,1);else addedge(i,T,1); 85 } 86 for(int i=0;i<m;++i) 87 { 88 int u,v; 89 scanf("%d%d",&u,&v); 90 addedge(u,v,1); 91 addedge(v,u,1); 92 } 93 printf("%d ",maxflow(S,T)); 94 return 0; 95 }
最大费用最大流（可行流）
1 /* 2 千万要注意初始化len=-1!!!!!!!! 3 */ 4 #include<bits/stdc++.h> 5 using namespace std; 6 const int maxn=200,maxm=5000,inf=1e7; 7 const double eps=1e-12; 8 struct wjmzbmr 9 { 10 int from,to,cap,flow; 11 double cost; 12 }edge[8*maxm]; 13 vector<int> g[maxn+50]; 14 double d[maxn+50]; 15 int a[maxn+50],b[maxn+50],last[maxn+50],f[maxn+50]; 16 int n,m,t,len,S,T; 17 bool v[maxn+50]; 18 void add(int from,int to,int cap,double cost) 19 { 20 edge[++len]=(wjmzbmr){from,to,cap,0,cost},g[from].push_back(len); 21 edge[++len]=(wjmzbmr){to,from,0,0,-cost},g[to].push_back(len); 22 } 23 bool spfa(int S,int T,int &flow,double &cost) 24 { 25 for(int i=0;i<=n+1;++i) d[i]=(double)-inf,f[i]=inf; 26 memset(v,0,sizeof(v)); 27 d[S]=0.0,v[S]=1,last[S]=0,f[S]=inf; 28 queue<int> q; 29 while(!q.empty()) q.pop(); 30 q.push(S); 31 while(!q.empty()) 32 { 33 int u=q.front(); 34 q.pop(); 35 v[u]=0; 36 for(int i=0;i<g[u].size();++i) 37 { 38 wjmzbmr& e=edge[g[u][i]]; 39 if(e.cap>e.flow&&(d[e.to]+eps<d[u]+e.cost)) 40 { 41 d[e.to]=d[u]+e.cost; 42 last[e.to]=g[u][i]; 43 f[e.to]=min(f[u],e.cap-e.flow); 44 if(!v[e.to]) 45 { 46 q.push(e.to); 47 v[e.to]=1; 48 } 49 } 50 } 51 } 52 if(d[T]==-inf) return false; 53 //if(d[T]*f[T]<0) return false;最大费用可行流 54 flow+=f[T]; 55 cost+=d[T]*f[T]; 56 int u=T; 57 while(u!=S) 58 { 59 edge[last[u]].flow+=f[T]; 60 edge[last[u]^1].flow-=f[T]; 61 u=edge[last[u]].from; 62 } 63 return true; 64 } 65 int main() 66 { 67 scanf("%d",&t); 68 while(t--) 69 { 70 scanf("%d %d",&n,&m); 71 len=-1; 72 for(int i=0;i<=n;++i) g[i].clear(); 73 for(int i=1;i<=n;++i) scanf("%d %d",&a[i],&b[i]); 74 for(int i=1;i<=m;++i) 75 { 76 int x,y,z;double p; 77 scanf("%d %d %d %lf",&x,&y,&z,&p); 78 p=log2(1.0-p); 79 if(z==0) continue; 80 else 81 { 82 add(x,y,1,0); 83 add(x,y,z-1,p); 84 } 85 } 86 S=0,T=n+1; 87 for(int i=1;i<=n;++i) 88 { 89 if(a[i]-b[i]==0) continue; 90 if(a[i]-b[i]>0) add(S,i,a[i]-b[i],0); 91 else add(i,T,b[i]-a[i],0); 92 } 93 int flow=0;double cost=0.0; 94 while(spfa(S,T,flow,cost)) ; 95 printf("%.2f ",1.0-pow(2.0,cost)); 96 } 97 return 0; 98 }
最小费用最大流(Dijkstra实现)
1 #include<bits/stdc++.h> 2 using namespace std; 3 #define mp make_pair 4 const int maxn=2000+5,inf=1e9; 5 pair<int,int> a[maxn+5]; 6 struct edge 7 { 8 int to,cap,flow,cost,rev; 9 }; 10 vector<edge> g[maxn+5]; 11 struct heapnode 12 { 13 int id,v; 14 bool operator < (const heapnode &t) const 15 { 16 return v>t.v; 17 } 18 }; 19 priority_queue<heapnode> q; 20 int n,m,S,T,flow,cost; 21 int h[maxn+5],dis[maxn+5],cur[maxn+5]; 22 bool ins[maxn+5]; 23 void addedge(int from,int to,int cap,int cost) 24 { 25 g[from].push_back({to,cap,0,cost,(int)g[to].size()}); 26 g[to].push_back({from,0,0,-cost,(int)g[from].size()-1}); 27 } 28 bool Dijkstra() 29 { 30 for (int i=0;i<=T;++i) 31 { 32 h[i]=min(h[i]+dis[i],inf); 33 dis[i]=i==S?0:inf; 34 } 35 q.push(heapnode{S,0}); 36 while (!q.empty()) 37 { 38 heapnode x=q.top(); 39 q.pop(); 40 if(x.v>dis[x.id]) continue; 41 for(int i=0;i<g[x.id].size();++i) 42 { 43 edge e=g[x.id][i]; 44 if(e.cap>e.flow&&x.v+h[x.id]+e.cost-h[e.to]<dis[e.to]) 45 { 46 dis[e.to]=x.v+h[x.id]+e.cost-h[e.to]; 47 q.push(heapnode{e.to, dis[e.to]}); 48 } 49 } 50 } 51 return dis[T]<inf; 52 } 53 int dfs(int x,int a) 54 { 55 if(x == T) return a; 56 int ans = 0; 57 ins[x] = 1; 58 for(int &i=cur[x];i<g[x].size();++i) 59 { 60 edge &e=g[x][i]; 61 if(!ins[e.to]&&dis[x]+h[x]+e.cost-h[e.to]==dis[e.to]&&e.cap>e.flow) 62 { 63 int now=dfs(e.to,min(a,e.cap-e.flow)); 64 e.flow += now; 65 g[e.to][e.rev].flow-=now; 66 ans+=now; 67 a-=now; 68 if (!a) break; 69 } 70 } 71 ins[x] = 0; 72 return ans; 73 } 74 int main() 75 { 76 freopen("ce.in","r",stdin); 77 scanf("%d%d",&n,&m); 78 for(int i=1;i<=n;++i) scanf("%d%d",&a[i].first,&a[i].second); // bottle 79 for(int i=n+1;i<=n+m;++i) scanf("%d%d",&a[i].first,&a[i].second); // people 80 scanf("%d%d",&a[n+m+1].first,&a[n+m+1].second); 81 S=0,T=n+m+2; 82 for(int i=n+1;i<=n+m;++i) addedge(S,i,1,0); 83 addedge(S,n+m+1,n-1,0); 84 for(int i=1;i<=n;++i) addedge(i,T,1,0); 85 for(int i=n+1;i<=n+m+1;++i) 86 for(int j=1;j<=n;++j) 87 addedge(i,j,1,abs(a[i].first-a[j].first)+abs(a[i].second-a[j].second)); 88 while(Dijkstra()) 89 { 90 memset(cur,0,sizeof(cur)); 91 int now=dfs(S,inf); 92 flow+=now; 93 cost+=now*(dis[T]+h[T]-h[S]); 94 } 95 for(int i=1;i<=n;++i) cost+=abs(a[i].first-a[n+m+1].first)+abs(a[i].second-a[n+m+1].second); 96 printf("%d",cost); 97 return 0; 98 }
KM算法
1 #include<bits/stdc++.h> 2 using namespace std; 3 typedef int ll; 4 const int maxn=407,inf=1e9; 5 int nl,nr,m; 6 struct KuhnMunkres 7 { 8 int n;//左边1~n个点，右边1~n个点 9 int a[maxn+5][maxn+5]; 10 int lx[maxn+5],ly[maxn+5],sla[maxn+5];//lx是左边顶标 ly是右边顶标 11 int fl[maxn+5],fr[maxn+5];//fl[i]表示左边第i个点匹配右边哪个点 fr[i]表示右边第i个点匹配哪个点 12 int vx[maxn+5],vy[maxn+5],pre[maxn+5]; 13 int q[maxn+5],tp; 14 void match(int x) 15 { 16 while(x) 17 { 18 fr[x]=pre[x]; 19 int y=fl[fr[x]]; 20 fl[fr[x]]=x; 21 x=y; 22 } 23 } 24 void find(int x) 25 { 26 fill(vx,vx+n+1,0); 27 fill(vy,vy+n+1,0); 28 fill(sla,sla+n+1,inf); 29 q[tp=1]=x;vx[x]=1; 30 while(1) 31 { 32 for(int i=1;i<=tp;i++) 33 { 34 int x=q[i]; 35 for(int y=1;y<=n;y++) 36 { 37 int t=lx[x]+ly[y]-a[x][y]; 38 if(vy[y]||t>sla[y])continue; 39 pre[y]=x; 40 if(t==0) 41 { 42 if(!fr[y]){match(y);return;} 43 q[++tp]=fr[y];vy[y]=1;vx[fr[y]]=1; 44 } 45 else sla[y]=t; 46 } 47 } 48 int d=inf;tp=0; 49 for(int i=1;i<=n;i++)if(!vy[i]&&d>sla[i])d=sla[i],x=i; 50 for(int i=1;i<=n;i++) 51 { 52 if(vx[i])lx[i]-=d; 53 if(vy[i])ly[i]+=d; 54 else sla[i]-=d; 55 } 56 if(!fr[x]){match(x);return;} 57 q[++tp]=fr[x];vy[x]=vx[fr[x]]=1; 58 } 59 } 60 void solve() 61 { 62 memset(lx,0,sizeof(lx)); 63 memset(ly,0,sizeof(ly)); 64 memset(fl,0,sizeof(fl)); 65 memset(fr,0,sizeof(fr)); 66 for(int i=1;i<=n;++i) lx[i]=*max_element(a[i]+1,a[i]+n+1); 67 for(int i=1;i<=n;++i) find(i); 68 } 69 }km; 70 int main() 71 { 72 scanf("%d%d%d",&nl,&nr,&m); 73 while(m--) 74 { 75 int u,v; 76 scanf("%d%d",&u,&v); 77 scanf("%d",&km.a[u][v]); 78 } 79 km.n=max(nl,nr); 80 km.solve(); 81 long long ans=0; 82 for(int i=1;i<=nl;++i)ans+=km.a[i][km.fl[i]]; 83 printf("%lld ",ans); 84 for(int i=1;i<=nl;++i) 85 if(km.a[i][km.fl[i]]==0) printf("0 "); 86 else 87 printf("%d ",km.fl[i]); 88 return 0; 89 }
有向图强连通分量（两次DFS）
1 #include<cstring> 2 #include<cstdio> 3 #include<algorithm> 4 #include<vector> 5 using namespace std; 6 const int maxn=1e5,inf=1e9; 7 vector < int > g[maxn+50],g1[maxn+50],g2[maxn+50]; 8 bool v[maxn+50]; 9 int n,m,color[maxn+50],t[maxn+50],len,c,d[maxn+50]; 10 void dfs(int k) 11 { 12 v[k]=1; 13 for(int i=0;i<g[k].size();++i) 14 if(!v[g[k][i]]) dfs(g[k][i]); 15 ++len; 16 t[len]=k; 17 } 18 void dfs1(int k) 19 { 20 v[k]=1; 21 color[k]=c; 22 for(int i=0;i<g1[k].size();++i) 23 if(!v[g1[k][i]]) dfs1(g1[k][i]); 24 } 25 bool check2(int a,int b) 26 { 27 if(a==b) return 0; 28 for(int i=0;i<g2[a].size();++i) 29 if(g2[a][i]==b) return 0; 30 return 1; 31 } 32 int main() 33 { 34 35 while(scanf("%d %d",&n,&m)!=EOF) 36 { 37 c=0; 38 len=0; 39 memset(d,0,sizeof(d)); 40 for(int i=0;i<=n;++i) g[i].clear(),g1[i].clear(),g2[i].clear(); 41 for(int i=1;i<=m;++i) 42 { 43 int x,y; 44 scanf("%d %d",&x,&y); 45 g[x].push_back(y); 46 g1[y].push_back(x); 47 } 48 memset(v,0,sizeof(v)); 49 memset(t,0,sizeof(t)); 50 for(int i=1;i<=n;++i) if(!v[i]) dfs(i); 51 memset(v,0,sizeof(v)); 52 for(int i=len;i>=1;--i) if(!v[t[i]]) ++c,dfs1(t[i]); 53 for(int i=1;i<=n;++i) 54 for(int j=0;j<g[i].size();++j) 55 if(check2(color[i],color[g[i][j]])) 56 g2[color[i]].push_back(color[g[i][j]]),++d[color[i]]; 57 int s=0,x=0; 58 for(int i=1;i<=c;++i) if(d[i]==0) ++s,x=i; 59 if(s>1) printf("0 "); 60 else 61 { 62 int ans=0; 63 for(int i=1;i<=n;++i) 64 if(color[i]==x) ++ans; 65 printf("%d ",ans); 66 } 67 } 68 return 0; 69 } 70 /*bitset优化*/ 71 void dfs(int k) 72 { 73 if(vis[k]==0) return; 74 vis.reset(k); 75 bitset<maxn+5> nx=vis&g[k]; 76 int u=nx._Find_first(); 77 while(u<=n) 78 { 79 dfs(u); 80 u=nx._Find_next(u); 81 } 82 a.push_back(k); 83 } 84 void dfs1(int k) 85 { 86 if(vis[k]==0) return; 87 vis.reset(k); 88 bitset<maxn+5> nx=vis&g1[k]; 89 int u=nx._Find_first(); 90 while(u<=n) 91 { 92 dfs1(u); 93 u=nx._Find_next(u); 94 } 95 }
无向图点双连通分量
1 #include<bits/stdc++.h> 2 using namespace std; 3 const int maxn=1000,maxm=1000000; 4 vector<int> g[maxn+50],bcc[maxn+50]; 5 int dfstime[maxn+50],low[maxn+50],color[maxn+50],head[maxn+50]; 6 bool p[maxn+5][maxn+5],map[maxn+5][maxn+5],flag; 7 int n,m,top,c,t; 8 bool ins[maxn+50]; 9 struct wjmzbmr 10 { 11 int u,v; 12 }s[maxm]; 13 void tarjan(int k,int fa) 14 { 15 low[k]=dfstime[k]=++t; 16 for(int i=0;i<g[k].size();++i) 17 { 18 int u=g[k][i]; 19 if(u==fa) continue; 20 if(dfstime[u]) 21 if(dfstime[u]<dfstime[k]) low[k]=min(low[k],dfstime[u]),s[++top]={k,u}; 22 else; 23 else 24 { 25 s[++top]={k,u}; 26 tarjan(u,k); 27 low[k]=min(low[k],low[u]); 28 if(low[u]>=dfstime[k]) 29 { 30 ++c; 31 while(1) 32 { 33 wjmzbmr e=s[top--]; 34 if(color[e.u]!=c) 35 { 36 bcc[c].push_back(e.u); 37 color[e.u]=c; 38 } 39 if(color[e.v]!=c) 40 { 41 bcc[c].push_back(e.v); 42 color[e.v]=c; 43 } 44 if(e.u==k&&e.v==u) break; 45 } 46 } 47 } 48 } 49 } 50 int main() 51 { 52 scanf("%d %d",&n,&m); 53 while(!(n==0&&m==0)) 54 { 55 for(int i=0;i<=n;++i) g[i].clear(),bcc[i].clear(); 56 memset(p,0,sizeof(p)); 57 for(int i=1;i<=m;++i) 58 { 59 int x,y; 60 scanf("%d %d",&x,&y); 61 p[x][y]=p[y][x]=1; 62 } 63 for(int i=1;i<=n;++i) 64 for(int j=1;j<=n;++j) 65 if(i!=j&&!p[i][j]) g[i].push_back(j); 66 memset(dfstime,0,sizeof(dfstime)); 67 memset(s,0,sizeof(s)); 68 memset(color,0,sizeof(color)); 69 top=c=t=0; 70 for(int i=1;i<=n;++i) 71 if(!dfstime[i]) tarjan(i,-1); 72 scanf("%d %d",&n,&m); 73 } 74 return 0; 75 }
无向图边双连通分量
1 #include<bitsstdc++.h> 2 using namespace std; 3 const int maxn=4e5; 4 struct wjmzbmr 5 { 6 int x,y,pos; 7 }; 8 wjmzbmr e[2*maxn+50]; 9 vector<int> g[maxn+50]; 10 int dfstime[maxn+50],low[maxn+50],s[maxn+50],color[maxn+50],num[maxn+50]; 11 int n,m,top,c,t,ans1=0; 12 bool ins[maxn+50]; 13 void tarjan(int k,int fa) 14 { 15 low[k]=dfstime[k]=++t; 16 s[++top]=k; 17 ins[k]=1; 18 for(int i=0;i<g[k].size();++i) 19 { 20 int u=e[g[k][i]].y; 21 if(fa==u) continue; 22 if(!dfstime[u]) 23 { 24 tarjan(u,k); 25 //if(low[u]>dfstime[k]) p[g[k][i]]=1; 判断此边是否为桥 26 low[k]=min(low[k],low[u]); 27 } 28 else 29 if(ins[u]) low[k]=min(low[k],low[u]); 30 } 31 if(dfstime[k]==low[k]) 32 { 33 ++c; 34 while(1) 35 { 36 int v=s[top--]; 37 ins[v]=0,color[v]=c,++num[c]; 38 if(v==k) break; 39 } 40 } 41 } 42 int main() 43 { 44 scanf("%d %d",&n,&m); 45 for(int i=0;i<=n;++i) g[i].clear(); 46 for(int i=1;i<=m;++i) 47 { 48 int x,y; 49 scanf("%d %d",&x,&y); 50 g[x].push_back(2*i-1),g[y].push_back(2*i); 51 e[2*i-1]={x,y,i},e[2*i]={y,x,i}; 52 } 53 memset(dfstime,0,sizeof(dfstime)); 54 for(int i=1;i<=n;++i) low[i]=n+1; 55 memset(color,0,sizeof(color)); 56 memset(num,0,sizeof(num)); 57 memset(s,0,sizeof(s)); 58 memset(ins,0,sizeof(ins)); 59 top=c=t=0; 60 tarjan(1,-1); 61 return 0; 62 }
图的遍历（非递归）
1 void dfs(int k,int last) 2 { 3 /* L[k]=++t; 4 deep[k]=deep[last]+1; 5 fa[k][0]=last; 6 for(int i=0;i<g[k].size();++i) 7 if(g[k][i]!=last) dfs(g[k][i],k); 8 R[k]=t;*/ 9 while(!s.empty()) s.pop(); 10 memset(head,0,sizeof(head));//head[i]表示第i个点当前遍历到了第几个相邻点 11 s.push(0); 12 s.push(1); 13 while(s.size()>1) 14 { 15 int k=s.top(); 16 s.pop(); 17 int last=s.top(); 18 s.push(k); 19 if(!head[k]) 20 { 21 L[k]=++t; 22 deep[k]=deep[last]+1; 23 fa[k][0]=last; 24 } 25 if(head[k]<g[k].size()) 26 if(g[k][head[k]]==last) ++head[k]; 27 if(head[k]==g[k].size()) 28 { 29 R[k]=t; 30 s.pop(); 31 } 32 else 33 s.push(g[k][head[k]++]); 34 } 35 }
字符串
EXKMP
1 #include<cstring> 2 #include<algorithm> 3 #include<cstdio> 4 using namespace std; 5 const int maxn=1e6; 6 char S[maxn+50],T[maxn+50];//S是母串，T是子串 7 int len1,len2; 8 int next[maxn+50],extend[maxn+50];//extend[i]表示S[i..len1-1]和T的最长公共前缀的长度，next[i]表示T[i..len2-1]和T的最长公共前缀的长度 9 void getnext() 10 { 11 next[0]=len2; 12 int j=0; 13 while(j+1<len2&&T[j]==T[j+1]) ++j; 14 next[1]=j; 15 int k=1; 16 for(int i=2;i<len2;++i) 17 { 18 int p=k+next[k]-1,l=next[i-k]; 19 if(i+l<p+1) next[i]=l; 20 else 21 { 22 j=max(p-i+1,0); 23 while(i+j<len2&&T[i+j]==T[j]) ++j; 24 next[i]=j; 25 k=i; 26 } 27 } 28 } 29 void ekmp() 30 { 31 int j=0; 32 while(j<len1&&j<len2&&S[j]==T[j]) ++j; 33 extend[0]=j; 34 int k=0; 35 for(int i=1;i<len1;++i) 36 { 37 int p=k+extend[k]-1,l=next[i-k];//p表示到达的最远位置，k是对应最远位置的i 38 if(i+l<p+1) extend[i]=l; 39 else 40 { 41 j=max(p-i+1,0); 42 while(i+j<len1&&j<len2&&S[i+j]==T[j]) ++j; 43 extend[i]=j; 44 k=i; 45 } 46 } 47 } 48 int main() 49 { 50 scanf("%s%s",S,T); 51 len1=strlen(S),len2=strlen(T); 52 getnext(); 53 ekmp(); 54 for(int i=0;i<len1;++i) printf("%d ",extend[i]); 55 return 0; 56 }
Manacher
1 #include<cstring> 2 #include<algorithm> 3 #include<cstdio> 4 using namespace std; 5 const int maxn=10000; 6 char s[maxn*2+50]; 7 int len=1,mx=0,id=0; 8 int p[2*maxn+50];//p[i]表示以i为中心，向两边扩展的最长长度 9 char c; 10 int main() 11 { 12 s[0]='$',s[1]='#'; 13 while(scanf("%c",&c)==1) s[++len]=c,s[++len]='#';//在原字符串每个中间插上#，包括头尾，使得回文串长度为奇数，同时为了防止越界，第一个字符设为$ 14 for(int i=0;i<=len;++i) 15 { 16 if(mx>i) p[i]=min(p[2*id-i],mx-i);else p[i]=1; 17 while(s[i+p[i]]==s[i-p[i]]) ++p[i]; 18 if(i+p[i]>mx) 19 { 20 mx=i+p[i]; 21 id=i; 22 } 23 }//O(n)求p数组 24 for(int i=0;i<=len;++i) printf("%d ",p[i]); 25 return 0; 26 }
后缀数组
1 #include<bits/stdc++.h> 2 using namespace std; 3 const int maxn=2e5; 4 char s[maxn+50]; 5 int sa[maxn+50],rk[maxn+50],height[maxn+50],d[maxn+50][30]; 6 int t[maxn+50],t2[maxn+50],c[maxn+50]; 7 int len,k; 8 void getsa(int m)//m表示最大字符的编码 9 { 10 memset(t,-1,sizeof(t)); 11 memset(t2,-1,sizeof(t2)); 12 int *x=t,*y=t2; 13 for(int i=0;i<m;++i) c[i]=0; 14 for(int i=0;i<len;++i) c[x[i]=s[i]]++; 15 for(int i=1;i<m;++i) c[i]+=c[i-1]; 16 for(int i=len-1;i>=0;--i) sa[--c[x[i]]]=i; 17 for(int k=1;k<=len;k<<=1) 18 { 19 int p=0; 20 for(int i=len-k;i<len;++i) y[p++]=i; 21 for(int i=0;i<len;++i) if(sa[i]>=k) y[p++]=sa[i]-k; 22 for(int i=0;i<m;++i) c[i]=0; 23 for(int i=0;i<len;++i) c[x[y[i]]]++; 24 for(int i=0;i<m;++i) c[i]+=c[i-1]; 25 for(int i=len-1;i>=0;--i) sa[--c[x[y[i]]]]=y[i]; 26 swap(x,y); 27 p=1,x[sa[0]]=0; 28 for(int i=1;i<len;++i) 29 if(y[sa[i-1]]==y[sa[i]]&&y[sa[i-1]+k]==y[sa[i]+k]) x[sa[i]]=p-1;else x[sa[i]]=p++; 30 if(p>=len) break; 31 m=p; 32 } 33 } 34 void getheight() 35 { 36 int k=0; 37 for(int i=0;i<len;++i) rk[sa[i]]=i; 38 for(int i=0;i<len;++i) 39 { 40 if(k) --k; 41 if(rk[i]==0) continue; 42 int j=sa[rk[i]-1]; 43 while(s[i+k]==s[j+k]) ++k; 44 height[rk[i]]=k; 45 } 46 } 47 void rmq_init() 48 { 49 for(int i=0;i<len;i++) d[i][0]=height[i]; 50 for(int j=1;(1<<j)-1<=len;j++) 51 for(int i=0;i+(1<<j)-1<len;i++) 52 d[i][j]=min(d[i][j-1],d[i+(1<<(j-1))][j-1]); 53 } 54 int lcp(int l,int r) 55 { 56 if(l<0||r>=len) return 0; 57 l=rk[l],r=rk[r]; 58 if(l>r) swap(l,r); 59 ++l; 60 int k=0; 61 while((1<<(k+1))<=r-l+1) k++; 62 return min(d[l][k],d[r-(1<<k)+1][k]); 63 } 64 int main() 65 { 66 scanf("%s",s); 67 len=strlen(s); 68 getsa('z'+1); 69 getheight(); 70 rmq_init(); 71 int ans=0; 72 for(int l=1;l<=len;++l) 73 for(int i=0;i<len;i+=l) 74 { 75 int m=lcp(i,i+l); 76 ans=max(ans,m/l+1); 77 ans=max(ans,lcp(i-l+m%l,i+m%l)/l+1); 78 } 79 printf("%d ",ans); 80 // for(int i=0;i<n;++i) printf("%d ",sa[i]);printf(" "); 81 // for(int i=0;i<n;++i) printf("%d ",rk[i]);printf(" "); 82 // for(int i=0;i<n;++i) printf("%d ",height[i]);printf(" "); 83 return 0; 84 85 }
AC自动机
1 #include<bits/stdc++.h> 2 using namespace std; 3 const int maxn=1e5; 4 int ch[maxn+50][26]; 5 int sz=0,n,root=0; 6 char s[maxn+50],S[maxn+50]; 7 int danger[maxn+50]; 8 int nx[maxn+50]; 9 int last[maxn+50]; 10 queue<int> q; 11 bool vis[maxn+50]; 12 void buildtrie(char *s) 13 { 14 int u=root; 15 int len=strlen(s); 16 for(int i=0;i<len;++i) 17 { 18 int id=s[i]-'a'; 19 if(!ch[u][id]) ch[u][id]=++sz; 20 u=ch[u][id]; 21 } 22 ++danger[u]; 23 } 24 void buildfail() 25 { 26 while(!q.empty()) q.pop(); 27 for(int i=0;i<26;++i) if(ch[root][i]) q.push(ch[root][i]),nx[ch[root][i]]=root,last[ch[root][i]]=0; 28 while(!q.empty()) 29 { 30 int u=q.front(); 31 q.pop(); 32 for(int i=0;i<26;++i) 33 if(ch[u][i]) 34 { 35 int k=nx[u]; 36 while(!ch[k][i]&&k) k=nx[k]; 37 nx[ch[u][i]]=ch[k][i]; 38 //danger[ch[u][i]]|=danger[ch[k][i]]; 39 if(danger[ch[k][i]]) last[ch[u][i]]=ch[k][i];else last[ch[u][i]]=last[ch[k][i]]; 40 q.push(ch[u][i]); 41 } 42 else ch[u][i]=u==0?0:ch[nx[u]][i]; 43 } 44 } 45 bool query(char *s) 46 { 47 int len=strlen(s); 48 int u=root; 49 int ans=0; 50 for(int i=0;i<len;++i) 51 { 52 int id=s[i]-'a'; 53 u=ch[u][id]; 54 int v=u; 55 while(v) 56 { 57 if(vis[v]) break; 58 vis[v]=1; 59 if(danger[v]) ans+=danger[v],danger[v]=0; 60 v=last[v]; 61 } 62 } 63 return ans==n; 64 } 65 void init() 66 { 67 for(int i=0;i<=sz;++i) memset(ch[i],0,sizeof(ch[i])); 68 for(int i=0;i<=sz;++i) danger[i]=0; 69 for(int i=0;i<=sz;++i) last[i]=0; 70 for(int i=0;i<=sz;++i) nx[i]=0; 71 for(int i=0;i<=sz;++i) vis[i]=0; 72 sz=0; 73 } 74 int main() 75 { 76 int T; 77 scanf("%d",&T); 78 while(T--) 79 { 80 scanf("%d",&n); 81 init(); 82 int mx=0; 83 for(int i=1;i<=n;++i) 84 { 85 scanf("%s",s); 86 int l=strlen(s); 87 if(l>mx) 88 { 89 mx=l; 90 strcpy(S,s); 91 } 92 buildtrie(s); 93 } 94 buildfail(); 95 if(query(S)) printf("%s ",S);else printf("No "); 96 } 97 return 0; 98 }
SAM（后缀自动机）
1 #include<bits/stdc++.h> 2 using namespace std; 3 const int maxn=1e6; 4 char s[maxn+5]; 5 int sz=0,len,last=0; 6 int maxlen[2*maxn+50],minlen[2*maxn+50],trans[2*maxn+50][26],slink[2*maxn+50],endpos[2*maxn+50]; 7 int sa[maxn*2+5],cnt[maxn*2+5]; 8 long long ans[maxn*2+50]; 9 int suf[maxn*2+5]; 10 int g[maxn*2+5][26]; 11 bool vis[maxn*2+5]; 12 int build(int _maxlen,int _minlen,int* _trans,int _slink) 13 { 14 maxlen[++sz]=_maxlen; 15 minlen[sz]=_minlen; 16 for(int i=0;i<26;++i) 17 if(_trans==NULL) trans[sz][i]=-1;else trans[sz][i]=_trans[i]; 18 slink[sz]=_slink; 19 return sz; 20 } 21 int add(char ch,int u) 22 { 23 int c=ch-'a'; 24 if(trans[u][c]!=-1&&maxlen[u]+1==maxlen[trans[u][c]]) return trans[u][c]; //顺着Trie树走 25 int z=build(maxlen[u]+1,-1,NULL,-1); 26 int v=u; 27 while(v!=-1&&trans[v][c]==-1) trans[v][c]=z,v=slink[v]; 28 if(v==-1)//最简单的情况，suffix-path(u->S)上都没有对应字符ch的转移 29 { 30 minlen[z]=1; 31 slink[z]=0; 32 return z; 33 } 34 int x=trans[v][c]; 35 if(maxlen[v]+1==maxlen[x])//较简单的情况，不用拆分x 36 { 37 minlen[z]=maxlen[x]+1; 38 slink[z]=x; 39 return z; 40 } 41 int y=build(maxlen[v]+1,-1,trans[x],slink[x]); //最复杂的情况，拆分x，y表示<=maxlen[v]+1的那段 42 slink[y]=slink[x]; 43 minlen[x]=maxlen[y]+1; 44 slink[x]=y; 45 minlen[z]=maxlen[y]+1; 46 slink[z]=y; 47 int w=v; 48 while(w!=-1&&trans[w][c]==x) trans[w][c]=y,w=slink[w]; 49 minlen[y]=maxlen[slink[y]]+1; 50 return z; 51 } 52 void addedge(int u,int v,int w) 53 { 54 g[u][w]=v; 55 } 56 void goup(int k,int start,int step) 57 { 58 if(k==0||vis[k]) return; 59 vis[k]=1; 60 int father=slink[k]; 61 step-=maxlen[k]-minlen[k]+1; 62 addedge(father,k,s[start+step]-'a'); 63 goup(father,start,step); 64 } 65 int main() 66 { 67 int T; 68 scanf("%d",&T); 69 while(T--) 70 { 71 scanf("%s",s); 72 len=strlen(s); 73 for(int i=0;i<=2*len;++i) 74 { 75 for(int j=0;j<26;++j) trans[i][j]=-1; 76 slink[i]=-1; 77 maxlen[i]=minlen[i]=0; 78 memset(g[i],0,sizeof(g[i])); 79 suf[i]=-1; 80 } 81 for(int i=0;i<26;++i) trans[0][i]=slink[0]=-1; 82 maxlen[0]=minlen[0]=0; 83 sz=0; 84 last=0; 85 memset(endpos,0,sizeof(endpos)); 86 for(int i=0;i<len;++i) 87 { 88 last=add(s[i],last);//后缀树要倒过来建 89 suf[last]=i; 90 endpos[last]=1; 91 } 92 for(int i=1;i<=sz;++i) ++cnt[maxlen[i]]; 93 for(int i=2;i<=sz;++i) cnt[i]+=cnt[i-1]; 94 for(int i=1;i<=sz;++i) sa[cnt[maxlen[i]]--]=i; 95 for(int i=sz;i>=1;--i) 96 { 97 int u=sa[i]; 98 endpos[slink[u]]+=endpos[u]; 99 } 100 101 for(int i=0;i<=sz;++i) vis[i]=0; 102 for(int i=1;i<=sz;++i) 103 if(suf[i]!=-1) goup(i,suf[i],len-suf[i]); 104 105 106 memset(ans,0,sizeof(ans)); 107 for(int i=1;i<=sz;++i) ans[maxlen[i]]=max(ans[maxlen[i]],(long long)endpos[i]); 108 for(int i=len-1;i>=1;--i) ans[i]=max(ans[i],ans[i+1]); 109 for(int i=1;i<=len;++i) printf("%lld ",ans[i]); 110 } 111 return 0; 112 }
回文树
1 #include<bits/stdc++.h> 2 using namespace std; 3 typedef long long ll; 4 const ll mod=1e9+7; 5 const int maxn=2000000; 6 ll ans; 7 ll pw[maxn+5]; 8 char s[maxn+5]; 9 struct Palindromic_Tree 10 { 11 int ch[maxn+5][11];//next指针，next指针和字典树类似，指向的串为当前串两端加上同一个字符构成 12 int fail[maxn+5];//fail指针，失配后跳转到fail指针指向的节点 13 int cnt[maxn+5];//某个点表示的回文子串出现的次数 14 int num[maxn+5];//表示以节点i表示的最长回文串的最右端点为回文串结尾的回文串个数 15 int len[maxn+5];//len[i]表示节点i表示的回文串的长度 16 int last;//指向上一个字符所在的节点，方便下一次add 17 int n;//字符数组指针 18 int p;//节点指针 19 int newnode(int l) 20 { //新建节点 21 for(int i=0;i<=10;++i) ch[p][i]=0; 22 cnt[p]=0; 23 num[p]=0; 24 len[p]=l; 25 return p++; 26 } 27 void init() 28 { //初始化 29 p=0; 30 newnode(0); 31 newnode(-1); 32 last=0; 33 n=0; 34 s[n]=-1;//开头放一个字符集中没有的字符，减少特判 35 fail[0]=1; 36 } 37 int get_fail(int x) 38 { //和KMP一样，失配后找一个尽量最长的 39 while(s[n-len[x]-1]!=s[n]) x=fail[x]; 40 return x; 41 } 42 void add(int c) 43 { 44 c-='0'; 45 s[++n]=c; 46 int cur=get_fail(last);//通过上一个回文串找这个回文串的匹配位置 47 if(!ch[cur][c]) 48 { //如果这个回文串没有出现过，说明出现了一个新的本质不同的回文串 49 int now=newnode(len[cur]+2);//新建节点 50 fail[now]=ch[get_fail(fail[cur])][c];//和AC自动机一样建立fail指针，以便失配后跳转 51 ch[cur][c]=now; 52 num[now]=num[fail[now]]+1; 53 } 54 last=ch[cur][c]; 55 cnt[last]++; 56 } 57 void count () 58 { 59 for ( int i = p - 1 ; i >= 0 ; -- i ) cnt[fail[i]] += cnt[i] ; 60 //父亲累加儿子的cnt，因为如果fail[v]=u，则u一定是v的子回文串！ 61 } 62 void dfs(int k,ll num,int len) 63 { 64 for(int i=0;i<=10;++i) 65 if(ch[k][i]) 66 { 67 ans=(ans+(num*10%mod+i+i*pw[len+1]%mod)%mod)%mod; 68 dfs(ch[k][i],(num*10%mod+i+i*pw[len+1]%mod)%mod,len+2); 69 } 70 } 71 void work() 72 { 73 dfs(0,0LL,0); 74 for(int i=0;i<=10;++i) 75 if(ch[1][i]) 76 { 77 ans+=i; 78 dfs(ch[1][i],1LL*i,1); 79 } 80 } 81 }PT; 82 int main() 83 { 84 85 pw[0]=1; 86 for(int i=1;i<=maxn;++i) pw[i]=pw[i-1]*10%mod; 87 scanf("%s",s+1); 88 int len=strlen(s+1); 89 PT.init(); 90 for(int i=1;i<=len;++i) PT.add(s[i]); 91 PT.work(); 92 printf("%lld ",ans); 93 return 0; 94 }
数论
莫比乌斯反演（筛积性函数+数论分块）
1 #include<bits/stdc++.h> 2 using namespace std; 3 const int maxn=1e5; 4 int mu[maxn+50],prime[maxn+50],sum[maxn+50]; 5 bool f[maxn+50]; 6 int a,b,c,d,k; 7 long long solve(int n,int m) 8 { 9 if(n>m) swap(n,m); 10 long long ans=0; 11 for(int i=1,la=0;i<=n;i=la+1) 12 { 13 la=min(n/(n/i),m/(m/i)); 14 ans+=(long long)(sum[la]-sum[i-1])*(n/i)*(m/i); 15 }//对于n/i * m/i 采用分块求和的根号n+根号m的做法 16 return ans; 17 } 18 int main() 19 { 20 mu[1]=1; 21 memset(f,0,sizeof(f)); 22 f[1]=1; 23 for(int i=2;i<=maxn;++i) 24 { 25 if(!f[i]) 26 { 27 prime[++prime[0]]=i; 28 mu[i]=-1; 29 } 30 for(int j=1;j<=prime[0];++j) 31 { 32 if(i*prime[j]>maxn) break; 33 f[i*prime[j]]=1; 34 if(i%prime[j]==0) 35 { 36 mu[i*prime[j]]=0; 37 break; 38 } 39 else mu[i*prime[j]]=-mu[i]; 40 } 41 }//筛积性函数 42 memset(sum,0,sizeof(sum)); 43 sum[1]=mu[1]; 44 for(int i=2;i<=maxn;++i) sum[i]=sum[i-1]+mu[i]; 45 int T; 46 scanf("%d",&T); 47 for(int cas=1;cas<=T;++cas) 48 { 49 printf("Case %d: ",cas); 50 scanf("%d%d%d%d%d",&a,&b,&c,&d,&k); 51 if(k==0) printf("0 ");else 52 { 53 long long ans=solve(b/k,d/k)-solve((a-1)/k,d/k)-solve(b/k,(c-1)/k)+solve((a-1)/k,(c-1)/k); 54 a=max(a,c),b=min(b,d); 55 long long ans1=solve(b/k,b/k)-solve((a-1)/k,b/k)-solve(b/k,(a-1)/k)+solve((a-1)/k,(a-1)/k); 56 printf("%lld ",ans-ans1/2); 57 } 58 } 59 return 0; 60 }
CRT
1 void gcd(ll a,ll b,ll &d,ll &x,ll &y) 2 { 3 if(!b) 4 { 5 d=a,x=1,y=0; 6 return; 7 } 8 gcd(b,a%b,d,y,x); 9 y-=x*(a/b); 10 } 11 ll inv(ll a,ll n) 12 { 13 ll d,x,y; 14 gcd(a,n,d,x,y); 15 return d==1?(x+n)%n:-1; 16 } 17 ll CRT(int n,ll *a,ll *m) 18 { 19 /*n个方程 x=a[i] (mod m[i]) (0<=i<n)*/ 20 ll M=1,d,y,x=0; 21 for(int i=0;i<n;++i) M*=m[i]; 22 for(int i=0;i<n;++i) 23 { 24 ll w=M/m[i]; 25 gcd(m[i],w,d,d,y); 26 x=(x+y*w*a[i])%M; 27 } 28 return (x+M)%M; 29 }
类欧几里得
1 typedef long long ll; 2 ll solve(ll a,ll b,ll c,ll n) 3 { 4 /* 5 计算 sum_{i=0..n} (ax+b)/c 的值 6 */ 7 if(n<0) return 0; 8 if(c<0) a=-a,b=-b,c=-c; 9 if(a<0) b+=a*n,a=-a; 10 if(b<0) 11 { 12 ll p=(-b)/c+1; 13 //(a*x+b+p*c)/c-p 14 return (solve(a,b+p*c,c,n)-p*(n+1))%mod; 15 } 16 if(a==0) 17 return (b/c)*(n+1)%mod; 18 if(a>=c) 19 return (n*(n+1)/2%mod*(a/c)%mod+solve(a%c,b,c,n))%mod; 20 if(b>=c) 21 return ((b/c)*(n+1)%mod+solve(a,b%c,c,n))%mod; 22 //LL m=(a*n+b)/c%MOD; 23 ll m=(n/c*a+(n%c*a+b)/c); 24 return ((m%mod)*(n%mod)%mod-solve(c,c-b-1,a,m-1)+mod)%mod; 25 } 26 ll cal(ll a,ll b,ll c,ll n) 27 { 28 /* 29 计算y=(ax+b)/c下方整点个数，x=0..n 30 要求a>=0,c>0 31 */ 32 if(b>=0) return (solve(a,b,c,n)+n%mod+1)%mod; 33 ll x=(-b)/a; 34 if(x*a+b<0) ++x; 35 if(n<x) return 0; 36 b+=a*x; 37 return (solve(a,b,c,n-x)+(n-x)%mod+1)%mod; 38 }
杜教筛
1 #include<bits/stdc++.h> 2 #define debug cout<<"debug "<<++debug_num<<" :" 3 #define pb push_back 4 #define mp make_pair 5 #define lson l,m,rt<<1 6 #define rson m+1,r,rt<<1|1 7 #define bit(a,b) ((a>>b)&1) //from 0 8 #define all(x) (x).begin(),(x).end() 9 using namespace std; 10 typedef long long ll; 11 typedef pair<int,int> PII; 12 int debug_num=0; 13 14 const int maxn=5e6+10; 15 bool valid[maxn]; 16 ll phi[maxn]; 17 int mu[maxn]; 18 int ans[maxn/10]; 19 int tot; 20 int up; 21 int m; 22 const int maxm=(1LL<<32)/maxn; 23 ll help1[maxm]; 24 int help2[maxm]; 25 bool vis[maxm]; 26 27 void get_prime(int n) 28 { 29 memset(valid,true,sizeof(valid)); 30 tot=0; 31 phi[1]=mu[1]=1; 32 for(int i=2;i<=n;++i){ 33 if(valid[i]){ 34 ans[++tot]=i; 35 mu[i]=-1; 36 phi[i]=i-1; 37 } 38 for(int j=1;j<=tot && i*ans[j]<=n;++j){ 39 int tp=i*ans[j]; 40 valid[tp]=false; 41 if(i%ans[j]==0){ 42 mu[tp]=0; 43 phi[tp]=phi[i]*ans[j]; 44 break; 45 } 46 else{ 47 mu[tp]=-mu[i]; 48 phi[tp]=phi[i]*(ans[j]-1); 49 } 50 } 51 } 52 for(int i=1;i<=n;++i){ 53 phi[i]=phi[i-1]+phi[i]; 54 mu[i]+=mu[i-1]; 55 } 56 } 57 58 ll get_phi(ll n) 59 { 60 return (n<=up)? phi[n] : help1[m/n]; 61 } 62 63 ll get_mu(ll n) 64 { 65 return (n<=up)? mu[n] : help2[m/n]; 66 } 67 68 void solve(ll n) 69 { 70 int t=m/n; 71 if(n<=up || vis[t]) return ; 72 vis[t]=true; 73 help1[t]=n*(n+1)/2;//单位函数前缀和 74 help2[t]=1;//恒等函数前缀和 75 for(ll l=2,r;l<=n;l=r+1){ 76 r=n/(n/l); 77 solve(n/r); 78 help1[t]=help1[t]-(r-l+1)*get_phi(n/r); 79 help2[t]=help2[t]-(r-l+1)*get_mu(n/r); 80 } 81 } 82 83 int main() 84 { 85 //freopen("in.txt","r",stdin); 86 up=maxn-10; 87 get_prime(up); 88 //cout<<clock()<<endl; 89 //cout<<tot<<endl; 90 //cout<<phi[up]<<endl; 91 //for(int i=1;i<=up;++i) if(mu[i]==0) cout<<"fuck"<<endl; 92 int t,n; 93 cin>>t; 94 while(t--) 95 { 96 cin>>n; 97 m=n; 98 if(n<=up) cout<<phi[n]<<" "<<mu[n]<<endl; 99 else{ 100 memset(vis,0,sizeof(vis));//注意清空 101 solve(n); 102 cout<<help1[1]<<" "<<help2[1]<<endl; 103 } 104 } 105 return 0; 106 }
O(1)快速乘
1 inline ll multi(ll x,ll y,ll mod) 2 { 3 ll tmp=(x*y-(ll)((long double)x/mod*y+1.0e-8)*mod); 4 return tmp<0 ? tmp+mod : tmp; 5 }
min25筛
1 //f(1)=1 2 //f(p)=p^2-p 3 //f(p^e)=(p-1)p^(2e-1) 4 #include<bits/stdc++.h> 5 using namespace std; 6 typedef long long ll; 7 const ll mod=1e9+7,inv6=166666668; 8 ll n,M; 9 //pre[d][i]预处理的是2~i中的p^d的和 p是素数 10 //suf[d][i]预处理2~n/i的p^d的和 p是素数 11 //pre[0][i]预处理的是2~i中p^2-p的和 p是素数 12 //suf[0][i]预处理的是2~n/i中p^2-p的和 p是素数 13 vector<ll> pre[3],suf[3],prime; 14 //res是n/枚举的数 15 ll mul(ll a,ll b) 16 { 17 return a*b%mod; 18 } 19 ll add(ll a,ll b) 20 { 21 return (a+b)%mod; 22 } 23 ll sub(ll a,ll b) 24 { 25 return (a-b)%mod; 26 } 27 ll dfs(ll res,int last,ll f) 28 { 29 ll t=(res>M?suf[0][n/res]:pre[0][res])-pre[0][prime[last]-1]; 30 t%=mod; 31 ll ans=mul(t,f);//需要修改 32 for(int i=last;i<(int)prime.size();++i) 33 { 34 int p=prime[i]; 35 if(1LL*p*p>res) break; 36 for(ll q=p,nres=res,nf=f*p%mod*(p-1)%mod;q*p<=res;q*=p)//nf需要修改 37 { 38 ans=add(ans,dfs(nres/=p,i+1,nf)); //枚举更大的数 39 nf=mul(nf,mul(p,p)); //需要修改，继续枚举当前素数，指数大于1的时候，指数每+1，nf*=p^2 40 ans=add(ans,nf); //指数大于1的时候记上贡献 41 } 42 } 43 return ans; 44 } 45 ll f(ll x) 46 { 47 return x*(1+x)/2%mod; 48 } 49 ll ff(ll x) 50 { 51 return x*(x+1)%mod*(2*x+1)%mod*inv6%mod; 52 } 53 ll solve(ll n) 54 { 55 M=sqrt(n); 56 for(int i=0;i<3;++i) pre[i].clear(),pre[i].resize(M+1); 57 for(int i=0;i<3;++i) suf[i].clear(),suf[i].resize(M+1); 58 prime.clear(); 59 for(int i=1;i<=M;++i) 60 { 61 pre[1][i]=f(i)-1; 62 suf[1][i]=f(n/i)-1; 63 pre[2][i]=ff(i)-1; 64 suf[2][i]=ff(n/i)-1; 65 } 66 for(int p=2;p<=M;++p) 67 { 68 if(pre[1][p]==pre[1][p-1]) continue; 69 prime.push_back(p); 70 const ll q=1LL*p*p,m=n/p,pnt[3]={0,pre[1][p-1],pre[2][p-1]}; 71 const int mid=M/p; 72 const int End=min((ll)M,n/q); 73 for(int i=1;i<=mid;++i) 74 { 75 suf[1][i]=sub(suf[1][i],(suf[1][i*p]-pnt[1])*p%mod); 76 suf[2][i]=sub(suf[2][i],(suf[2][i*p]-pnt[2])*p%mod*p%mod); 77 } 78 for(int i=mid+1;i<=End;++i) 79 { 80 suf[1][i]=sub(suf[1][i],(pre[1][m/i]-pnt[1])*p%mod); 81 suf[2][i]=sub(suf[2][i],(pre[2][m/i]-pnt[2])*p%mod*p%mod); 82 } 83 for(int i=M;i>=q;--i) 84 { 85 pre[1][i]=sub(pre[1][i],(pre[1][i/p]-pnt[1])*p%mod); 86 pre[2][i]=sub(pre[2][i],(pre[2][i/p]-pnt[2])*p%mod*p%mod); 87 } 88 } 89 for(int i=1;i<=M;++i) 90 { 91 pre[0][i]=(pre[2][i]-pre[1][i])%mod; 92 suf[0][i]=(suf[2][i]-suf[1][i])%mod; 93 } 94 prime.push_back(M+1); 95 return n>1?1+dfs(n,0,1):1; 96 } 97 int main() 98 { 99 scanf("%lld",&n); 100 printf("1 %lld ",(solve(n)+mod)%mod); 101 return 0; 102 }
矩阵类
1 struct matrix 2 { 3 int n,a[N][N]; 4 matrix operator *(const matrix&b) const 5 { 6 matrix c;c.n=n;memset(c.a,0,sizeof(c.a)); 7 for (int i=0;i<n;i++) 8 for (int j=0;j<b.n;j++) 9 for (int k=0;k<b.n;k++) 10 c.a[i][j]=(c.a[i][j]+1ll*a[i][k]*b.a[k][j])%P; 11 return c; 12 } 13 }
多项式
FFT
1 #include<bits/stdc++.h> 2 using namespace std; 3 const int maxn=400000; 4 const double pi=acos(-1.0); 5 int l,m,n,num; 6 int a[maxn+50],b[maxn+50]; 7 struct wjmzbmr 8 { 9 double r,i; 10 wjmzbmr(double real=0.0,double image=0.0){r=real;i=image;} 11 wjmzbmr operator + (const wjmzbmr o) 12 { 13 return wjmzbmr(r+o.r,i+o.i); 14 } 15 wjmzbmr operator - (const wjmzbmr o) 16 { 17 return wjmzbmr(r-o.r,i-o.i); 18 } 19 wjmzbmr operator * (const wjmzbmr o) 20 { 21 return wjmzbmr(r*o.r-i*o.i,r*o.i+i*o.r); 22 } 23 }; 24 wjmzbmr x1[maxn+50],x2[maxn+50]; 25 void brc(wjmzbmr *y,int l) 26 { 27 for(int i=1,j=l/2;i<l-1;i++) 28 { 29 if(i<j) swap(y[i],y[j]); 30 int k=l/2; 31 while(j>=k)j-=k,k/=2; 32 if(j<k) j+=k; 33 } 34 } 35 void fft(wjmzbmr *y,int l,double on) 36 { 37 wjmzbmr u,t; 38 brc(y,l); 39 for(int h=2;h<=l;h<<=1) 40 { 41 wjmzbmr wn(cos(on*2*pi/h),sin(on*2*pi/h)); 42 for(int j=0;j<l;j+=h) 43 { 44 wjmzbmr w(1,0); 45 for(int k=j;k<j+h/2;k++) 46 { 47 u=y[k]; 48 t=w*y[k+h/2]; 49 y[k]=u+t; 50 y[k+h/2]=u-t; 51 w=w*wn; 52 } 53 } 54 } 55 if(on==-1)for(int i=0;i<l;i++) y[i].r/=l; 56 } 57 void init(int *a,int n,wjmzbmr *x1,int *b,int m,wjmzbmr *x2) 58 { 59 /*将a[0]~a[n-1]放到x1中，将b[0]~b[m-1]放到x2中*/ 60 l=1; 61 while(l<max(n,m)*2) l<<=1; 62 for(int i=0;i<n;++i) 63 { 64 x1[i].r=a[i]; 65 x1[i].i=0.0; 66 } 67 for(int i=n;i<l;++i)x1[i].r=x1[i].i=0.0; 68 for(int i=0;i<m;++i) 69 { 70 x2[i].r=b[i]; 71 x2[i].i=0.0; 72 } 73 for(int i=m;i<l;i++) x2[i].r=x2[i].i=0.0; 74 } 75 int main() 76 { 77 scanf("%d %d",&n,&m); 78 ++n,++m; 79 for(int i=0;i<n;++i) scanf("%d",&a[i]); 80 for(int i=0;i<m;++i) scanf("%d",&b[i]); 81 init(a,n,x1,b,m,x2); 82 fft(x1,l,1); 83 fft(x2,l,1); 84 for(int i=0;i<l;++i) x1[i]=x1[i]*x2[i]; 85 fft(x1,l,-1); 86 for(int i=0;i<n+m-1;++i) printf("%lld ",(long long)(x1[i].r+0.5)); 87 return 0; 88 }
任意模数FFT
1 #include<bits/stdc++.h> 2 using namespace std; 3 const int maxn=524300,mod=1e9+7; 4 int n; 5 int A[maxn+5],B[maxn+5],C[maxn+5]; 6 int len; 7 int pos[maxn+5]; 8 namespace FFT 9 { 10 /* 11 模任意质数mod的FFT 12 */ 13 const int M=1000; 14 struct comp 15 { 16 long double r,i; 17 comp(long double _r=0,long double _i=0){r=_r;i=_i;} 18 comp operator+(const comp x){return comp(r+x.r,i+x.i);} 19 comp operator-(const comp x){return comp(r-x.r,i-x.i);} 20 comp operator*(const comp x){return comp(r*x.r-i*x.i,r*x.i+i*x.r);} 21 comp conj(){return comp(r,-i);} 22 }A[maxn+5],B[maxn+5]; 23 int a0[maxn+5],b0[maxn+5],a1[maxn+5],b1[maxn+5]; 24 const long double pi=acos(-1.0); 25 void FFT(comp a[],int n,int t) 26 { 27 for(int i=1;i<n;i++) 28 if(i<pos[i])swap(a[i],a[pos[i]]); 29 for(int d=0;(1<<d)<n;d++) 30 { 31 int m=1<<d,m2=m<<1; 32 long double o=pi*2/m2*t;comp _w(cos(o),sin(o)); 33 for(int i=0;i<n;i+=m2) 34 { 35 comp w(1,0); 36 for(int j=0;j<m;j++) 37 { 38 comp&A=a[i+j+m],&B=a[i+j],t=w*A; 39 A=B-t;B=B+t;w=w*_w; 40 } 41 } 42 } 43 if(t==-1)for(int i=0;i<n;i++)a[i].r/=n; 44 } 45 void mul(int*a,int*b,int*c,int len) 46 {//c=a*b 47 for(int i=0;i<len;i++)A[i]=comp(a[i],b[i]); 48 FFT(A,len,1); 49 for(int i=0;i<len;i++) 50 { 51 int j=(len-i)&(len-1); 52 B[i]=(A[i]*A[i]-(A[j]*A[j]).conj())*comp(0,-0.25); 53 } 54 FFT(B,len,-1); 55 for(int i=0;i<len;i++)c[i]=((long long)(B[i].r+0.5))%mod; 56 } 57 //输入两个多项式，求a*b mod mod，保存在c中，c不能为a或b，长度为0~len-1 58 void mulmod(int*a,int*b,int*c,int len) 59 { 60 for(int i=0;i<len;i++)a0[i]=a[i]/M,b0[i]=b[i]/M; 61 mul(a0,b0,a0,len); 62 for(int i=0;i<len;i++) 63 { 64 c[i]=1LL*a0[i]*M*M%mod; 65 a1[i]=a[i]%M,b1[i]=b[i]%M; 66 } 67 mul(a1,b1,a1,len); 68 for(int i=0;i<len;i++) 69 { 70 c[i]=(a1[i]+c[i])%mod,a0[i]=(a0[i]+a1[i])%mod; 71 a1[i]=a[i]/M+a[i]%M,b1[i]=b[i]/M+b[i]%M; 72 } 73 mul(a1,b1,a1,len); 74 for(int i=0;i<len;i++)c[i]=(1LL*M*(a1[i]-a0[i]+mod)+c[i])%mod; 75 } 76 } 77 int main() 78 { 79 for(len=1;len<=n;len<<=1);len<<=1; 80 int x=__builtin_ctz(len)-1; 81 for(int i=0;i<len;i++)pos[i]=pos[i>>1]>>1|((i&1)<<x); 82 for(int i=2;i<=n;i++) A[i]=1LL*pa[n-i]*inv[n-i]%mod; 83 for(int i=2;i<=n;i++) B[i]=1LL*pb[n-i]*inv[n-i]%mod; 84 FFT::mulmod(A,B,C,len); 85 return 0; 86 }
NTT
1 #include<bits/stdc++.h> 2 using namespace std; 3 typedef long long ll; 4 const int maxn=262144*4; 5 //const long long P=50000000001507329LL; // 190734863287 * 2 ^ 18 + 1 6 //const int P=1004535809; // 479 * 2 ^ 21 + 1 7 const ll mod=998244353; // 119 * 2 ^ 23 + 1 8 const ll G=3; 9 10 ll len=0; 11 ll pw[maxn+5],pwinv[maxn+5]; 12 ll A[maxn+5],B[maxn+5]; 13 ll f[maxn+5]; 14 int n,m; 15 16 ll Pow(ll a,ll b,ll mod) 17 { 18 ll ans=1; 19 while(b) 20 { 21 if(b&1) ans=ans*a%mod; 22 a=a*a%mod; 23 b>>=1; 24 } 25 return ans; 26 } 27 ll Inv(ll x) 28 { 29 return Pow(x,mod-2,mod); 30 } 31 void Init() 32 { 33 ll inv=Inv(G); 34 for (int i=1;i<=maxn;i<<=1) 35 { 36 pw[i]=Pow(G,(mod-1)/i,mod); 37 pwinv[i]=Pow(inv,(mod-1)/i,mod); 38 } 39 } 40 void rader(ll *a) 41 { 42 for(int i=0,j=0;i<len;i++) 43 { 44 if(i>j) swap(a[i],a[j]); 45 int k=len; 46 do{k>>=1;j^=k;}while(j<k); 47 } 48 } 49 void ntt(ll *a,int f) 50 { 51 rader(a); 52 for(int i=2;i<=len;i<<=1) 53 { 54 int m=i>>1; 55 for(int j=0;j<len;j+=i) 56 { 57 ll w=1,wn=pw[i]; 58 if(f==-1) wn=pwinv[i]; 59 for(int k=0;k<m;k++) 60 { 61 ll x=a[j+k+m]*w%mod; 62 a[j+k+m]=(a[j+k]-x+mod)%mod; 63 a[j+k]=(a[j+k]+x)%mod; 64 w=w*wn%mod; 65 } 66 } 67 } 68 if(f==-1) 69 { 70 ll inv=Inv(len); 71 for(int i=0;i<len;i++) a[i]=a[i]*inv%mod; 72 } 73 } 74 void con(ll *A,int n,ll *B,int m) 75 { 76 /*A[0..n-1]与B[0..m-1]卷积*/ 77 for(len=1;len<max(n,m);len<<=1); 78 len<<=1; 79 for(int i=n;i<len;++i) A[i]=0; 80 for(int i=m;i<len;++i) B[i]=0; 81 ntt(A,1); 82 ntt(B,1); 83 for(int i=0;i<len;++i) A[i]=A[i]*B[i]%mod; 84 ntt(A,-1); 85 } 86 int main() 87 { 88 Init(); 89 con(A,n,B,m); 90 return 0; 91 }
FWT
1 #include<bits/stdc++.h> 2 using namespace std; 3 const int maxn=1024,mod=1e9+7,rev=(mod+1)>>1; 4 int a[maxn+50],b[maxn+50]; 5 long long fc[maxn*maxn]; 6 int T,n,m; 7 void fwt(int *a,int n) 8 { 9 for(int d=1;d<n;d<<=1) 10 for(int m=d<<1,i=0;i<n;i+=m) 11 for(int j=0;j<d;j++) 12 { 13 int x=a[i+j],y=a[i+j+d]; 14 a[i+j]=(x+y)%mod,a[i+j+d]=(x-y+mod)%mod; 15 //xor:a[i+j]=x+y,a[i+j+d]=(x-y+mod)%mod; 16 //and:a[i+j]=x+y; 17 //or:a[i+j+d]=x+y; 18 } 19 } 20 void ufwt(int *a,int n) 21 { 22 for(int d=1;d<n;d<<=1) 23 for(int m=d<<1,i=0;i<n;i+=m) 24 for(int j=0;j<d;j++) 25 { 26 int x=a[i+j],y=a[i+j+d]; 27 a[i+j]=1LL*(x+y)*rev%mod,a[i+j+d]=(1LL*(x-y)*rev%mod+mod)%mod; 28 //xor:a[i+j]=(x+y)/2,a[i+j+d]=(x-y)/2; 29 //and:a[i+j]=x-y; 30 //or:a[i+j+d]=y-x; 31 } 32 } 33 void solve(int *a,int *b,int n)//下标0..n-1的数组a和b求异或卷积，O(nlogn)，返回值在a中 34 { 35 fwt(a,n); 36 fwt(b,n); 37 for(int i=0;i<n;++i) a[i]=1LL*a[i]*b[i]%mod; 38 ufwt(a,n); 39 } 40 int main() 41 { 42 fc[0]=1; 43 for(int i=1;i<=1000000;++i) fc[i]=fc[i-1]*i%mod; 44 scanf("%d",&T); 45 while(T--) 46 { 47 scanf("%d%d",&n,&m); 48 memset(a,0,sizeof(a)); 49 memset(b,0,sizeof(b)); 50 for(int i=1;i<=n;++i) a[i]=1; 51 for(int i=1;i<=m;++i) b[i]=1; 52 solve(a,b,maxn); 53 long long ans=1; 54 for(int i=0;i<maxn;++i) ans=ans*fc[a[i]]%mod; 55 printf("%lld ",ans); 56 } 57 return 0; 58 }
BM算法（模数只能是质数）
1 #include<bits/stdc++.h> 2 using namespace std; 3 #define rep(i,a,n) for (int i=a;i<n;i++) 4 #define per(i,a,n) for (int i=n-1;i>=a;i--) 5 #define pb push_back 6 #define mp make_pair 7 #define all(x) (x).begin(),(x).end() 8 #define fi first 9 #define se second 10 #define SZ(x) ((int)(x).size()) 11 typedef vector<int> VI; 12 typedef long long ll; 13 typedef pair<int,int> PII; 14 const ll mod=1000000007; 15 ll powmod(ll a,ll b) {ll res=1;a%=mod; assert(b>=0); for(;b;b>>=1){if(b&1)res=res*a%mod;a=a*a%mod;}return res;} 16 // head 17 18 ll n; 19 namespace linear_seq { 20 const int N=10010; 21 ll res[N],base[N],_c[N],_md[N]; 22 23 vector<int> Md; 24 void mul(ll *a,ll *b,int k) { 25 rep(i,0,k+k) _c[i]=0; 26 rep(i,0,k) if (a[i]) rep(j,0,k) _c[i+j]=(_c[i+j]+a[i]*b[j])%mod; 27 for (int i=k+k-1;i>=k;i--) if (_c[i]) 28 rep(j,0,SZ(Md)) _c[i-k+Md[j]]=(_c[i-k+Md[j]]-_c[i]*_md[Md[j]])%mod; 29 rep(i,0,k) a[i]=_c[i]; 30 } 31 int solve(ll n,VI a,VI b) { // a 系数 b 初值 b[n+1]=a[0]*b[n]+... 32 // printf("%d ",SZ(b)); 33 ll ans=0,pnt=0; 34 int k=SZ(a); 35 assert(SZ(a)==SZ(b)); 36 rep(i,0,k) _md[k-1-i]=-a[i];_md[k]=1; 37 Md.clear(); 38 rep(i,0,k) if (_md[i]!=0) Md.push_back(i); 39 rep(i,0,k) res[i]=base[i]=0; 40 res[0]=1; 41 while ((1ll<<pnt)<=n) pnt++; 42 for (int p=pnt;p>=0;p--) { 43 mul(res,res,k); 44 if ((n>>p)&1) { 45 for (int i=k-1;i>=0;i--) res[i+1]=res[i];res[0]=0; 46 rep(j,0,SZ(Md)) res[Md[j]]=(res[Md[j]]-res[k]*_md[Md[j]])%mod; 47 } 48 } 49 rep(i,0,k) ans=(ans+res[i]*b[i])%mod; 50 if (ans<0) ans+=mod; 51 return ans; 52 } 53 VI BM(VI s) { 54 VI C(1,1),B(1,1); 55 int L=0,m=1,b=1; 56 rep(n,0,SZ(s)) { 57 ll d=0; 58 rep(i,0,L+1) d=(d+(ll)C[i]*s[n-i])%mod; 59 if (d==0) ++m; 60 else if (2*L<=n) { 61 VI T=C; 62 ll c=mod-d*powmod(b,mod-2)%mod; 63 while (SZ(C)<SZ(B)+m) C.pb(0); 64 rep(i,0,SZ(B)) C[i+m]=(C[i+m]+c*B[i])%mod; 65 L=n+1-L; B=T; b=d; m=1; 66 } else { 67 ll c=mod-d*powmod(b,mod-2)%mod; 68 while (SZ(C)<SZ(B)+m) C.pb(0); 69 rep(i,0,SZ(B)) C[i+m]=(C[i+m]+c*B[i])%mod; 70 ++m; 71 } 72 } 73 return C; 74 } 75 int gao(VI a,ll n) { 76 VI c=BM(a); 77 c.erase(c.begin()); 78 rep(i,0,SZ(c)) c[i]=(mod-c[i])%mod; 79 return solve(n,c,VI(a.begin(),a.begin()+SZ(c))); 80 } 81 }; 82 vector<int> a; 83 int main() { 84 int T; 85 scanf("%d",&T); 86 a.clear(); 87 a.pb(31),a.pb(197),a.pb(1255),a.pb(7997),a.pb(50959),a.pb(324725); 88 while(T--) 89 { 90 scanf("%lld",&n); 91 printf("%d ",linear_seq::gao(a,n-1)); 92 } 93 return 0; 94 }
k^2logn求线性递推
1 #include<bits/stdc++.h> 2 using namespace std; 3 const int maxn=4000,mod=1000000007; 4 int a[maxn+5],p[maxn+5],ans[maxn+5],num[maxn+5]; 5 int h[maxn+5],tmp[maxn+5]; 6 int n,k; 7 void mul(int *a,int *b,int *ans) 8 { 9 for(int i=0;i<=2*k;++i) tmp[i]=0; 10 for(int i=0;i<k;++i) 11 for(int j=0;j<k;++j) 12 tmp[i+j]=(tmp[i+j]+1LL*a[i]*b[j])%mod; 13 for(int i=2*k-2;i>=k;--i) 14 { 15 for(int j=k-1;j>=0;--j) 16 tmp[i-k+j]=(tmp[i-k+j]-1LL*tmp[i]*p[j])%mod,tmp[i-k+j]=(tmp[i-k+j]+mod)%mod; 17 tmp[i]=0; 18 } 19 for(int i=0;i<k;++i) ans[i]=tmp[i]; 20 } 21 int main() 22 { 23 scanf("%d%d",&n,&k); 24 for(int i=1;i<=k;++i) scanf("%d",&a[i]); 25 for(int i=0;i<k;++i) scanf("%d",&h[i]); 26 p[k]=1; 27 for(int i=1;i<=k;++i) p[k-i]=mod-a[i]; 28 for(int i=k;i<2*k;++i) 29 for(int j=1;j<=k;++j) 30 { 31 h[i]=h[i]+1LL*h[i-j]*a[j]%mod; 32 h[i]%=mod; 33 } 34 if(n<2*k) return 0*printf("%d ",h[n]); 35 int b=n-k+1; 36 num[1]=1,ans[0]=1; 37 while(b) 38 { 39 if(b&1) mul(ans,num,ans); 40 mul(num,num,num); 41 b>>=1; 42 } 43 long long res=0; 44 for(int i=0;i<k;++i) res=(res+1LL*ans[i]*h[i+k-1])%mod; 45 printf("%lld ",(res+mod)%mod); 46 return 0; 47 }
lagrange插值
1 void mul(int *a,int d,int c) 2 { 3 /* 4 a*(x-c) 5 */ 6 int b[maxn+5]; 7 // memset(b,0,sizeof(b)); 8 for(int i=0;i<=d;++i) b[i]=a[i]; 9 for(int i=d+1;i>=1;--i) a[i]=a[i-1]; 10 a[0]=0; 11 for(int i=0;i<=d;++i) 12 { 13 a[i]=(a[i]-1LL*b[i]*c)%mod; 14 if(a[i]<0) a[i]+=mod; 15 } 16 } 17 void div(int *a,int d,int c,int *ans) 18 { 19 /* 20 ans=a/(x-c) 21 */ 22 ans[d-1]=a[d]; 23 for(int i=d-2;i>=0;--i) 24 ans[i]=(a[i+1]+1LL*c*ans[i+1])%mod; 25 } 26 void lagrange(int *y,int d,int *p) 27 { 28 /* 29 (0,y[0]),(1,y[1]),...,(d,y[d]) 30 插出的系数放p中 31 */ 32 int tmp[maxn+5],now[maxn+5]; 33 memset(tmp,0,sizeof(tmp)); 34 memset(now,0,sizeof(now)); 35 tmp[0]=1; 36 for(int i=0;i<=d;++i) 37 { 38 mul(tmp,i,i); 39 } 40 for(int i=0;i<=d;++i) 41 { 42 for(int j=0;j<=d;++j) now[j]=0; 43 div(tmp,d+1,i,now); 44 long long q=1; 45 for(int j=0;j<=d;++j) 46 if(i!=j) 47 { 48 q=q*(i-j)%mod; 49 if(q<0) q+=mod; 50 } 51 q=inv(q); 52 q=q*y[i]%mod; 53 for(int j=0;j<=d;++j) now[j]=1LL*now[j]*q%mod; 54 for(int j=0;j<=d;++j) p[j]=(p[j]+now[j])%mod; 55 } 56 }
高斯消元
1 /* 2 n个未知数，m个方程 3 */ 4 #include<bits/stdc++.h> 5 using namespace std; 6 const int maxn=1000; 7 const double eps=1e-6; 8 double a[maxn+5][maxn+5]; 9 int n,m; 10 int main() 11 { 12 scanf("%d%d",&n,&m); 13 for(int i=1;i<=m;++i) 14 for(int j=1;j<=n+1;++j) 15 scanf("%lf",&a[i][j]); 16 for(int i=1;i<=n;++i)//枚举列 17 { 18 if(fabs(a[i][i]-0.0)<=eps) 19 { 20 int j; 21 bool flag=0; 22 for(j=i+1;j<=m;++j) 23 if(fabs(a[j][i]-0.0)>eps) 24 { 25 flag=1; 26 break; 27 } 28 if(!flag) return 0*printf("Many solutions");//如果第i~m行的第i列都是0，那么多解 29 for(int k=1;k<=n+1;++k) swap(a[i][k],a[j][k]); 30 } 31 for(int j=1;j<=n+1;++j) if(i!=j)a[i][j]/=a[i][i];a[i][i]=1.0; 32 for(int j=1;j<=m;++j) 33 if(j!=i) 34 { 35 for(int k=1;k<=n+1;++k) if(k!=i) a[j][k]-=a[i][k]*a[j][i];a[j][i]=0.0; 36 } 37 } 38 for(int i=n+1;i<=m;++i) 39 if(fabs(a[i][n+1]-0.0)>eps) return 0*printf("No solutions");//n+1~m这些方程应该所有变量都被消干净了，所以如果不为0，则说明无解 40 for(int i=1;i<=n;++i) printf("%d ",(int)(a[i][n+1]+0.5)); 41 return 0; 42 }
求行列式
1 const int maxn=1000; 2 ll a[maxn+5][maxn+5]; 3 int turn,n; 4 void gcd(ll a,ll b,ll &d,ll &x,ll &y) 5 { 6 if(!b) d=a,x=1,y=0; 7 else 8 { 9 ++turn; 10 gcd(b,a%b,d,y,x); 11 y-=x*(a/b); 12 } 13 } 14 ll det(ll n) 15 { 16 //求行列式a[0..n-1][0..n-1] 17 ll tmp1[maxn+5],tmp2[maxn+5]; 18 ll ans=1; 19 for(int i=0;i<n;++i) 20 { 21 for(int j=i+1;j<n;++j) 22 if(a[j][i]!=0) 23 { 24 ll A=a[i][i],B=a[j][i],d,x,y; 25 turn=0; 26 gcd(A,B,d,x,y); 27 for(int k=0;k<n;++k) tmp1[k]=a[i][k],tmp2[k]=a[j][k]; 28 for(int k=0;k<n;++k) a[i][k]=(x*tmp1[k]+y*tmp2[k])%mod; 29 A/=d,B/=d; 30 if(turn&1) x=B,y=-A,ans=-ans%mod;else x=-B,y=A; 31 for(int k=0;k<n;++k) a[j][k]=(x*tmp1[k]+y*tmp2[k])%mod; 32 } 33 ans=ans*a[i][i]%mod; 34 } 35 if(ans<0) ans+=mod; 36 return ans; 37 }
计算几何
求凸包
1 #include<cstring> 2 #include<algorithm> 3 #include<cstdio> 4 using namespace std; 5 const int maxn=50000; 6 struct wjmzbmr 7 { 8 int x,y; 9 bool operator < (const wjmzbmr& a) const 10 { 11 return (x<a.x)||(x==a.x&&y<a.y); 12 } 13 }a[maxn+50]; 14 int n,s[maxn+50],q[maxn+50],m=0,len=0; 15 int cross(int x1,int y1,int x2,int y2) 16 { 17 return (x1*y2-y1*x2); 18 } 19 int main() 20 { 21 scanf("%d",&n); 22 for(int i=1;i<=n;++i) scanf("%d %d",&a[i].x,&a[i].y); 23 sort(a+1,a+n+1); 24 m=2; 25 s[1]=1,s[2]=2; 26 for(int i=3;i<=n;++i) 27 { 28 while(m>=2&&cross(a[s[m]].x-a[s[m-1]].x,a[s[m]].y-a[s[m-1]].y,a[i].x-a[s[m]].x,a[i].y-a[s[m]].y)<=0) --m; 29 s[++m]=i; 30 } 31 s[++m]=n-1; 32 int k=m; 33 for(int i=n-2;i>=1;--i) 34 { 35 while(m>=k&&cross(a[s[m]].x-a[s[m-1]].x,a[s[m]].y-a[s[m-1]].y,a[i].x-a[s[m]].x,a[i].y-a[s[m]].y)<=0) --m; 36 s[++m]=i; 37 } 38 int ans=0; 39 for(int i=1;i<=m-1;++i) 40 for(int j=i+1;j<=m;++j) 41 if((a[s[i]].x-a[s[j]].x)*(a[s[i]].x-a[s[j]].x)+(a[s[i]].y-a[s[j]].y)*(a[s[i]].y-a[s[j]].y)>ans) ans=(a[s[i]].x-a[s[j]].x)*(a[s[i]].x-a[s[j]].x)+(a[s[i]].y-a[s[j]].y)*(a[s[i]].y-a[s[j]].y); 42 printf("%d",ans); 43 return 0; 44 }
自适应辛普森
1 #include<bits/stdc++.h> 2 using namespace std; 3 const double eps=1e-6; 4 double a,b,l,r; 5 double F(double x) 6 { 7 //Simpson公式用到的函数 8 return 2*sqrt(b*b-b*b*x*x/a/a); 9 } 10 double simpson(double l,double r)//三点Simpson法，这里要求F是一个全局函数 11 { 12 double mid=(l+r)*0.5; 13 return (F(l)+4*F(mid)+F(r))*(r-l)/6; 14 } 15 double asr(double l,double r)//自适应Simpson公式（递归过程） 16 { 17 double mid=(l+r)*0.5; 18 double x=simpson(l,mid),y=simpson(mid,r); 19 double z=simpson(l,r); 20 if(fabs(x+y-z)<=15*eps) return x+y+(x+y-z)/15.0; 21 return asr(l,mid)+asr(mid,r); 22 } 23 int main() 24 { 25 int T; 26 scanf("%d",&T); 27 while(T--) 28 { 29 scanf("%lf%lf%lf%lf",&a,&b,&l,&r); 30 printf("%.3f ",asr(l,r)); 31 } 32 return 0; 33 }
chelly的计算几何模板
1 typedef double db; 2 const db eps=1e-12; 3 struct Point 4 { 5 /*点类*/ 6 db x,y; 7 Point(){} 8 Point(db _x,db _y):x(_x),y(_y){} 9 Point operator + (const Point &t)const 10 { 11 return Point(x+t.x,y+t.y); 12 } 13 Point operator - (const Point &t)const 14 { 15 return Point(x-t.x,y-t.y); 16 } 17 Point operator * (const db &t)const 18 { 19 return Point(x*t,y*t); 20 } 21 Point operator / (const db &t)const 22 { 23 return Point(x/t,y/t); 24 } 25 db operator * (const Point &t)const 26 { 27 return x*t.y-y*t.x; 28 } 29 db len()const 30 { 31 return sqrtl(x*x+y*y); 32 } 33 Point rot90()const 34 { 35 return Point(-y,x); 36 } 37 }; 38 struct Complex 39 { 40 /*复数类，用于旋转*/ 41 db x,y; 42 Complex(db x=0.0,db y=0.0):x(x),y(y){} 43 Complex operator + (const Complex &t)const 44 { 45 return Complex(x+t.x,y+t.y); 46 } 47 Complex operator - (const Complex &t)const 48 { 49 return Complex(x-t.x,y-t.y); 50 } 51 Complex operator * (const Complex &t)const 52 { 53 return Complex(x*t.x-y*t.y,x*t.y+y*t.x); 54 } 55 Complex operator / (const Complex &t)const 56 { 57 return Complex((x*t.x+y*t.y)/(t.x*t.x+t.y*t.y),(y*t.x-x*t.y)/(t.x*t.x+t.y*t.y)); 58 } 59 }; 60 struct Circle 61 { 62 /*圆类*/ 63 Point o; 64 db r; 65 Circle(){} 66 Circle(Point _o,db _r):o(_o),r(_r){} 67 friend vector<Point> operator & (const Circle &c1,const Circle &c2) 68 { 69 /* 70 若两圆相离，返回空 71 若两圆相切，返回两个相同的点 72 若两圆相交，返回两个不同的点 73 */ 74 db d=(c1.o-c2.o).len(); 75 if(d>c1.r+c2.r+eps||d<fabs(c1.r-c2.r)-eps) return vector<Point>(); 76 db dt=(c1.r*c1.r-c2.r*c2.r)/d,d1=(d+dt)/2; 77 Point dir=(c2.o-c1.o)/d,pcrs=c1.o+dir*d1; 78 dt=sqrt(max(0.0L,c1.r*c1.r-d1*d1)),dir=dir.rot90(); 79 return vector<Point>{pcrs+dir*dt,pcrs-dir*dt}; 80 } 81 }; 82 void calcentre() 83 { 84 /* 85 计算凸多边形的重心 86 将凸多边形划分成若干小三角形，求出三角形的重心，然后根据小三角形面积求加权平均 87 */ 88 db area=0.0; 89 for(int i=2;i<n;++i) area+=(a[i+1]-a[i])*(a[1]-a[i]); 90 g=Point(0,0); 91 for(int i=2;i<n;++i) 92 { 93 db s=(a[i+1]-a[i])*(a[1]-a[i]); 94 g=g+(a[1]+a[i]+a[i+1])/3.0*(s/area); 95 } 96 }
ABerror计算几何
namespace Geometry {
using namespace std;
const double INF = 1e100;
const double EPS = 1e-12;
const double PI = acos(-1.0);
typedef complex<double> Pt;
typedef vector<Pt> Poly;
Pt IJ = Pt(0, 1);
int dcmp(double db) {
if (db < -EPS)return -1;
else return db > EPS;
}
int dcmp(double A, double B) {
return dcmp(A - B);
}
double cross(const Pt &A, const Pt &B) { return imag(conj(A) * B); }
double dot(const Pt &A, const Pt &B) { return real(conj(A) * B); }
double angle(const Pt &A) { return atan2(A.imag(), A.real()); }
Pt rotate(const Pt &A, double angle) { return A * exp(Pt(0, angle)); }
typedef struct Seg {
Pt A, B;
Seg(
double x1 = 0,
double y1 = 0,
double x2 = 0,
double y2 = 0
) :A(x1, y1), B(x2, y2) {
;
}
Seg(Pt A, Pt B) :A(A), B(B) { ; }
Pt vec()const { return B - A; }
bool operator<(const Seg & other)const {
return dcmp(angle(vec()) - angle(other.vec())) < 0;
}
}Line;
///////////////////////////////////////////////////////////////
//函数 :bool SegSegInterset(Seg fst, Seg sec, Pt &P)
//功能 :计算前两条直线的交点并返回是否有交点
// 如果需要计算交点，就在P处返回
//注意 :需要提前处理共线情况
///////////////////////////////////////////////////////////////
bool SegSegInterset(Seg fst, Seg sec) {
assert(dcmp(cross(fst.vec(), sec.vec())) != 0);
double p1 = cross(fst.A - sec.A, sec.vec());
double p2 = cross(fst.B - sec.A, sec.vec());
double p3 = cross(sec.A - fst.A, fst.vec());
double p4 = cross(sec.B - fst.A, fst.vec());
return dcmp(p1 * p2) <= 0 && dcmp(p3 * p4) <= 0;
}
bool SegSegInterset(Seg fst, Seg sec, Pt &P) {
if (SegSegInterset(fst, sec)) {
Pt A = fst.A, B = sec.A;
Pt a = fst.B - fst.A, b = sec.B - sec.A;
double t1 = cross(A - B, b) / cross(b, a);
P = A + a * t1;
return true;
}
else return false;
}
bool LineLineInterset(Line fst, Line sec, Pt& P) {
if (dcmp(cross(fst.vec(), sec.vec())) != 0) {
Pt A = fst.A, B = sec.A;
Pt a = fst.B - fst.A, b = sec.B - sec.A;
double t1 = cross(A - B, b) / cross(b, a);
P = A + a * t1;
return true;
}
else
exit(0xf0f00000);
}
///////////////////////////////////////////////////////////////
//函数 : Poly& convex(vector<Pt> &pts)
//功能 : 计算pts里的凸包并返回（逆时针排序)
//注意 : 返回值保存在pts里面
// 这个函数假设了一定存在凸包，对于点全部共线或者重合未定义
///////////////////////////////////////////////////////////////
Pt cmpPt; Seg cmpSeg;
bool cmp(Pt A, Pt B) {
A = A - cmpPt; B = B - cmpPt;
if (dcmp(cross(A, B)) != 0)
return dcmp(cross(A, B)) > 0;
else
return dcmp(abs(A) - abs(B)) < 0;
}
Poly& convex(vector<Pt> &pts) {
int pos = 0;
rep(i, 0, (int)pts.size()) {
Pt &e = pts[i];
if (dcmp(e.real() - pts[pos].real()) < 0 ||
dcmp(e.real() - pts[pos].real()) == 0 &&
dcmp(e.imag() - pts[pos].imag()) < 0
) {
pos = i;
}
}
vector<Pt> ans;
ans.push_back(pts[pos]);
cmpPt = pts[pos];
pts.erase(pts.begin() + pos);
sort(pts.begin(), pts.end(), cmp);
rep(j, 0, pts.size()) {
Pt &i = pts[j];
while (ans.size() >= 2u &&
cross(
ans[ans.size() - 1] - ans[ans.size() - 2],
i - ans[ans.size() - 2]
) <= 0) {
ans.pop_back();
}
ans.push_back(i);
}
return pts = ans;
}
///////////////////////////////////////////////////
//函数： Poly HalfPlaneIntersect(vector<Seg> seg)
//功能： 半平面交
//注意： 1. 这个程序不能判断交集为一个点的情况
// 2. 如果真的遇见这种情况，可以把所有边右移 1~2 个EPS
///////////////////////////////////////////////////
inline bool OnLeft(const Line & line, Pt pt) {
return dcmp(cross(line.vec(), pt - line.A)) > 0;
}
Poly HalfPlaneIntersect(vector<Seg> seg) {
vector<Seg> s(seg.size() + 5);
vector<Pt> p(seg.size() + 5);
int N = (int)seg.size(), l = 0, r = 0;
sort(seg.begin(), seg.end());
s[r++] = seg[0];
rep(i, 1, N) {
//tips:这里会弹出在直线上方的点，因此不能处理最后的解为一个点的情况
while (r - l >= 2 && !OnLeft(seg[i], p[r - 1]))r--;
while (r - l >= 2 && !OnLeft(seg[i], p[l + 1]))l++;
s[r++] = seg[i];
if (dcmp(cross(s[r - 1].vec(), s[r - 2].vec())) == 0) {
//exit(0xf0f00011);//debug
//如果最后一条线在倒数第二条线的左边，则挑选这条线
if (OnLeft(s[r - 2], s[r - 1].A))
s[r - 2] = s[r - 1];
r--;
}
//如果符合了上述情况，下面这个语句会发生重新求交
//如果不符合，这个语句会求新的交点
//注意由于上面那个语句，使得半平面的数量少于2，因此要作判断
if (r - l >= 2)
LineLineInterset(s[r - 2], s[r - 1], p[r - 1]);
}
//注意，下面这个部分对于判断交集是否为空有重大作用
//如果要起到这个作用，交点在s[l]上的时候也应该被弹出
while (r - l >= 2 && !OnLeft(s[l], p[r - 1]))r--;
if (r - l <= 2)return Poly();//无交集返回空集
LineLineInterset(s[r - 1], s[l], p[l]);
return Poly(p.begin() + l, p.begin() + r);
}
///////////////////////////////////////////////////
//函数： bool PtOnSeg(Pt p, Seg seg)
//功能： 判断点在线段上
//注意： 在该程序中，在线段的两个端点上也算在线段上
///////////////////////////////////////////////////
bool PtOnSeg(Pt p, Seg seg) {
//判断共直线
Pt l = seg.vec();
assert(dcmp(abs(l)) != 0);
if (dcmp(cross(l, p - seg.A)) != 0)return false;
double d = dot(p - seg.A, l);
return dcmp(d) >= 0 && dcmp(d - dot(l, l)) <= 0;
}
///////////////////////////////////////////////////
//函数： int PtInPolygon(const Pt P, const Poly &poly)
//功能： 询问点是否在多边形内部
///////////////////////////////////////////////////
int PtInPolygon(const Pt P, const Poly &poly) {
int N = (int)poly.size();
double sum = 0;
rep(i, 0, N) {
//判断在交点上
if (dcmp(abs(P - poly[i])) == 0)
return 1;
//判断在线段上
if (PtOnSeg(P, Seg(poly[i], poly[(i + 1) % N])))
return 2;
Pt A = poly[i] - P;
Pt B = poly[(i + 1) % N] - P;
sum += angle(rotate(B, -angle(A)));
}
//判断在多边形内部,无论的逆时针还是顺时针都可以进行判断
return 3 * (dcmp(sum - 2 * PI) == 0 || dcmp(sum + 2 * PI) == 0);
}
///////////////////////////////////////////////////
//函数：double PtLineDistance(Pt P, Line L)
//功能： 计算点到直线上的距离
///////////////////////////////////////////////////
double PtLineDistance(Pt P, Line L) {
Pt v = L.vec();
return fabs(cross(P - L.A, v)) / abs(v);
}
///////////////////////////////////////////////////
//函数： Pt PtLineProjection(Pt P, Line L)
//功能： 计算点在直线上的投影
///////////////////////////////////////////////////
Pt PtLineProjection(Pt P, Line L) {
Pt v = L.vec();
return L.A + v * dot(P - L.A, v) / dot(v, v);
}
///////////////////////////////////////////////////
//函数: double PtSegDistance(Pt P, Seg s)
//功能: 计算点到线段的距离
///////////////////////////////////////////////////
double PtSegDistance(Pt P, Seg s) {
double AB = abs(s.vec());
double d = dot(P - s.A, s.vec() / AB); //AP在于AB方向单位向量点积
if (dcmp(d) >= 0 && dcmp(d - AB) <= 0) {//等于AP在AB方向投影长度
return PtLineDistance(P, s); //如果这个长度在0到|AB|之间
}
else return min(abs(P - s.A), abs(P - s.B));
}
/*************************************************************/
/////////////////////////////////////////////////////
//圆的部分（圆和圆弧)
//数据结构: Cir/Circle 圆
// Arc 圆弧
/////////////////////////////////////////////////////
typedef struct Cir {
Pt P; //圆心
double R; //半径
Cir(Pt P = Pt(), double R = 0) :P(P), R(R) { ; }
Pt pt(double angle) {
return P + R * exp(angle * IJ);
}
}Circle;
struct Arc :public Cir {
double l, r;//角度在[l,r]范围内
Arc(double l = 0, double r = 0, Cir C = Cir()) :
l(l), r(r), Cir(C) {
;
}
};
/////////////////////////////////////////////////////
//函数： LineCirIntersect(Line L, Cir C)
//功能：
// 直线和圆的交点
// 以Poly的形式返回直线和圆的交点的集合
// 有几个就返回几个
/////////////////////////////////////////////////////
Poly LineCirIntersect(Line L, Cir C) {
Poly ret;
double dis = PtLineDistance(C.P, L);
Pt Pjt = PtLineProjection(C.P, L);
Pt v = L.vec();
switch (dcmp(dis - C.R)) {
//相交
case -1:
v = v / abs(v) * sqrt(C.R * C.R - dis * dis);
ret.push_back(Pjt - v);
ret.push_back(Pjt + v);
//相切,就放进去投影点
case 0:ret.push_back(Pjt); break;
//相离，直接退出
case 1:break;
default:exit(0xf0f00002);
}
return ret;
}
///////////////////////////////////////////////////
//函数: bool PtOnArc(const Pt p,const Arc &arc)
//功能: 判断点在圆弧上
//注意： 1. 这个函数默认点在圆上，用于判定是否在圆弧上
// 2. 如果在圆弧端点上，认为在圆弧上
///////////////////////////////////////////////////
bool PtOnArc(const Pt p, const Arc &arc) {
double ang = angle(p - arc.P);
assert(dcmp(arc.l - arc.r) != 0);
if (arc.l < arc.r)
return dcmp(ang - arc.l) >= 0 && dcmp(ang - arc.r) <= 0;
else
return dcmp(ang - arc.l) >= 0 || dcmp(ang - arc.r) <= 0;
}
///////////////////////////////////////////////////
//函数: Poly SegArcIntersect(Seg S, Arc A)
//功能: 计算线段和圆弧的交点
//注意： 有几个交点就返回几个交点
///////////////////////////////////////////////////
Poly SegArcIntersect(Seg S, Arc A) {
Poly ret = LineCirIntersect(S, A);
for (int i = (int)ret.size() - 1; i >= 0; i--) {
if (!PtOnSeg(ret[i], S) || !PtOnArc(ret[i], A))
ret.erase(ret.begin() + i);
}
return ret;
}
///////////////////////////////////////////////////
//函数: Poly CirCirIntersect(Cir A, Cir B)
//功能: 两圆相交，计算交点，这里只是计算交点，不讨论类型
//注意： 某些情况下，要判断5种情况,设两圆半径r <= R,圆心距d
// 1.相离，R + r - d < 0
// 2.外切，R + r - d = 0
// 3.相交，R + r - d > 0 && R - r - d < 0
// 4.内切，R - r - d = 0
// 5.内含，R - r - d > 0
// 6.重合，R = r,d = 0
///////////////////////////////////////////////////
Poly CirCirIntersect(Cir A, Cir B) {
Poly ret;
if (A.R < B.R)swap(A, B);
Pt AB = B.P - A.P;
double R = A.R, r = B.R, d = abs(AB);
double cmp1 = dcmp(R + r - d), cmp2 = dcmp(R - r - d);
assert(!(dcmp(R - r) == 0 && dcmp(d) == 0));
if (cmp1 == 0 || cmp2 == 0) {
ret.push_back(A.P + AB * R / d);
}
else if (cmp1 > 0 && cmp2 < 0) {
double theta = angle(AB);
double alpha = acos((R * R + d * d - r * r) / (2 * R * d));
ret.push_back(A.P + R * exp((theta + alpha) * IJ));
ret.push_back(A.P + R * exp((theta - alpha) * IJ));
}
return ret;
}
///////////////////////////////////////////////////
//函数: Poly ArcArIntersect(Arc A, Arc B)
//功能: 计算两个圆弧的交点，有几个输出几个
//注意：
///////////////////////////////////////////////////
Poly ArcArcIntersect(Arc A, Arc B) {
Poly ret = CirCirIntersect(A, B);
for (int i = (int)ret.size() - 1; i >= 0; i--) {
if (!PtOnArc(ret[i], A) || !PtOnArc(ret[i], B))
ret.erase(ret.begin() + i);
}
return ret;
}
};
其它科技
手写bitset
1 struct Bitset 2 { 3 /*32位压成一个unsigned int，最多压8个*/ 4 unsigned int a[8]; 5 void clear() 6 { 7 for(int i=0;i<8;++i) a[i]=0; 8 } 9 void set(int x) 10 { 11 /*第x位置1*/ 12 a[x>>5]|=1U<<(x&31); 13 /* 14 x>>5即x/32，确定x在哪个数组里 15 x&31即x%32，确定x在对应数组里的第几小位 16 */ 17 } 18 void flip(int x) 19 { 20 /*将第x位反转*/ 21 a[x>>5]^=1U<<(x&31); 22 } 23 int get(int x) 24 { 25 /*返回第x位的值0/1*/ 26 return a[x>>5]&(1U<<(x&31)); 27 } 28 }num; 29 //手写bitset的强大之处就是可以快速取出所有是1的位置，不过其实bitset里也封装了_Find_first()和_Find_next(x) 30 for(int i=0;i<8;++i) 31 { 32 while(true) 33 { 34 if(!num.a[i]) break; 35 int p=__builtin_ctz(num.a[i]);//__builtin_ctz(unsigned int x) 是通过O(1)的时间返回数字x最右边连续0的个数 36 printf("%d ",i<<5|p); 37 num.flip(i<<5|p); 38 } 39 } 40 //手写bitset还可以取出连续的一段区间，但是封装的却不行 41 void make(int l,int r) 42 { 43 int shift=l&31; 44 int y=(r-l)>>5; 45 int j=l>>5; 46 for(int i=0;i<y;++i) 47 { 48 u num=(a.num[j]>>shift); 49 if(shift!=0) num|=a.num[j+1]<<(32-shift);//注意x<<32并不是0，而是x本身 50 ans.num[i]^=num; 51 ++j; 52 } 53 for(int i=l+32*y;i<=r;++i) 54 if(a.get(i)) 55 ans.flip(i-l); 56 }
决策单调性优化
1 #include<bits/stdc++.h> 2 using namespace std; 3 const int maxn=200,inf=1e7; 4 int a[maxn+50],f[maxn+50][50+5]; 5 int kk[maxn+50][50+5]; 6 int m,s,n; 7 struct wjmzbmr 8 { 9 int l,r,p; 10 }q[maxn+50]; 11 int cal(int j,int i,int p) 12 { 13 return f[p][j]+a[i]-a[p+1]+1; 14 } 15 int find(int j,int l,int r,int p,int pp) 16 { 17 int mid; 18 bool flag=0; 19 while(l<r) 20 { 21 mid=(l+r)>>1; 22 if(cal(j-1,mid,p)>cal(j-1,mid,pp)) l=mid+1;else r=mid; 23 } 24 if(cal(j-1,r,p)<cal(j-1,r,pp)) return r;else return r+1; 25 } 26 int main() 27 { 28 scanf("%d%d%d",&m,&s,&n); 29 for(int i=1;i<=n;++i) scanf("%d",&a[i]); 30 sort(a+1,a+n+1); 31 for(int i=0;i<=n;++i) 32 for(int j=0;j<=m;++j) 33 f[i][j]=inf; 34 memset(f[0],0,sizeof(f[0])); 35 for(int i=1;i<=n;++i) f[i][1]=a[i]-a[1]+1; 36 for(int j=2;j<=m;++j) 37 { 38 int head=1,tail=1; 39 q[1]={1,n,0}; 40 for(int i=1;i<=n;++i) 41 { 42 while(head<tail&&q[head].r<i) ++head; 43 f[i][j]=cal(j-1,i,q[head].p); 44 while(head<tail&&cal(j-1,q[tail].l,i)<cal(j-1,q[tail].l,q[tail].p)) --tail; 45 int position=find(j,q[tail].l,q[tail].r,i,q[tail].p); 46 if(position<=n) 47 { 48 q[tail+1]={position,n,i}; 49 q[tail].r=position-1; 50 if(q[tail].l>q[tail].r) ++head; 51 ++tail; 52 } 53 } 54 } 55 printf("%d ",f[n][m]); 56 return 0; 57 }
数位DP
1 //数位dp的每个数字都在它的前缀处统计，所以每次dp的数字其实是一个前缀 2 #include<bits/stdc++.h> 3 using namespace std; 4 #define mp make_pair 5 typedef long long ll; 6 const ll mod=998244353; 7 pair<ll,ll> dp[19][1024]; //first记录满足条件的数字个数，second记录满足条件的数字之和 8 int a[19]; 9 int d; 10 ll pw[20]; 11 long long n; 12 void update(ll &a,ll b) 13 { 14 a=(a+b)%mod; 15 } 16 pair<ll,ll> dfs(int len,int ok,bool flag,bool zero) //ok表示之前的位置已经用了哪些数字，flag表示是否有取值限制，zero表示是否有前导0 17 { 18 if(len==-1) return mp(1,0); 19 if(!flag&&!zero&&dp[len][ok].first!=-1) return dp[len][ok]; 20 pair<ll,ll> ans=mp(0,0); 21 int limit=flag?a[len]:9; 22 for(int i=0;i<=limit;++i) 23 { 24 int mask=ok|(1<<i); 25 if(i==0&&zero) mask=ok; 26 if(__builtin_popcount(mask)>d) continue; 27 pair<ll,ll> to=dfs(len-1,mask,flag&&i==a[len],zero&&i==0); 28 update(ans.first,to.first); 29 update(ans.second,to.second); 30 update(ans.second,1LL*to.first*pw[len]%mod*i%mod); 31 } 32 if(!flag) dp[len][ok]=ans; 33 return ans; 34 } 35 ll cal(ll n) 36 { 37 int len=-1; 38 while(n) 39 { 40 a[++len]=n%10; 41 n/=10; 42 } 43 for(int i=0;i<19;++i) 44 for(int j=0;j<1024;++j) 45 dp[i][j]=mp(-1,-1); 46 return dfs(len,0,1,1).second; 47 } 48 int main() 49 { 50 pw[0]=1; 51 for(int i=1;i<20;++i) pw[i]=pw[i-1]*10%mod; 52 ll l,r; 53 cin>>l>>r>>d; 54 printf("%lld ",(cal(r)-cal(l-1)+mod)%mod); 55 return 0; 56 }
FastIO（不支持读负数、EOF）
1 typedef long long LL; 2 namespace fastIO { 3 #define BUF_SIZE 100000 4 //fread -> read 5 bool IOerror = 0; 6 inline char nc() { 7 static char buf[BUF_SIZE], *p1 = buf + BUF_SIZE, *pend = buf + BUF_SIZE; 8 if(p1 == pend) { 9 p1 = buf; 10 pend = buf + fread(buf, 1, BUF_SIZE, stdin); 11 if(pend == p1) { 12 IOerror = 1; 13 return -1; 14 } 15 } 16 return *p1++; 17 } 18 inline bool blank(char ch) { 19 return ch == ' ' || ch == ' ' || ch == ' ' || ch == ' '; 20 } 21 inline void read(int &x) { 22 char ch; 23 while(blank(ch = nc())); 24 if(IOerror) 25 return; 26 for(x = ch - '0'; (ch = nc()) >= '0' && ch <= '9'; x = x * 10 + ch - '0'); 27 } 28 #undef BUF_SIZE 29 }; 30 using namespace fastIO;
__int128读入输出
1 void scan(__int128 &x)//输入 2 { 3 x = 0; 4 int f = 1; 5 char ch; 6 if((ch = getchar()) == ‘-‘) f = -f; 7 else x = x*10 + ch-‘0‘; 8 while((ch = getchar()) >= ‘0‘ && ch <= ‘9‘) 9 x = x*10 + ch-‘0‘; 10 x *= f; 11 } 12 void print(__int128 x)//输出 13 { 14 if(x < 0) 15 { 16 x = -x; 17 putchar(‘-‘); 18 } 19 if(x > 9) print(x/10); 20 putchar(x%10 + ‘0‘); 21 }
单纯形法
1 #include<bits/stdc++.h> 2 using namespace std; 3 const int maxn=500,maxm=500; 4 const double eps=1e-8; 5 int n,m,op,tot,q[maxn+5],idx[maxn+5],idy[maxm+5]; 6 double a[maxm+5][maxn+5]; 7 double x[maxn+5]; 8 void pivot(int x,int y) 9 { 10 swap(idy[x],idx[y]); 11 double tmp = a[x][y]; a[x][y] = 1/a[x][y]; 12 for (int i = 0;i <= n;++i) 13 if (y != i) a[x][i] /= tmp; 14 tot = 0; 15 for (int i = 0;i <= n;++i) 16 if (i != y&&(a[x][i] > eps||a[x][i] < -eps)) q[++tot] = i; 17 for (int i = 0;i <= m;++i) 18 { 19 if ((x == i)||(a[i][y] < eps&&a[i][y] > -eps)) continue; 20 for (int j = 1;j <= tot;++j) a[i][q[j]] -= a[x][q[j]]*a[i][y]; 21 a[i][y] = -a[i][y]/tmp; 22 } 23 } 24 int simplex(double &ans,double *x) 25 { 26 /* 27 -1:无界 28 0:无解 29 1:有最优解 30 有最优解时候答案放ans中，解放x[]中 31 */ 32 for(int i = 1;i <= n;++i) idx[i] = i; 33 for(int i = 1;i <= m;++i) idy[i] = i+n; 34 while(true) 35 { 36 int x = 0,y = 0; 37 for (int i = 1;i <= m;++i) 38 if (a[i][0] < -eps&&((!x)||(rand()&1))) x = i; 39 if (!x) break; 40 for (int i = 1;i <= n;++i) 41 if (a[x][i] < -eps&&((!y)||(rand()&1))) y = i; 42 if (!y) return 0; 43 pivot(x,y); 44 } 45 while (true) 46 { 47 int x = 0,y = 0; double mn = 1e15; 48 for (int i = 1;i <= n;++i) 49 if (a[0][i] > eps) { y = i; break; } 50 if (!y) break; 51 for (int i = 1;i <= m;++i) 52 if (a[i][y] > eps && a[i][0]/a[i][y] < mn) mn = a[i][0]/a[i][y],x = i; 53 if (!x) return -1; 54 pivot(x,y); 55 } 56 ans=-a[0][0]; 57 for (int i = 1;i <= m;++i) 58 if (idy[i] <= n) x[idy[i]] = a[i][0]; 59 return 1; 60 } 61 int main() 62 { 63 scanf("%d%d%d",&n,&m,&op); 64 srand(233); 65 for (int i = 1;i <= n;++i) scanf("%lf",&a[0][i]); //第0行放目标函数的系数 66 for (int i = 1;i <= m;++i) 67 { 68 for (int j = 1;j <= n;++j) scanf("%lf",&a[i][j]); 69 scanf("%lf",&a[i][0]); //第0列放每个限制右边的上界 70 } 71 double ans; 72 int tmp=simplex(ans,x); 73 if(tmp==0) puts("Infeasible"); 74 else 75 if(tmp==-1) puts("Unbounded"); 76 else 77 { 78 printf("%.9f ",ans); 79 if(op==1) for(int i=1;i<=n;++i) printf("%.9f ",x[i]); 80 } 81 return 0; 82 }
BigInt
1 typedef long long LL; 2 struct BigInt 3 { 4 static const int BASE = 100000000; //基数 5 static const int WIDTH = 32; //基宽 6 vector<int> s; //存储 7 8 //构造函数 -- 使用LL 9 BigInt(long long num = 0) 10 { 11 *this = num; 12 } 13 14 //赋值运算符 -- 使用LL 15 BigInt operator = (long long num) 16 { 17 s.clear(); 18 do 19 { 20 s.push_back(num % BASE); 21 num /= BASE; 22 } 23 while(num > 0); 24 return *this; 25 } 26 //赋值运算符 -- 使用string 27 BigInt operator = (const string& str) 28 { 29 s.clear(); 30 int x, len = (str.length() - 1) / WIDTH + 1; 31 for(int i = 0; i < len; i++) 32 { 33 int end = str.length() - i*WIDTH; 34 int start = max(0, end - WIDTH); 35 sscanf(str.substr(start, end-start).c_str(), "%d", &x); 36 s.push_back(x); 37 } 38 return *this; 39 } 40 41 //比较运算符 42 bool operator < (const BigInt& b) const 43 { 44 if(s.size() != b.s.size()) return s.size() < b.s.size(); 45 for(int i = s.size()-1; i >= 0; i--) 46 if(s[i] != b.s[i]) return s[i] < b.s[i]; 47 return false; //相等 48 } 49 bool operator > (const BigInt& b) const 50 { 51 return b < *this; 52 } 53 bool operator <= (const BigInt& b) const 54 { 55 return !(b < *this); 56 } 57 bool operator >= (const BigInt& b) const 58 { 59 return !(*this < b); 60 } 61 bool operator != (const BigInt& b) const 62 { 63 return b < *this || *this < b; 64 } 65 bool operator == (const BigInt& b) const 66 { 67 return !(b < *this) && !(*this < b); 68 } 69 70 //重载输入输出 71 friend ostream& operator << (ostream &out, const BigInt& x) 72 { 73 out << x.s.back(); 74 for(int i = x.s.size()-2; i >= 0; i--) 75 { 76 char buf[20]; 77 sprintf(buf, "%08d", x.s[i]); 78 for(int j = 0; j < strlen(buf); j++) out << buf[j]; 79 } 80 return out; 81 } 82 friend istream& operator >> (istream &in, BigInt& x) 83 { 84 string s; 85 if(!(in >> s)) return in; 86 x = s; 87 return in; 88 } 89 90 //加法 91 BigInt operator + (const BigInt& b) const 92 { 93 BigInt c; 94 c.s.clear(); 95 for(int i = 0, g = 0; ; i++) 96 { 97 if(g == 0 && i >= s.size() && i >= b.s.size()) break; 98 int x = g; 99 if(i < s.size()) x += s[i]; 100 if(i < b.s.size()) x += b.s[i]; 101 c.s.push_back(x % BASE); 102 g = x / BASE; 103 } 104 return c; 105 } 106 107 //加等于 108 BigInt operator += (const BigInt& b) 109 { 110 *this = *this + b; 111 return *this; 112 } 113 114 115 116 //减法 -- 大减小 -- 其他的在外部处理 117 BigInt operator - (const BigInt& b) const 118 { 119 BigInt c; 120 c.s.clear(); 121 if((*this) == b) return c = 0; 122 for(int i = 0, g = 0; ; ++i) 123 { 124 if(g == 0 && i >= s.size() && i >= b.s.size()) break; 125 int x = -g; 126 g = 0; 127 if(i < s.size()) x += s[i]; 128 if(i < b.s.size()) x -= b.s[i]; 129 if(x < 0) x += BASE, ++g; 130 c.s.push_back(x); 131 } 132 for(int i = c.s.size()-1; i >= 0 && c.s[i] == 0; --i) c.s.pop_back(); 133 return c; 134 } 135 136 137 //乘法 -- 未使用FFT 138 BigInt operator * (const BigInt& b) const 139 { 140 int s1 = s.size(), s2 = b.s.size(); 141 vector<LL> v(s1+s2+1,0); 142 BigInt c; 143 c.s.resize(s1+s2, 0); 144 for(int i = 0; i < s1; ++i) //相乘 145 for(int j = 0; j < s2; ++j) 146 { 147 v[i+j] += (LL)s[i]*b.s[j]; 148 } 149 for(int i = 0; i < s1+s2; ++i) //进位 150 { 151 v[i+1] += v[i]/BASE; 152 v[i] %= BASE; 153 c.s[i] = v[i]; 154 } 155 for(int i = c.s.size()-1; i >= 0 && c.s[i] == 0; --i) c.s.pop_back(); 156 return c; 157 } 158 159 //除法 -- 二分 160 BigInt operator / (const BigInt& b) const 161 { 162 int s1 = s.size(), s2 = b.s.size(); 163 int len = s1-s2; 164 BigInt c, x = *this, y = b; 165 if(len < 0) return c = 0; 166 167 c.s.resize(len+1, 0); 168 for(int i = len; i >= 0; --i) 169 { 170 //先将除数对齐 171 y.s.resize(s2+i,0); 172 for(int j = s2+i-1; j >= 0; --j) 173 { 174 if(j >= i) y.s[j] = b.s[j-i]; 175 else y.s[j] = 0; 176 } 177 //再二分找商 178 int L = 0, R = BASE; 179 BigInt t; 180 while(L < R) 181 { 182 int mid = (L+R)>>1; 183 t = mid; 184 if(t*y > x) R = mid; 185 else L = mid+1; 186 } 187 c.s[i] = L-1; 188 //更新被除数 189 t = L-1; 190 x = x - t*y; 191 } 192 for(int i = c.s.size()-1; i >= 0 && c.s[i] == 0; --i) c.s.pop_back(); 193 return c; 194 } 195 196 //取模 197 BigInt operator % (const BigInt& b) const 198 { 199 return (*this) - (*this)/b * b; 200 } 201 202 //开方 -- 二分 -- 浮点数可放大1e4精确一位 203 BigInt Sqrt() 204 { 205 BigInt L = 1, R = *this; 206 while(L < R) 207 { 208 BigInt mid = (L+R)/2; 209 if(mid*mid > *this) R = mid; 210 else L = mid+1; 211 } 212 return L-1; 213 } 214 };
Java大数类
1 import java.math.BigDecimal; 2 import java.math.BigInteger; 3 import java.util.Scanner; 4 5 Scanner cin=new Scanner(System.in); 6 7 BigInteger num1=new BigInteger("12345"); 8 BigInteger num2=cin.nextBigInteger(); 9 10 BigDecimal num3=new BigDecimal("123.45"); 11 BigDecimal num4=cin.nextBigDecimal(); 12 13 //BigInteger 14 public class Main { 15 public static void main(String[] args) { 16 BigInteger num1=new BigInteger("12345"); 17 BigInteger num2=new BigInteger("45"); 18 //加法 19 System.out.println(num1.add(num2)); 20 //减法 21 System.out.println(num1.subtract(num2)); 22 //乘法 23 System.out.println(num1.multiply(num2)); 24 //除法(相除取整) 25 System.out.println(num1.divide(num2)); 26 //取余 27 System.out.println(num1.mod(num2)); 28 //最大公约数GCD 29 System.out.println(num1.gcd(num2)); 30 //取绝对值 31 System.out.println(num1.abs()); 32 //取反 33 System.out.println(num1.negate()); 34 //取最大值 35 System.out.println(num1.max(num2)); 36 //取最小值 37 System.out.println(num1.min(num2)); 38 //是否相等 39 System.out.println(num1.equals(num2)); 40 } 41 } 42 43 //BigDecimal 44 public class Main { 45 public static void main(String[] args) { 46 BigDecimal num1=new BigDecimal("123.45"); 47 BigDecimal num2=new BigDecimal("4.5"); 48 //加法 49 System.out.println(num1.add(num2)); 50 //减法 51 System.out.println(num1.subtract(num2)); 52 //乘法 53 System.out.println(num1.multiply(num2)); 54 //除法（在divide的时候就设置好要精确的小数位数和舍入模式） 55 System.out.println(num1.divide(num2,10,BigDecimal.ROUND_HALF_DOWN)); 56 //取绝对值 57 System.out.println(num1.abs()); 58 //取反 59 System.out.println(num1.negate()); 60 //取最大值 61 System.out.println(num1.max(num2)); 62 //取最小值 63 System.out.println(num1.min(num2)); 64 //是否相等 65 System.out.println(num1.equals(num2)); 66 //判断大小( > 返回1, < 返回-1) 67 System.out.println(num2.compareTo(num1)); 68 } 69 }
Montgomery乗算
1 #include <bits/stdc++.h> 2 using namespace std; 3 #define rep(i,a,n) for (int i=a;i<n;i++) 4 // head 5 typedef unsigned long long u64; 6 typedef long long ll; 7 typedef __int128_t i128; 8 typedef __uint128_t u128; 9 int k; 10 u64 A0,A1,M0,M1,C,M; 11 12 struct Mod64 { 13 Mod64():n_(0) {} 14 Mod64(u64 n):n_(init(n)) {} 15 static u64 init(u64 w) { return reduce(u128(w) * r2); } 16 static void set_mod(u64 m) { 17 mod=m; assert(mod&1); 18 inv=m; rep(i,0,5) inv*=2-inv*m; 19 r2=-u128(m)%m; 20 } 21 static u64 reduce(u128 x) { 22 u64 y=u64(x>>64)-u64((u128(u64(x)*inv)*mod)>>64); 23 return ll(y)<0?y+mod:y; 24 } 25 Mod64& operator += (Mod64 rhs) { n_+=rhs.n_-mod; if (ll(n_)<0) n_+=mod; return *this; } 26 Mod64 operator + (Mod64 rhs) const { return Mod64(*this)+=rhs; } 27 Mod64& operator -= (Mod64 rhs) { n_-=rhs.n_; if (ll(n_)<0) n_+=mod; return *this; } 28 Mod64 operator - (Mod64 rhs) const { return Mod64(*this)-=rhs; } 29 Mod64& operator *= (Mod64 rhs) { n_=reduce(u128(n_)*rhs.n_); return *this; } 30 Mod64 operator * (Mod64 rhs) const { return Mod64(*this)*=rhs; } 31 u64 get() const { return reduce(n_); } 32 static u64 mod,inv,r2; 33 u64 n_; 34 }; 35 u64 Mod64::mod,Mod64::inv,Mod64::r2; 36 int main() { 37 int T; 38 scanf("%d",&T); 39 while(T--) 40 { 41 scanf("%llu%llu%llu%llu%llu%llu%d",&A0,&A1,&M0,&M1,&C,&M,&k); 42 Mod64::set_mod(M); 43 Mod64 a0(A0),a1(A1),m0(M0),m1(M1),c(C),ans(1),a2(0); 44 for (int i=0;i<=k;i++) { 45 ans=ans*a0; 46 a2=m0*a1+m1*a0+c; 47 a0=a1; a1=a2; 48 } 49 printf("%llu ",ans.get()); 50 } 51 }
一维模拟退火
1 /* 2 F(x)=6*x^7+8*x^6+7*x^3+5*x^2-y*x (0 <= x <=100) 3 求F(x)最小值 4 */ 5 #include<bits/stdc++.h> 6 using namespace std; 7 const double EPS=1e-6; 8 const double r=0.98; 9 const int dx[2]={-1,1}; 10 double y; 11 12 double F(double x) 13 { 14 return 6*pow(x,7)+8*pow(x,6)+7*pow(x,3)+5*pow(x,2)-y*x; 15 } 16 17 int main() 18 { 19 int T; 20 cin>>T; 21 while(T--) 22 { 23 cin>>y; 24 double step=100; 25 double x=0; 26 while(step>EPS) 27 { 28 for(int i=0;i<2;i++) 29 { 30 double next_x=x+dx[i]*step; 31 if(next_x<0||next_x>100) 32 continue; 33 if(F(next_x)<F(x)) 34 x=next_x; 35 } 36 step*=r; 37 } 38 printf("%.4f ",F(x)); 39 } 40 return 0; 41 }
二维模拟退火
1 /* 2 给你一个搜索范围 ( 0,0 ) 到 ( X ,Y ）这个矩形范围,给你N个点, 3 要求找出一个点使这个点到这 N 个点的最大距离最小。 4 可以想象的就是把一维转变成二维，那么单点搜索我们就可以变成随机多点并行搜索， 5 原来的双向移动，转变为360度随机移动。最后在多点的搜索结果中取最优。 6 显然移动是以圆的范围移动的，那么温度上界 T 取矩形对角线， 即圆的最大半径， 7 精度为小数点后1位， 精确度EPS取1e-3, 降温速率 r 取0.8，对于每个点我随机了36次移动 8 */ 9 #include<bits/stdc++.h> 10 using namespace std; 11 12 const double EPS=1e-3; 13 const double PI=acos(-1.0); 14 const double INF=0x3f3f3f3f; 15 const double r=0.8; 16 const int maxn=1e4+5; 17 const int num=20; 18 19 struct Point{ 20 double x,y; 21 }; 22 23 Point a[maxn]; 24 Point m[maxn]; 25 int n; 26 double d[maxn]; 27 28 double dist(Point A, Point B) 29 { 30 return sqrt((A.x-B.x)*(A.x-B.x)+(A.y-B.y)*(A.y-B.y)); 31 } 32 33 double MAX(Point t) 34 { 35 double res=-1; 36 for(int i=0;i<n;i++) 37 res=max(res,dist(t,a[i])); 38 return res; 39 } 40 41 int main() 42 { 43 double xx,yy; 44 srand(time(0)); 45 while(cin>>xx>>yy>>n) 46 { 47 for(int i=0;i<n;i++) 48 cin>>a[i].x>>a[i].y; 49 for(int i=0;i<num;i++) 50 { 51 m[i].x=(rand()%1000+1)/1000.0*xx; 52 m[i].y=(rand()%1000+1)/1000.0*yy; 53 d[i]=MAX(m[i]); 54 } 55 double T=sqrt(xx*xx+yy*yy); 56 Point next,now; 57 while(T>EPS) 58 { 59 for(int i=0;i<num; i++) 60 { 61 now.x=m[i].x; 62 now.y=m[i].y; 63 for(int j=0;j<36;j++) 64 { 65 double ang=(rand()%1000+1)/1000.0*2*PI; 66 next.x=now.x+cos(ang)*T; 67 next.y=now.y+sin(ang)*T; 68 if(next.x<0||next.x>xx||next.y<0||next.y>yy) 69 continue; 70 if(MAX(next)<d[i]) 71 { 72 m[i].x=next.x; 73 m[i].y=next.y; 74 d[i]=MAX(next); 75 } 76 } 77 } 78 T*=r; 79 } 80 int res; 81 double ans=INF; 82 for(int i=0;i<num;i++) 83 { 84 if(d[i]<ans) 85 { 86 ans=d[i]; 87 res=i; 88 } 89 } 90 printf("(%.1f,%.1f). ",m[res].x,m[res].y); 91 printf("%.1f ",ans); 92 } 93 return 0; 94 }
更新日志：
2018.9.19 创建版本1.0.0
2018.9.25 更新SAM模板拓扑排序部分（同时支持自动机和parent树的拓扑序）
2018.9.25 更新杜教筛写法，实现了O(1)记忆化（鸣谢Feynman1999）
2018.10.9 加入“最小费用最大流（dijkstrar实现）”，这个费用流板子跑得飞快
2018.10.10 更新单纯形模板
2018.10.12 加入“min25筛”模板
2018.10.19 加入“KM算法”模板
2018.10.23 更新了“类欧几里得”模板
2018.10.25 加入“Montgomery乗算”模板
2018.10.26 更新了“数位DP”模板
2018.11.23 加入了“一维模拟退火”模板
2018.11.23 加入了“二维模拟退火”模板