【快速傅里叶变换】
〖相关资料〗
《从多项式乘法到快速傅里叶变换》 《Fast Fourier Transform》
《【快速傅里叶变换】【FFT】【WikiOI】【P3132】【高精度练习之超大整数乘法】》
〖模板代码〗
[FFT]
1 const double pi=acos(-1); 2 struct cpx{double r,i;cpx(double _r=0,double _i=0):r(_r),i(_i){};}; 3 cpx operator + (cpx a,cpx b){return cpx(a.r+b.r,a.i+b.i);} 4 cpx operator - (cpx a,cpx b){return cpx(a.r-b.r,a.i-b.i);} 5 cpx operator * (cpx a,cpx b){return cpx(a.r*b.r-a.i*b.i,a.r*b.i+a.i*b.r);} 6 void fft(cpx *a,int n,int f) 7 { 8 int k=0;while((1<<k)<n)k++; 9 for(int i=0;i<n;i++) 10 { 11 int t=0; 12 for(int j=0;j<k;j++) 13 if(i&(1<<j))t|=(1<<(k-j-1)); 14 if(i<t)swap(a[i],a[t]); 15 } 16 for(int l=2;l<=n;l<<=1) 17 { 18 int m=l>>1; 19 cpx nw=cpx(cos(2*pi/l),sin(2*pi/l)*f); 20 for(cpx *p=a;p!=a+n;p+=l) 21 { 22 cpx w=cpx(1,0); 23 for(int i=0;i<m;i++) 24 { 25 cpx t=p[m+i]*w;w=w*nw; 26 p[m+i]=p[i]-t;p[i]=p[i]+t; 27 } 28 } 29 } 30 } 31 void mul(cpx *c1,int n1,cpx *c2,int n2,int *ans) 32 { 33 int n=1;while(n<n1+n2)n<<=1; 34 fft(c1,n,1);fft(c2,n,1); 35 for(int i=0;i<n;i++)c1[i]=c1[i]*c2[i]; 36 fft(c1,n,-1); 37 for(int i=0;i<n1+n2-1;i++)ans[i]=(int)(c1[i].r/n+0.5); 38 }
[MTT]
1 #include<cstdio> 2 #include<algorithm> 3 #include<cstring> 4 #include<cstdlib> 5 #include<cmath> 6 #define double long double 7 #define LL long long 8 using namespace std; 9 const int N=524288+5; 10 const double pi=acos(-1); 11 int n,n1,n2,mod,qmod=32768; 12 struct cpx{double r,i;cpx(double _r=0,double _i=0):r(_r),i(_i){};}; 13 cpx operator + (cpx a,cpx b){return cpx(a.r+b.r,a.i+b.i);} 14 cpx operator - (cpx a,cpx b){return cpx(a.r-b.r,a.i-b.i);} 15 cpx operator * (cpx a,cpx b){return cpx(a.r*b.r-a.i*b.i,a.r*b.i+a.i*b.r);} 16 cpx a1[N],a2[N],b1[N],b2[N]; 17 cpx c1[N],c2[N],c3[N]; 18 cpx w,omg[N],omgi[N]; 19 struct mtt{LL s[N];}a,b,ans,zero; 20 int read() 21 { 22 int x=0,f=1;char c=getchar(); 23 while(c<'0'||c>'9'){if(c=='-')f=-1;c=getchar();} 24 while(c>='0'&&c<='9'){x=x*10+c-'0';c=getchar();} 25 return x*f; 26 } 27 void fft(cpx *a,int n,int f) 28 { 29 int k=0;while((1<<k)<n)k++; 30 for(int i=0;i<n;i++) 31 { 32 int t=0; 33 for(int j=0;j<k;j++) 34 if(i&(1<<j))t|=(1<<(k-j-1)); 35 if(i<t)swap(a[i],a[t]); 36 } 37 for(int l=2;l<=n;l<<=1) 38 { 39 int m=l>>1; 40 for(cpx *p=a;p!=a+n;p+=l) 41 for(int i=0;i<m;i++) 42 { 43 if(f==1)w=omg[n/l*i]; 44 else w=omgi[n/l*i]; 45 cpx t=p[m+i]*w; 46 p[m+i]=p[i]-t;p[i]=p[i]+t; 47 } 48 } 49 if(f==-1)for(int i=0;i<n;i++)a[i].r/=n; 50 } 51 mtt operator * (mtt x,mtt y) 52 { 53 for(int i=0;i<n;i++) 54 { 55 a1[i].r=x.s[i]/qmod;a1[i].i=0; 56 a2[i].r=x.s[i]%qmod;a2[i].i=0; 57 b1[i].r=y.s[i]/qmod;b1[i].i=0; 58 b2[i].r=y.s[i]%qmod;b2[i].i=0; 59 } 60 fft(a1,n,1);fft(a2,n,1); 61 fft(b1,n,1);fft(b2,n,1); 62 for(int i=0;i<n;i++) 63 { 64 c1[i]=a1[i]*b1[i]; 65 c2[i]=a2[i]*b1[i]+a1[i]*b2[i]; 66 c3[i]=a2[i]*b2[i]; 67 } 68 fft(c1,n,-1);fft(c2,n,-1);fft(c3,n,-1); 69 mtt ret=zero; 70 for(int i=0;i<n;i++) 71 ret.s[i]=((LL)(c1[i].r+0.5)%mod*qmod%mod*qmod%mod 72 +(LL)(c2[i].r+0.5)%mod*qmod%mod 73 +(LL)(c3[i].r+0.5)%mod)%mod; 74 return ret; 75 } 76 int main() 77 { 78 n1=read()+1;n2=read()+1;mod=read(); 79 n=1;while(n<n1+n2)n<<=1; 80 for(int i=0;i<n;i++) 81 { 82 omg[i]=cpx(cos(2*pi/n*i),sin(2*pi/n*i)); 83 omgi[i]=cpx(cos(2*pi/n*i),sin(-2*pi/n*i)); 84 } 85 for(int i=0;i<n1;i++)a.s[i]=read(); 86 for(int i=0;i<n2;i++)b.s[i]=read(); 87 ans=a*b; 88 for(int i=0;i<n1+n2-1;i++)printf("%lld ",ans.s[i]); 89 return 0; 90 }
〖相关题目〗
1.【bzoj2179】FFT快速傅立叶
题意：给出两个n位10进制整数x和y，输出x*y。
分析：FFT裸题
1 #include<cstdio> 2 #include<algorithm> 3 #include<cmath> 4 #include<cstdlib> 5 #include<cstring> 6 #include<iostream> 7 using namespace std; 8 const double pi=acos(-1); 9 const int N=131072+5; 10 struct cpx{double r,i;cpx(double _r=0,double _i=0):r(_r),i(_i){};}a[N],b[N]; 11 cpx operator + (cpx a,cpx b){return cpx(a.r+b.r,a.i+b.i);} 12 cpx operator - (cpx a,cpx b){return cpx(a.r-b.r,a.i-b.i);} 13 cpx operator * (cpx a,cpx b){return cpx(a.r*b.r-a.i*b.i,a.r*b.i+a.i*b.r);} 14 int n,ans[N]; 15 char ch[N]; 16 void fft(cpx *a,int n,int f) 17 { 18 int k=0;while((1<<k)<n)k++; 19 for(int i=0;i<n;i++) 20 { 21 int t=0; 22 for(int j=0;j<k;j++) 23 if(i&(1<<j))t|=(1<<(k-j-1)); 24 if(i<t)swap(a[i],a[t]); 25 } 26 for(int l=2;l<=n;l<<=1) 27 { 28 int m=l>>1; 29 cpx nw=cpx(cos(2*pi/l),sin(2*pi/l)*f); 30 for(cpx *p=a;p!=a+n;p+=l) 31 { 32 cpx w=cpx(1,0); 33 for(int i=0;i<m;i++) 34 { 35 cpx t=p[m+i]*w;w=w*nw; 36 p[m+i]=p[i]-t;p[i]=p[i]+t; 37 } 38 } 39 } 40 } 41 void mul(cpx *c1,int n1,cpx *c2,int n2,int *ans) 42 { 43 int n=1;while(n<n1+n2)n<<=1; 44 fft(c1,n,1);fft(c2,n,1); 45 for(int i=0;i<n;i++)c1[i]=c1[i]*c2[i]; 46 fft(c1,n,-1); 47 for(int i=0;i<n1+n2-1;i++)ans[i]=(int)(c1[i].r/n+0.5); 48 } 49 int main() 50 { 51 scanf("%d",&n); 52 scanf("%s",ch); 53 for(int i=0;i<n;i++)a[i].r=ch[n-i-1]-'0'; 54 scanf("%s",ch); 55 for(int i=0;i<n;i++)b[i].r=ch[n-i-1]-'0'; 56 mul(a,n,b,n,ans);int m=n*2-2; 57 for(int i=0;i<=m;i++) 58 if(ans[i]>=10) 59 { 60 ans[i+1]+=ans[i]/10;ans[i]%=10; 61 if(i==m)m++; 62 } 63 for(int i=m;i>=0;i--)printf("%d",ans[i]); 64 return 0; 65 } 66
2.【bzoj2194】快速傅立叶之二
题意：计算C[k]=sigma(a[i]*b[i-k]) 其中 k < = i < n。
分析：Heret1cの博客
1 #include<cstdio> 2 #include<algorithm> 3 #include<cmath> 4 #include<cstdlib> 5 #include<cstring> 6 #include<iostream> 7 using namespace std; 8 const double pi=acos(-1); 9 const int N=262144+5; 10 struct cpx{double r,i;cpx(double _r=0,double _i=0):r(_r),i(_i){};}a[N],b[N]; 11 cpx operator + (cpx a,cpx b){return cpx(a.r+b.r,a.i+b.i);} 12 cpx operator - (cpx a,cpx b){return cpx(a.r-b.r,a.i-b.i);} 13 cpx operator * (cpx a,cpx b){return cpx(a.r*b.r-a.i*b.i,a.r*b.i+a.i*b.r);} 14 int n,ans[N]; 15 int read() 16 { 17 int x=0,f=1;char c=getchar(); 18 while(c<'0'||c>'9'){if(c=='-')f=-1;c=getchar();} 19 while(c>='0'&&c<='9'){x=x*10+c-'0';c=getchar();} 20 return x*f; 21 } 22 void fft(cpx *a,int n,int f) 23 { 24 int k=0;while((1<<k)<n)k++; 25 for(int i=0;i<n;i++) 26 { 27 int t=0; 28 for(int j=0;j<k;j++) 29 if(i&(1<<j))t|=(1<<(k-j-1)); 30 if(i<t)swap(a[i],a[t]); 31 } 32 for(int l=2;l<=n;l<<=1) 33 { 34 int m=l>>1; 35 cpx nw=cpx(cos(2*pi/l),sin(2*pi/l)*f); 36 for(cpx *p=a;p!=a+n;p+=l) 37 { 38 cpx w=cpx(1,0); 39 for(int i=0;i<m;i++) 40 { 41 cpx t=p[m+i]*w;w=w*nw; 42 p[m+i]=p[i]-t;p[i]=p[i]+t; 43 } 44 } 45 } 46 } 47 void mul(cpx *c1,int n1,cpx *c2,int n2,int *ans) 48 { 49 int n=1;while(n<n1+n2)n<<=1; 50 fft(c1,n,1);fft(c2,n,1); 51 for(int i=0;i<n;i++)c1[i]=c1[i]*c2[i]; 52 fft(c1,n,-1); 53 for(int i=0;i<n1+n2-1;i++)ans[i]=(int)(c1[i].r/n+0.5); 54 } 55 int main() 56 { 57 n=read(); 58 for(int i=0;i<n;i++)a[i]=read(),b[n-i-1]=read(); 59 mul(a,n,b,n,ans); 60 for(int i=n-1;i<n*2-1;i++)printf("%d ",ans[i]); 61 return 0; 62 }
3.【bzoj3527】[Zjoi2014]力
题意：见原题
分析：Mychaelの博客
1 #include<cstdio> 2 #include<algorithm> 3 #include<cmath> 4 #include<cstdlib> 5 #include<cstring> 6 #include<iostream> 7 using namespace std; 8 const double pi=acos(-1); 9 const int N=262144+5; 10 struct cpx{double r,i;cpx(double _r=0,double _i=0):r(_r),i(_i){};}f1[N],f2[N],g[N]; 11 cpx operator + (cpx a,cpx b){return cpx(a.r+b.r,a.i+b.i);} 12 cpx operator - (cpx a,cpx b){return cpx(a.r-b.r,a.i-b.i);} 13 cpx operator * (cpx a,cpx b){return cpx(a.r*b.r-a.i*b.i,a.r*b.i+a.i*b.r);} 14 int n;double x,ans1[N],ans2[N]; 15 void fft(cpx *a,int n,int f) 16 { 17 int k=0;while((1<<k)<n)k++; 18 for(int i=0;i<n;i++) 19 { 20 int t=0; 21 for(int j=0;j<k;j++) 22 if(i&(1<<j))t|=(1<<(k-j-1)); 23 if(i<t)swap(a[i],a[t]); 24 } 25 for(int l=2;l<=n;l<<=1) 26 { 27 int m=l>>1; 28 cpx nw=cpx(cos(2*pi/l),sin(2*pi/l)*f); 29 for(cpx *p=a;p!=a+n;p+=l) 30 { 31 cpx w=cpx(1,0); 32 for(int i=0;i<m;i++) 33 { 34 cpx t=p[m+i]*w;w=w*nw; 35 p[m+i]=p[i]-t;p[i]=p[i]+t; 36 } 37 } 38 } 39 } 40 int main() 41 { 42 scanf("%d",&n); 43 for(int i=0;i<n;i++) 44 { 45 scanf("%lf",&x); 46 f1[i].r=f2[n-i-1].r=x; 47 } 48 for(int i=1;i<n;i++)g[i].r=1.0/i/i; 49 int mx=1;while(mx<n*2)mx<<=1; 50 fft(f1,mx,1);fft(f2,mx,1);fft(g,mx,1); 51 for(int i=0;i<mx;i++)f1[i]=f1[i]*g[i]; 52 for(int i=0;i<mx;i++)f2[i]=f2[i]*g[i]; 53 fft(f1,mx,-1);fft(f2,mx,-1); 54 for(int i=0;i<n*2-1;i++)ans1[i]=f1[i].r/mx; 55 for(int i=0;i<n*2-1;i++)ans2[i]=f2[i].r/mx; 56 for(int i=0;i<n;i++) 57 printf("%.3lf ",ans1[i]-ans2[n-i-1]); 58 return 0; 59 }
4.【bzoj3160】万径人踪灭
题意：给定一个由'a'和'b'构成的字符串，求不连续回文子序列的个数。
分析：PoPoQQQの博客
1 #include<cstdio> 2 #include<algorithm> 3 #include<cstring> 4 #include<cmath> 5 #define LL long long 6 using namespace std; 7 const int N=262144; 8 const int M=1e6+5; 9 const int mod=1e9+7; 10 const double pi=acos(-1); 11 int len,ans; 12 int bin[M],tmp[N],p[N]; 13 char ch[M],s[N]; 14 struct cpx{double r,i;cpx(double _r=0,double _i=0):r(_r),i(_i){};}a[N],b[N]; 15 cpx operator + (cpx a,cpx b){return cpx(a.r+b.r,a.i+b.i);} 16 cpx operator - (cpx a,cpx b){return cpx(a.r-b.r,a.i-b.i);} 17 cpx operator * (cpx a,cpx b){return cpx(a.r*b.r-a.i*b.i,a.r*b.i+a.i*b.r);} 18 void fft(cpx *a,int n,int f) 19 { 20 int k=0;while((1<<k)<n)k++; 21 for(int i=0;i<n;i++) 22 { 23 int t=0; 24 for(int j=0;j<k;j++) 25 if(i&(1<<j))t|=(1<<(k-j-1)); 26 if(i<t)swap(a[i],a[t]); 27 } 28 for(int l=2;l<=n;l<<=1) 29 { 30 int m=l>>1; 31 cpx nw=cpx(cos(2*pi/l),sin(2*pi/l)*f); 32 for(cpx *p=a;p!=a+n;p+=l) 33 { 34 cpx w=cpx(1,0); 35 for(int i=0;i<m;i++) 36 { 37 cpx t=p[m+i]*w;w=w*nw; 38 p[m+i]=p[i]-t;p[i]=p[i]+t; 39 } 40 } 41 } 42 } 43 int manacher() 44 { 45 int n=0;s[++n]='#'; 46 for(int i=1;i<=len;i++)s[++n]=ch[i],s[++n]='#'; 47 int mxr=0,id=0,sum=0; 48 for(int i=1;i<=n;i++) 49 { 50 if(i<mxr)p[i]=min(p[2*id-i],mxr-i);else p[i]=1; 51 for(;i+p[i]<=n&&s[i-p[i]]==s[i+p[i]];p[i]++); 52 if(i+p[i]>mxr)mxr=i+p[i],id=i; 53 sum=(sum+p[i]/2)%mod; 54 } 55 return sum; 56 } 57 int main() 58 { 59 scanf("%s",ch+1);len=strlen(ch+1); 60 bin[0]=1;for(int i=1;i<=M-5;i++)bin[i]=bin[i-1]*2%mod; 61 int n=1;while(n<len*2)n<<=1; 62 for(int i=1;i<=len;i++)a[i].r=(ch[i]=='a'); 63 fft(a,n,1);for(int i=0;i<n;i++)b[i]=a[i]*a[i]; 64 for(int i=0;i<n;i++)a[i].r=a[i].i=0; 65 for(int i=1;i<=len;i++)a[i].r=(ch[i]=='b'); 66 fft(a,n,1);for(int i=0;i<n;i++)b[i]=b[i]+a[i]*a[i]; 67 fft(b,n,-1);for(int i=0;i<n;i++)tmp[i]=(LL)(b[i].r/n+0.5); 68 for(int i=0;i<n;i++)ans=(ans+bin[(tmp[i]+1)>>1]-1)%mod; 69 printf("%d",(ans-manacher()+mod)%mod); 70 return 0; 71 }
5.【bzoj3992】[SDOI2015]序列统计
题意：见原题
1 #include<cstdio> 2 #include<cstring> 3 #include<cmath> 4 #include<cstdlib> 5 #include<algorithm> 6 #define LL long long 7 using namespace std; 8 const double pi=acos(-1); 9 const int N=2e4+5; 10 const int mod=1004535809; 11 int n,ci,m,x,num,tmp,qmod; 12 int lg[N],b[N]; 13 int read() 14 { 15 int x=0,f=1;char c=getchar(); 16 while(c<'0'||c>'9'){if(c=='-')f=-1;c=getchar();} 17 while(c>='0'&&c<='9'){x=x*10+c-'0';c=getchar();} 18 return x*f; 19 } 20 struct cpx{double r,i;cpx(double _r=0,double _i=0):r(_r),i(_i){};} 21 w,a1[N],b1[N],a2[N],b2[N],c1[N],c2[N],c3[N],omg[N],omgi[N]; 22 cpx operator + (cpx a,cpx b){return cpx(a.r+b.r,a.i+b.i);} 23 cpx operator - (cpx a,cpx b){return cpx(a.r-b.r,a.i-b.i);} 24 cpx operator * (cpx a,cpx b){return cpx(a.r*b.r-a.i*b.i,a.r*b.i+a.i*b.r);} 25 void fft(cpx *a,int n,int f) 26 { 27 int k=0;while((1<<k)<n)k++; 28 for(int i=0;i<n;i++) 29 { 30 int t=0; 31 for(int j=0;j<k;j++) 32 if(i&(1<<j))t|=(1<<(k-j-1)); 33 if(i<t)swap(a[i],a[t]); 34 } 35 for(int l=2;l<=n;l<<=1) 36 { 37 int m=l>>1; 38 for(cpx *p=a;p!=a+n;p+=l) 39 { 40 for(int i=0;i<m;i++) 41 { 42 if(f==1)w=omg[n/l*i]; 43 else w=omgi[n/l*i]; 44 cpx t=p[m+i]*w; 45 p[m+i]=p[i]-t;p[i]=p[i]+t; 46 } 47 } 48 } 49 if(f==-1)for(int i=0;i<n;i++)a[i].r/=n; 50 } 51 struct mtt{LL s[N];}zero,st,ans; 52 int power(int a,int b,int mod) 53 { 54 int ans=1; 55 while(b) 56 { 57 if(b&1)ans=1ll*ans*a%mod; 58 a=1ll*a*a%mod;b>>=1; 59 } 60 return ans; 61 } 62 int find(int n) 63 { 64 int mx=(int)(sqrt(n)+0.5),num=n-1,tot=0; 65 for(int i=2;i<=mx;i++) 66 if(num%i==0) 67 { 68 b[++tot]=i; 69 while(num%i==0)num/=i; 70 } 71 if(num!=1)b[++tot]=num; 72 for(int i=2;i<=n;i++) 73 { 74 bool ok=true; 75 for(int j=1;j<=tot;j++) 76 if(power(i,(n-1)/b[j],n)==1){ok=false;break;} 77 if(ok)return i; 78 } 79 return 0; 80 } 81 void pre() 82 { 83 int g=find(m),sum=1; 84 for(int i=0;i<m-1;i++) 85 { 86 lg[sum]=i; 87 sum=1ll*sum*g%m; 88 } 89 } 90 int tim=0; 91 mtt operator * (mtt x,mtt y) 92 { 93 for(int i=0;i<n;i++) 94 { 95 a1[i].r=x.s[i]/qmod;a1[i].i=0; 96 a2[i].r=x.s[i]%qmod;a2[i].i=0; 97 b1[i].r=y.s[i]/qmod;b1[i].i=0; 98 b2[i].r=y.s[i]%qmod;b2[i].i=0; 99 } 100 fft(a1,n,1);fft(a2,n,1); 101 fft(b1,n,1);fft(b2,n,1); 102 // printf("-----------%d------- ",n); 103 tim++; 104 for(int i=0;i<n;i++) 105 { 106 c1[i]=a1[i]*b1[i]; 107 c2[i]=a2[i]*b1[i]+a1[i]*b2[i]; 108 c3[i]=a2[i]*b2[i]; 109 // if(tim<=10)printf("%.3lf %.3lf %.3lf ",c1[i].r,c2[i].r,c3[i].r); 110 } 111 fft(c1,n,-1);fft(c2,n,-1);fft(c3,n,-1); 112 mtt ret=zero; 113 for(int i=0;i<n;i++) 114 { 115 LL temp=((LL)(c1[i].r+0.5)%mod*qmod%mod*qmod%mod 116 +(LL)(c2[i].r+0.5)%mod*qmod%mod 117 +(LL)(c3[i].r+0.5)%mod)%mod; 118 ret.s[i%(m-1)]=(ret.s[i%(m-1)]+temp)%mod; 119 } 120 return ret; 121 } 122 void work(mtt a,int b) 123 { 124 bool flag=false; 125 while(b) 126 { 127 if(b&1) 128 { 129 if(!flag)ans=a,flag=true; 130 else ans=ans*a; 131 } 132 a=a*a;b>>=1; 133 // puts(""); 134 // for(int i=0;i<n;i++)printf("%lld ",a.s[i]); 135 // puts(""); 136 } 137 } 138 int main() 139 { 140 ci=read();m=read();x=read();num=read(); 141 pre();n=1;while(n<m+m)n<<=1; 142 for(int i=0;i<n;i++) 143 { 144 omg[i]=cpx(cos(2*pi/n*i),sin(2*pi/n*i)); 145 omgi[i]=cpx(cos(2*pi/n*i),sin(-2*pi/n*i)); 146 } 147 for(int i=0;i<n;i++)zero.s[i]=0; 148 st=zero; 149 for(int i=1;i<=num;i++) 150 { 151 tmp=read(); 152 if(tmp)st.s[lg[tmp]]=1; 153 } 154 // for(int i=0;i<n;i++)printf("%d %d ",i,lg[i]); 155 // for(int i=0;i<n;i++)printf("%lld ",st.s[i]); 156 // puts(""); 157 qmod=sqrt(1.0*mod); 158 work(st,ci); 159 printf("%lld",(ans.s[lg[x]]%mod+mod)%mod); 160 return 0; 161 }
6.【bzoj4827】[Hnoi2017]礼物
题意：见原题
分析：cdcqの博客
1 #include<cstdio> 2 #include<algorithm> 3 #include<cstring> 4 #include<cmath> 5 #define LL long long 6 using namespace std; 7 const int N=262144+5; 8 const double pi=acos(-1); 9 int n,m,anss,sum,x[N],y[N],ans[N]; 10 int read() 11 { 12 int x=0,f=1;char c=getchar(); 13 while(c<'0'||c>'9'){if(c=='-')f=-1;c=getchar();} 14 while(c>='0'&&c<='9'){x=x*10+c-'0';c=getchar();} 15 return x*f; 16 } 17 struct cpx{double r,i;cpx(double _r=0,double _i=0):r(_r),i(_i){};}a[N],b[N]; 18 cpx operator + (cpx a,cpx b){return cpx(a.r+b.r,a.i+b.i);} 19 cpx operator - (cpx a,cpx b){return cpx(a.r-b.r,a.i-b.i);} 20 cpx operator * (cpx a,cpx b){return cpx(a.r*b.r-a.i*b.i,a.r*b.i+a.i*b.r);} 21 void fft(cpx *a,int n,int f) 22 { 23 int k=0;while((1<<k)<n)k++; 24 for(int i=0;i<n;i++) 25 { 26 int t=0; 27 for(int j=0;j<k;j++) 28 if(i&(1<<j))t|=(1<<(k-j-1)); 29 if(i<t)swap(a[i],a[t]); 30 } 31 for(int l=2;l<=n;l<<=1) 32 { 33 int m=l>>1; 34 cpx nw=cpx(cos(2*pi/l),sin(2*pi/l)*f); 35 for(cpx *p=a;p!=a+n;p+=l) 36 { 37 cpx w=cpx(1,0); 38 for(int i=0;i<m;i++) 39 { 40 cpx t=p[m+i]*w;w=w*nw; 41 p[m+i]=p[i]-t;p[i]=p[i]+t; 42 } 43 } 44 } 45 } 46 void mul(cpx *c1,int n1,cpx *c2,int n2,int *ans) 47 { 48 int n=1;while(n<n1+n2)n<<=1; 49 fft(c1,n,1);fft(c2,n,1); 50 for(int i=0;i<n;i++)c1[i]=c1[i]*c2[i]; 51 fft(c1,n,-1); 52 for(int i=0;i<n1+n2-1;i++)ans[i]=(int)(c1[i].r/n+0.5); 53 } 54 int main() 55 { 56 n=read();m=read(); 57 for(int i=0;i<n;i++)x[i]=read(); 58 for(int i=0;i<n;i++)a[i].r=x[n-i-1]; 59 for(int i=0;i<n;i++)y[i]=read(),b[i].r=b[i+n].r=y[i],sum+=x[i]-y[i]; 60 mul(a,n,b,n*2-1,ans); 61 for(int i=n-1;i<=2*n-2;i++)anss=max(anss,ans[i]); 62 anss=-anss*2;int mx=0x3f3f3f3f; 63 for(int c=0;c<=m;c++)mx=min(mx,sum*c*2+n*c*c); 64 anss+=mx; 65 for(int i=0;i<n;i++)anss+=x[i]*x[i]+y[i]*y[i]; 66 printf("%d",anss); 67 return 0; 68 }
7.【bzoj4259】残缺的字符串
题意：见原题
分析：Clarisの博客
1 #include<cstdio> 2 #include<cstring> 3 #include<algorithm> 4 #include<cmath> 5 using namespace std; 6 const int N=524288+5; 7 const double pi=acos(-1); 8 int n,n1,n2,cnt,out[N],a[N],b[N]; 9 char s1[N],s2[N]; 10 struct cpx{double r,i;cpx(double _r=0,double _i=0):r(_r),i(_i){};}; 11 cpx A[N],B[N],ans[N]; 12 cpx operator + (cpx a,cpx b){return cpx(a.r+b.r,a.i+b.i);} 13 cpx operator - (cpx a,cpx b){return cpx(a.r-b.r,a.i-b.i);} 14 cpx operator * (cpx a,cpx b){return cpx(a.r*b.r-a.i*b.i,a.r*b.i+a.i*b.r);} 15 void fft(cpx *a,int n,int f) 16 { 17 int k=0;while((1<<k)<n)k++; 18 for(int i=0;i<n;i++) 19 { 20 int t=0; 21 for(int j=0;j<k;j++) 22 if(i&(1<<j))t|=(1<<(k-j-1)); 23 if(i<t)swap(a[i],a[t]); 24 } 25 for(int l=2;l<=n;l<<=1) 26 { 27 int m=l>>1; 28 cpx nw=cpx(cos(2*pi/l),sin(2*pi/l)*f); 29 for(cpx *p=a;p!=a+n;p+=l) 30 { 31 cpx w=cpx(1,0); 32 for(int i=0;i<m;i++) 33 { 34 cpx t=p[m+i]*w;w=w*nw; 35 p[m+i]=p[i]-t;p[i]=p[i]+t; 36 } 37 } 38 } 39 if(f==-1)for(int i=0;i<n;i++)a[i].r/=n; 40 } 41 int main() 42 { 43 scanf("%d%d%s%s",&n1,&n2,s1,s2); 44 n=1;while(n<n1||n<n2)n<<=1; 45 for(int i=0;i<n1;i++)if(s1[i]!='*')a[n1-i-1]=s1[i]-'a'+1; 46 for(int i=0;i<n2;i++)if(s2[i]!='*')b[i]=s2[i]-'a'+1; 47 for(int i=0;i<n;i++)A[i]=cpx(a[i]*a[i]*a[i],0); 48 for(int i=0;i<n;i++)B[i]=cpx(b[i],0); 49 fft(A,n,1);fft(B,n,1); 50 for(int i=0;i<n;i++)ans[i]=A[i]*B[i]; 51 for(int i=0;i<n;i++)A[i]=cpx(a[i]*a[i],0); 52 for(int i=0;i<n;i++)B[i]=cpx(b[i]*b[i],0); 53 fft(A,n,1);fft(B,n,1); 54 for(int i=0;i<n;i++)ans[i]=ans[i]-A[i]*B[i]*cpx(2,0); 55 for(int i=0;i<n;i++)A[i]=cpx(a[i],0); 56 for(int i=0;i<n;i++)B[i]=cpx(b[i]*b[i]*b[i],0); 57 fft(A,n,1);fft(B,n,1); 58 for(int i=0;i<n;i++)ans[i]=ans[i]+A[i]*B[i]; 59 fft(ans,n,-1); 60 for(int i=n1-1;i<n2;i++) 61 if(ans[i].r<0.5)out[++cnt]=i-n1+2; 62 printf("%d ",cnt); 63 for(int i=1;i<=cnt;i++)printf("%d ",out[i]); 64 return 0; 65 }
【生成函数】
〖相关题目〗
1.【bzoj3028】食物
题意：给出n种重量为1的食物以及它们的限制，求出食物个数总和恰好为n的方案数。
分析：maijingの博客
1 #include<cstdio> 2 #include<cstring> 3 #include<algorithm> 4 using namespace std; 5 const int mod=1e4+7; 6 int n,ans; 7 char ch[505]; 8 int main() 9 { 10 scanf("%s",ch+1); 11 n=strlen(ch+1); 12 for(int i=1;i<=n;i++) 13 ans=(ans*10+ch[i]-'0')%mod; 14 printf("%d ",ans*(ans+1)%mod*(ans+2)%mod*1668%mod); 15 return 0; 16 }
【快速数论变换】
〖模板代码〗
[NTT]
1 const int N=131072+5; 2 const int mod=998244353; 3 int qmod(int a,int b) 4 { 5 int ans=1; 6 while(b) 7 { 8 if(b&1)ans=1ll*ans*a%mod; 9 a=1ll*a*a%mod;b>>=1; 10 } 11 return ans; 12 } 13 void ntt(int *a,int n,int f) 14 { 15 int k=0;while((1<<k)<n)k++; 16 for(int i=0;i<n;i++) 17 { 18 int t=0; 19 for(int j=0;j<k;j++) 20 if(i&(1<<j))t|=(1<<(k-j-1)); 21 if(i<t)swap(a[i],a[t]); 22 } 23 for(int l=2;l<=n;l<<=1) 24 { 25 int m=l>>1,nw=qmod(3,(mod-1)/l); 26 if(f==-1)nw=qmod(nw,mod-2); 27 for(int *p=a;p!=a+n;p+=l) 28 { 29 int w=1; 30 for(int i=0;i<m;i++) 31 { 32 int t=1ll*p[m+i]*w%mod; 33 p[m+i]=(p[i]-t+mod)%mod; 34 p[i]=(p[i]+t)%mod; 35 w=1ll*w*nw%mod; 36 } 37 } 38 } 39 if(f==-1) 40 { 41 int inv=qmod(n,mod-2); 42 for(int i=0;i<n;i++)a[i]=1ll*a[i]*inv%mod; 43 } 44 } 45 int main() 46 { 47 n=1;while(n<len1+len2)n<<=1; 48 ntt(a,n,1);ntt(b,n,1); 49 for(int i=0;i<n;i++)c[i]=1ll*a[i]*b[i]%mod; 50 ntt(c,n,-1); 51 return 0; 52 }
[多项式求逆]
1 int power(int a,int b) 2 { 3 int ans=1; 4 while(b) 5 { 6 if(b&1)ans=1ll*ans*a%mod; 7 a=1ll*a*a%mod;b>>=1; 8 } 9 return ans; 10 } 11 void ntt(int *a,int n,int f) 12 { 13 int k=0;while((1<<k)<n)k++; 14 for(int i=0;i<n;i++) 15 { 16 int t=0; 17 for(int j=0;j<k;j++) 18 if(i&(1<<j))t|=(1<<(k-j-1)); 19 if(i<t)swap(a[i],a[t]); 20 } 21 for(int l=2;l<=n;l<<=1) 22 { 23 int m=l>>1,nw=power(3,(mod-1)/l); 24 if(f==-1)nw=power(nw,mod-2); 25 for(int *p=a;p!=a+n;p+=l) 26 { 27 int w=1; 28 for(int i=0;i<m;i++) 29 { 30 int t=1ll*p[m+i]*w%mod; 31 p[m+i]=(p[i]-t+mod)%mod; 32 p[i]=(p[i]+t)%mod; 33 w=1ll*w*nw%mod; 34 } 35 } 36 } 37 if(f==-1) 38 { 39 int inv=power(n,mod-2); 40 for(int i=0;i<n;i++)a[i]=1ll*a[i]*inv%mod; 41 } 42 } 43 void pinv(int *a,int *b,int n) 44 { 45 if(n==1){b[0]=power(a[0],mod-2);return;} 46 pinv(a,b,n>>1);n<<=1; 47 for(int i=0;i<n/2;i++)c[i]=a[i],c[i+n/2]=0; 48 ntt(c,n,1);ntt(b,n,1); 49 for(int i=0;i<n;i++)b[i]=(2*b[i]%mod-1ll*c[i]*b[i]%mod*b[i]%mod+mod)%mod; 50 ntt(b,n,-1);for(int i=n/2;i<n;i++)b[i]=0; 51 }
[多项式开根]
1 int power(int a,int b) 2 { 3 int ans=1; 4 while(b) 5 { 6 if(b&1)ans=1ll*ans*a%mod; 7 a=1ll*a*a%mod;b>>=1; 8 } 9 return ans; 10 } 11 void ntt(int *a,int n,int f) 12 { 13 int k=0;while((1<<k)<n)k++; 14 for(int i=0;i<n;i++) 15 { 16 int t=0; 17 for(int j=0;j<k;j++) 18 if(i&(1<<j))t|=(1<<(k-j-1)); 19 if(i<t)swap(a[i],a[t]); 20 } 21 for(int l=2;l<=n;l<<=1) 22 { 23 int m=l>>1,nw=power(3,(mod-1)/l); 24 if(f==-1)nw=power(nw,mod-2); 25 for(int *p=a;p!=a+n;p+=l) 26 { 27 int w=1; 28 for(int i=0;i<m;i++) 29 { 30 int t=1ll*p[m+i]*w%mod; 31 p[m+i]=(p[i]-t+mod)%mod; 32 p[i]=(p[i]+t)%mod; 33 w=1ll*w*nw%mod; 34 } 35 } 36 } 37 if(f==-1) 38 { 39 int inv=power(n,mod-2); 40 for(int i=0;i<n;i++)a[i]=1ll*a[i]*inv%mod; 41 } 42 } 43 void pinv(int *a,int *b,int n) 44 { 45 if(n==1){b[0]=power(a[0],mod-2);return;} 46 pinv(a,b,n>>1);n<<=1; 47 for(int i=0;i<n/2;i++)c[i]=a[i],c[i+n/2]=0; 48 ntt(c,n,1);ntt(b,n,1); 49 for(int i=0;i<n;i++) 50 b[i]=(2*b[i]%mod-1ll*c[i]*b[i]%mod*b[i]%mod+mod)%mod; 51 ntt(b,n,-1);for(int i=n/2;i<n;i++)b[i]=0; 52 } 53 void psqrt(int *a,int *b,int *invb,int n) 54 { 55 if(n==1){b[0]=1;return;} 56 psqrt(a,b,invb,n>>1); 57 memset(invb,0,4*n);pinv(b,invb,n); 58 for(int i=0;i<n;i++)c[i]=a[i]; 59 for(int i=n;i<n+n;i++)c[i]=0; 60 n<<=1;ntt(c,n,1);ntt(b,n,1);ntt(invb,n,1); 61 for(int i=0;i<n;i++) 62 b[i]=(b[i]+1ll*invb[i]*c[i]%mod)%mod*inv%mod; 63 ntt(b,n,-1);for(int i=n/2;i<n;i++)b[i]=0; 64 }
〖相关题目〗
1.【bzoj2179】FFT快速傅立叶
题意：给出两个n位10进制整数x和y，输出x*y。
分析：NTT裸题
1 #include<cstdio> 2 #include<algorithm> 3 #include<cstring> 4 #define LL long long 5 using namespace std; 6 const int N=131072+5; 7 const int mod=998244353; 8 int n,len,a[N],b[N],c[N]; 9 char s[N]; 10 int qmod(int a,int b) 11 { 12 int ans=1; 13 while(b) 14 { 15 if(b&1)ans=1ll*ans*a%mod; 16 a=1ll*a*a%mod;b>>=1; 17 } 18 return ans; 19 } 20 void ntt(int *a,int n,int f) 21 { 22 int k=0;while((1<<k)<n)k++; 23 for(int i=0;i<n;i++) 24 { 25 int t=0; 26 for(int j=0;j<k;j++) 27 if(i&(1<<j))t|=(1<<(k-j-1)); 28 if(i<t)swap(a[i],a[t]); 29 } 30 for(int l=2;l<=n;l<<=1) 31 { 32 int m=l>>1,nw=qmod(3,(mod-1)/l); 33 if(f==-1)nw=qmod(nw,mod-2); 34 for(int *p=a;p!=a+n;p+=l) 35 { 36 int w=1; 37 for(int i=0;i<m;i++) 38 { 39 int t=1ll*p[m+i]*w%mod; 40 p[m+i]=(p[i]-t+mod)%mod; 41 p[i]=(p[i]+t)%mod; 42 w=1ll*w*nw%mod; 43 } 44 } 45 } 46 if(f==-1) 47 { 48 int inv=qmod(n,mod-2); 49 for(int i=0;i<n;i++)a[i]=1ll*a[i]*inv%mod; 50 } 51 } 52 int main() 53 { 54 scanf("%d",&len); 55 n=1;while(n<len*2)n<<=1; 56 scanf("%s",s+1); 57 for(int i=1;i<=len;i++)a[len-i]=s[i]-'0'; 58 scanf("%s",s+1); 59 for(int i=1;i<=len;i++)b[len-i]=s[i]-'0'; 60 ntt(a,n,1);ntt(b,n,1); 61 // for(int i=0;i<n;i++)printf("%d %d ",a[i],b[i]); 62 for(int i=0;i<n;i++)c[i]=1ll*a[i]*b[i]%mod; 63 ntt(c,n,-1); 64 for(int i=0;i<n;i++) 65 if(c[i]>=10) 66 { 67 c[i+1]+=c[i]/10;c[i]%=10; 68 if(i==n-1)n++; 69 } 70 n--;while(!c[n]&&n)n--; 71 for(int i=n;i>=0;i--)printf("%d",c[i]); 72 return 0; 73 }
2.【51nod1135】原根
题意：给出1个质数P，找出P最小的原根。
分析：无
1 #include<cstdio> 2 #include<cmath> 3 #include<cstring> 4 #include<algorithm> 5 using namespace std; 6 int mod,num,mx,tot,b[1005]; 7 int qmod(int a,int b) 8 { 9 int ans=1; 10 while(b) 11 { 12 if(b&1)ans=1ll*ans*a%mod; 13 a=1ll*a*a%mod;b>>=1; 14 } 15 return ans; 16 } 17 int main() 18 { 19 scanf("%d",&mod); 20 mx=(int)(sqrt(mod)+0.5);num=mod-1; 21 for(int i=2;i<=mx;i++) 22 if(num%i==0) 23 { 24 b[++tot]=i; 25 while(num%i==0)num/=i; 26 } 27 if(num!=1)b[++tot]=num; 28 for(int i=2;i<=mod;i++) 29 { 30 bool ok=true; 31 for(int j=1;j<=tot;j++) 32 if(qmod(i,(mod-1)/b[j])==1){ok=false;break;} 33 if(ok){printf("%d",i);return 0;} 34 } 35 return 0; 36 }
3.【bzoj3992】[SDOI2015]序列统计
题意：见原题
1 #include<cstdio> 2 #include<algorithm> 3 #include<cstring> 4 #include<cmath> 5 #include<cstdlib> 6 using namespace std; 7 const int N=2e4+5; 8 const int mod=1004535809; 9 int n,m,x,num,tmp; 10 int lg[N],b[N],c[N],ans[N],f[N]; 11 int read() 12 { 13 int x=0,f=1;char c=getchar(); 14 while(c<'0'||c>'9'){if(c=='-')f=-1;c=getchar();} 15 while(c>='0'&&c<='9'){x=x*10+c-'0';c=getchar();} 16 return x*f; 17 } 18 int power(int a,int b,int mod) 19 { 20 int ans=1; 21 while(b) 22 { 23 if(b&1)ans=1ll*ans*a%mod; 24 a=1ll*a*a%mod;b>>=1; 25 } 26 return ans; 27 } 28 int find(int n) 29 { 30 int mx=(int)(sqrt(n)+0.5),num=n-1,tot=0; 31 for(int i=2;i<=mx;i++) 32 if(num%i==0) 33 { 34 b[++tot]=i; 35 while(num%i==0)num/=i; 36 } 37 if(num!=1)b[++tot]=num; 38 for(int i=2;i<=n;i++) 39 { 40 bool ok=true; 41 for(int j=1;j<=tot;j++) 42 if(power(i,(n-1)/b[j],n)==1){ok=false;break;} 43 if(ok)return i; 44 } 45 return 0; 46 } 47 void pre() 48 { 49 int g=find(m),sum=1; 50 for(int i=0;i<m-1;i++) 51 { 52 lg[sum]=i; 53 sum=1ll*sum*g%m; 54 } 55 } 56 void ntt(int *a,int n,int f) 57 { 58 int k=0;while((1<<k)<n)k++; 59 for(int i=0;i<n;i++) 60 { 61 int t=0; 62 for(int j=0;j<k;j++) 63 if(i&(1<<j))t|=(1<<(k-j-1)); 64 if(i<t)swap(a[i],a[t]); 65 } 66 for(int l=2;l<=n;l<<=1) 67 { 68 int m=l>>1,nw=power(3,(mod-1)/l,mod); 69 if(f==-1)nw=power(nw,mod-2,mod); 70 for(int *p=a;p!=a+n;p+=l) 71 { 72 int w=1; 73 for(int i=0;i<m;i++) 74 { 75 int t=1ll*p[m+i]*w%mod; 76 p[m+i]=(p[i]-t+mod)%mod; 77 p[i]=(p[i]+t)%mod; 78 w=1ll*w*nw%mod; 79 } 80 } 81 } 82 if(f==-1) 83 { 84 int inv=power(n,mod-2,mod); 85 for(int i=0;i<n;i++)a[i]=1ll*a[i]*inv%mod; 86 } 87 } 88 void mul(int *a,int *b,int n) 89 { 90 for(int i=0;i<n;i++)c[i]=b[i]; 91 ntt(a,n,1);ntt(c,n,1); 92 for(int i=0;i<n;i++)a[i]=1ll*a[i]*c[i]%mod; 93 ntt(a,n,-1); 94 for(int i=m-1;i<n;i++) 95 a[i%(m-1)]=(a[i%(m-1)]+a[i])%mod,a[i]=0; 96 } 97 void work(int *a,int b) 98 { 99 int n=1;while(n<m+m)n<<=1; 100 ans[0]=1; 101 while(b) 102 { 103 if(b&1)mul(ans,a,n); 104 mul(a,a,n);b>>=1; 105 } 106 } 107 int main() 108 { 109 n=read();m=read();x=read();num=read(); 110 pre(); 111 for(int i=1;i<=num;i++) 112 { 113 tmp=read(); 114 if(tmp)f[lg[tmp]]=1; 115 } 116 work(f,n); 117 printf("%d",ans[lg[x]]); 118 return 0; 119 }
4.【bzoj4555】[Tjoi2016&Heoi2016]求和
题意：见原题
分析：zltttttの博客
1 #include<cstdio> 2 #include<algorithm> 3 #include<cstring> 4 #include<cmath> 5 #include<cstdlib> 6 using namespace std; 7 const int N=262144+5; 8 const int mod=998244353; 9 int n,m,ans; 10 int fac[N],inv[N],f[N],g[N],c[N]; 11 int power(int a,int b) 12 { 13 int ans=1; 14 while(b) 15 { 16 if(b&1)ans=1ll*ans*a%mod; 17 a=1ll*a*a%mod;b>>=1; 18 } 19 return ans; 20 } 21 void ntt(int *a,int n,int f) 22 { 23 int k=0;while((1<<k)<n)k++; 24 for(int i=0;i<n;i++) 25 { 26 int t=0; 27 for(int j=0;j<k;j++) 28 if(i&(1<<j))t|=(1<<(k-j-1)); 29 if(i<t)swap(a[i],a[t]); 30 } 31 for(int l=2;l<=n;l<<=1) 32 { 33 int m=l>>1,nw=power(3,(mod-1)/l); 34 if(f==-1)nw=power(nw,mod-2); 35 for(int *p=a;p!=a+n;p+=l) 36 { 37 int w=1; 38 for(int i=0;i<m;i++) 39 { 40 int t=1ll*p[m+i]*w%mod; 41 p[m+i]=(p[i]-t+mod)%mod; 42 p[i]=(p[i]+t)%mod; 43 w=1ll*w*nw%mod; 44 } 45 } 46 } 47 if(f==-1) 48 { 49 int inv=power(n,mod-2); 50 for(int i=0;i<n;i++)a[i]=1ll*a[i]*inv%mod; 51 } 52 } 53 void pinv(int *a,int *b,int n) 54 { 55 if(n==0){b[0]=power(a[0],mod-2);return;} 56 pinv(a,b,n>>1);n<<=1; 57 for(int i=0;i<n/2;i++)c[i]=a[i],c[i+n/2]=0; 58 ntt(c,n,1);ntt(b,n,1); 59 for(int i=0;i<n;i++)b[i]=(2*b[i]%mod-1ll*c[i]*b[i]%mod*b[i]%mod+mod)%mod; 60 ntt(b,n,-1);for(int i=n/2;i<n;i++)b[i]=0; 61 } 62 int main() 63 { 64 scanf("%d",&m); 65 n=1;while(n<=m)n<<=1; 66 fac[0]=1; 67 for(int i=1;i<=n;i++)fac[i]=1ll*fac[i-1]*i%mod; 68 inv[n]=power(fac[n],mod-2); 69 for(int i=n;i>=1;i--)inv[i-1]=1ll*inv[i]*i%mod; 70 for(int i=1;i<=n;i++)f[i]=1ll*(mod-inv[i])*2%mod; 71 f[0]=1;pinv(f,g,n); 72 for(int i=0;i<=m;i++)ans=(ans+1ll*g[i]*fac[i]%mod)%mod; 73 printf("%d",ans); 74 return 0; 75 }
5.【bzoj3456】城市规划
题意：求n个点的无向简单连通图个数。
分析：PoPoQQQの博客
1 #include<cstdio> 2 #include<cstring> 3 #include<algorithm> 4 #include<cmath> 5 using namespace std; 6 const int N=262144+5; 7 const int mod=1004535809; 8 int n,mx; 9 int fac[N],inv[N],c[N]; 10 int A[N],B[N],C[N]; 11 int power(int a,int b) 12 { 13 int ans=1; 14 while(b) 15 { 16 if(b&1)ans=1ll*ans*a%mod; 17 a=1ll*a*a%mod;b>>=1; 18 } 19 return ans; 20 } 21 void ntt(int *a,int n,int f) 22 { 23 int k=0;while((1<<k)<n)k++; 24 for(int i=0;i<n;i++) 25 { 26 int t=0; 27 for(int j=0;j<k;j++) 28 if(i&(1<<j))t|=(1<<(k-j-1)); 29 if(i<t)swap(a[i],a[t]); 30 } 31 for(int l=2;l<=n;l<<=1) 32 { 33 int m=l>>1,nw=power(3,(mod-1)/l); 34 if(f==-1)nw=power(nw,mod-2); 35 for(int *p=a;p!=a+n;p+=l) 36 { 37 int w=1; 38 for(int i=0;i<m;i++) 39 { 40 int t=1ll*p[m+i]*w%mod; 41 p[m+i]=(p[i]-t+mod)%mod; 42 p[i]=(p[i]+t)%mod; 43 w=1ll*w*nw%mod; 44 } 45 } 46 } 47 if(f==-1) 48 { 49 int inv=power(n,mod-2); 50 for(int i=0;i<n;i++)a[i]=1ll*a[i]*inv%mod; 51 } 52 } 53 void pinv(int *a,int *b,int n) 54 { 55 if(n==1){b[0]=power(a[0],mod-2);return;} 56 pinv(a,b,n>>1);n<<=1; 57 for(int i=0;i<n/2;i++)c[i]=a[i],c[i+n/2]=0; 58 ntt(c,n,1);ntt(b,n,1); 59 for(int i=0;i<n;i++) 60 b[i]=(2*b[i]%mod-1ll*c[i]*b[i]%mod*b[i]%mod+mod)%mod; 61 ntt(b,n,-1);for(int i=n/2;i<n;i++)b[i]=0; 62 } 63 int main() 64 { 65 scanf("%d",&n); 66 mx=1;while(mx<=n)mx<<=1; 67 fac[0]=1; 68 for(int i=1;i<=n;i++)fac[i]=1ll*fac[i-1]*i%mod; 69 inv[n]=power(fac[n],mod-2); 70 for(int i=n;i>=1;i--)inv[i-1]=1ll*inv[i]*i%mod; 71 for(int i=0;i<=n;i++) 72 C[i]=1ll*power(2,1ll*i*(i-1)/2%(mod-1))*inv[i]%mod; 73 pinv(C,A,mx); 74 // for(int i=1;i<=n;i++)printf("%d ",A[i]); 75 for(int i=1;i<=n;i++) 76 B[i]=1ll*power(2,1ll*i*(i-1)/2%(mod-1))*inv[i-1]%mod; 77 mx<<=1;ntt(A,mx,1);ntt(B,mx,1); 78 memset(C,0,sizeof(C)); 79 for(int i=0;i<mx;i++)C[i]=1ll*A[i]*B[i]%mod; 80 ntt(C,mx,-1); 81 printf("%lld",1ll*C[n]*fac[n-1]%mod); 82 return 0; 83 }
6.【bzoj3625】[Codeforces Round #250]小朋友和二叉树
题意：见原题
分析：wzq_QwQの博客
1 #include<cstdio> 2 #include<cstring> 3 #include<algorithm> 4 #include<cmath> 5 #define LL long long 6 using namespace std; 7 const int N=262144+5; 8 const int mod=998244353; 9 int n,num,m,x,inv; 10 int f[N],g[N],invg[N],c[N]; 11 int read() 12 { 13 int x=0,f=1;char c=getchar(); 14 while(c<'0'||c>'9'){if(c=='-')f=-1;c=getchar();} 15 while(c>='0'&&c<='9'){x=x*10+c-'0';c=getchar();} 16 return x*f; 17 } 18 int power(int a,int b) 19 { 20 int ans=1; 21 while(b) 22 { 23 if(b&1)ans=1ll*ans*a%mod; 24 a=1ll*a*a%mod;b>>=1; 25 } 26 return ans; 27 } 28 void ntt(int *a,int n,int f) 29 { 30 int k=0;while((1<<k)<n)k++; 31 for(int i=0;i<n;i++) 32 { 33 int t=0; 34 for(int j=0;j<k;j++) 35 if(i&(1<<j))t|=(1<<(k-j-1)); 36 if(i<t)swap(a[i],a[t]); 37 } 38 for(int l=2;l<=n;l<<=1) 39 { 40 int m=l>>1,nw=power(3,(mod-1)/l); 41 if(f==-1)nw=power(nw,mod-2); 42 for(int *p=a;p!=a+n;p+=l) 43 { 44 int w=1; 45 for(int i=0;i<m;i++) 46 { 47 int t=1ll*p[m+i]*w%mod; 48 p[m+i]=(p[i]-t+mod)%mod; 49 p[i]=(p[i]+t)%mod; 50 w=1ll*w*nw%mod; 51 } 52 } 53 } 54 if(f==-1) 55 { 56 int inv=power(n,mod-2); 57 for(int i=0;i<n;i++)a[i]=1ll*a[i]*inv%mod; 58 } 59 } 60 void pinv(int *a,int *b,int n) 61 { 62 if(n==1){b[0]=power(a[0],mod-2);return;} 63 pinv(a,b,n>>1);n<<=1; 64 for(int i=0;i<n/2;i++)c[i]=a[i],c[i+n/2]=0; 65 ntt(c,n,1);ntt(b,n,1); 66 for(int i=0;i<n;i++) 67 b[i]=(2*b[i]%mod-1ll*c[i]*b[i]%mod*b[i]%mod+mod)%mod; 68 ntt(b,n,-1);for(int i=n/2;i<n;i++)b[i]=0; 69 } 70 void psqrt(int *a,int *b,int *invb,int n) 71 { 72 if(n==1){b[0]=1;return;} 73 psqrt(a,b,invb,n>>1); 74 memset(invb,0,4*n);pinv(b,invb,n); 75 for(int i=0;i<n;i++)c[i]=a[i]; 76 for(int i=n;i<n+n;i++)c[i]=0; 77 n<<=1;ntt(c,n,1);ntt(b,n,1);ntt(invb,n,1); 78 for(int i=0;i<n;i++) 79 b[i]=(b[i]+1ll*invb[i]*c[i]%mod)%mod*inv%mod; 80 ntt(b,n,-1);for(int i=n/2;i<n;i++)b[i]=0; 81 } 82 int main() 83 { 84 num=read();m=read();inv=power(2,mod-2); 85 for(int i=1;i<=num;i++)x=read(),f[x]=1; 86 for(int i=1;i<=m;i++)f[i]=(-f[i]*4%mod+mod)%mod; 87 n=1;while(n<=m)n<<=1; 88 f[0]=1;psqrt(f,g,invg,n); 89 g[0]=(g[0]+1)%mod; 90 memset(f,0,sizeof(f)); 91 pinv(g,f,n); 92 for(int i=1;i<=m;i++)printf("%d ",f[i]*2%mod); 93 return 0; 94 }
【拉格朗日插值】
〖相关知识〗
给定 $n+1$ 个点 $(x_0,y_0),cdots (x_n,y_n)$ ，求穿过它们的 $n$ 次函数。
令 $ell _j(x)=prod _{i=0,i eq j}^{n}frac{x-x_i}{x_j-x_i}$ ，拉格朗日插值多项式 $L(x)=sum_{j=0}^{n}y_j ell _j(x)$。
〖相关题目〗
1.【bzoj4559】成绩比较
题意：见原题
分析：cdsszjjの博客
1 #include<cstdio> 2 #include<cstring> 3 #include<algorithm> 4 #include<cmath> 5 #define LL long long 6 using namespace std; 7 const int mod=1e9+7; 8 const int N=105; 9 int n,m,k; 10 int u[N],r[N],g[N],w[N]; 11 int c[N][N],f[N][N]; 12 int read() 13 { 14 int x=0,f=1;char c=getchar(); 15 while(c<'0'||c>'9'){if(c=='-')f=-1;c=getchar();} 16 while(c>='0'&&c<='9'){x=x*10+c-'0';c=getchar();} 17 return x*f; 18 } 19 int power(int a,int b) 20 { 21 int ans=1; 22 while(b) 23 { 24 if(b&1)ans=1ll*ans*a%mod; 25 a=1ll*a*a%mod;b>>=1; 26 } 27 return ans; 28 } 29 int lag(int u,int r) 30 { 31 memset(g,0,sizeof(g)); 32 for(int x=1;x<=n+1;x++) 33 { 34 for(int i=1;i<=x;i++) 35 g[x]=(g[x]+1ll*power(x-i,r-1)*power(i,n-r)%mod)%mod; 36 if(u==x)return g[x]; 37 } 38 for(int i=1;i<=n+1;i++) 39 { 40 w[i]=1; 41 for(int j=1;j<=n+1;j++) 42 if(i!=j)w[i]=(1ll*w[i]*(i-j)%mod+mod)%mod; 43 w[i]=power(w[i],mod-2); 44 } 45 int sum=0; 46 for(int i=1;i<=n+1;i++) 47 { 48 int tmp=1ll*w[i]*g[i]%mod; 49 for(int j=1;j<=n+1;j++) 50 if(i!=j)tmp=1ll*tmp*(u-j)%mod; 51 sum=(sum+tmp)%mod; 52 } 53 return sum; 54 } 55 int C(int n,int m) 56 { 57 if(n<0||m<0||n<m)return 0; 58 return c[n][m]; 59 } 60 int main() 61 { 62 n=read();m=read();k=read(); 63 for(int i=1;i<=m;i++)u[i]=read(); 64 for(int i=1;i<=m;i++)r[i]=read(); 65 for(int i=0;i<=n;i++)c[i][0]=1; 66 for(int i=1;i<=n;i++) 67 for(int j=1;j<=i;j++) 68 c[i][j]=(c[i-1][j]+c[i-1][j-1])%mod; 69 f[0][0]=1; 70 for(int i=1;i<=m;i++) 71 { 72 int tmp=lag(u[i],r[i]); 73 for(int j=0;j<=n-1;j++) 74 { 75 for(int w=0;w<=j;w++) 76 f[i][j]=(f[i][j]+1ll*C(n-1-w,j-w)*C(w,r[i]-1-j+w)%mod*f[i-1][w]%mod)%mod; 77 f[i][j]=1ll*f[i][j]*tmp%mod; 78 } 79 } 80 printf("%d",f[m][n-1-k]); 81 return 0; 82 }
2.【bzoj2655】calc
题意：见原题
分析：ez_ywwの博客
1 #include<cstdio> 2 #include<algorithm> 3 #include<cstring> 4 #include<cmath> 5 using namespace std; 6 const int N=505; 7 int m,n,mod,ans; 8 int w[N],f[N*2][N]; 9 int power(int a,int b) 10 { 11 int ans=1; 12 while(b) 13 { 14 if(b&1)ans=1ll*ans*a%mod; 15 a=1ll*a*a%mod;b>>=1; 16 } 17 return ans; 18 } 19 int main() 20 { 21 scanf("%d%d%d",&m,&n,&mod); 22 f[0][0]=1; 23 for(int i=1;i<=2*n;i++) 24 { 25 f[i][0]=f[i-1][0]; 26 for(int j=1;j<=n;j++) 27 f[i][j]=(1ll*f[i-1][j-1]*i%mod*j%mod+f[i-1][j])%mod; 28 } 29 if(m<=2*n){printf("%d",f[m][n]);return 0;} 30 for(int i=0;i<=2*n;i++) 31 { 32 w[i]=1; 33 for(int j=0;j<=2*n;j++) 34 if(i!=j)w[i]=(1ll*w[i]*(i-j)%mod+mod)%mod; 35 w[i]=power(w[i],mod-2); 36 } 37 for(int i=0;i<=2*n;i++) 38 { 39 int sum=1ll*f[i][n]*w[i]%mod; 40 for(int j=0;j<=2*n;j++) 41 if(i!=j)sum=1ll*sum*(m-j)%mod; 42 ans=(ans+sum)%mod; 43 } 44 printf("%d ",ans); 45 return 0; 46 }
3.【bzoj3453】tyvj 1858 XLkxc
题意：见原题
分析：DZYOの博客
1 #include<cstdio> 2 #include<cstring> 3 #include<algorithm> 4 #include<cmath> 5 #include<iostream> 6 using namespace std; 7 const int N=150; 8 const int mod=1234567891; 9 int T,k,a,n,d,m,sum,ans; 10 int g[N],w[N],f[N]; 11 int read() 12 { 13 int x=0,f=1;char c=getchar(); 14 while(c<'0'||c>'9'){if(c=='-')f=-1;c=getchar();} 15 while(c>='0'&&c<='9'){x=x*10+c-'0';c=getchar();} 16 return x*f; 17 } 18 int power(int a,int b) 19 { 20 int ans=1; 21 while(b) 22 { 23 if(b&1)ans=1ll*ans*a%mod; 24 a=1ll*a*a%mod;b>>=1; 25 } 26 return ans; 27 } 28 void work() 29 { 30 k=read();a=read();n=read();d=read(); 31 for(int i=1;i<=k+3;i++)g[i]=power(i,k); 32 for(int i=2;i<=k+3;i++)g[i]=(1ll*g[i-1]+g[i])%mod; 33 for(int i=2;i<=k+3;i++)g[i]=(1ll*g[i-1]+g[i])%mod; 34 for(int i=1;i<=k+3;i++) 35 { 36 w[i]=1; 37 for(int j=1;j<=k+3;j++) 38 if(i!=j)w[i]=(1ll*w[i]*(i-j)%mod+mod)%mod; 39 w[i]=power(w[i],mod-2); 40 } 41 for(int x=0;x<=k+4;x++) 42 { 43 f[x]=0;m=(a+1ll*x*d%mod)%mod; 44 for(int i=1;i<=k+3;i++) 45 { 46 sum=1ll*g[i]*w[i]%mod; 47 for(int j=1;j<=k+3;j++) 48 if(i!=j)sum=(1ll*sum*(m-j)%mod+mod)%mod; 49 f[x]=(1ll*f[x]+sum)%mod; 50 } 51 } 52 for(int i=1;i<=k+4;i++)f[i]=(1ll*f[i-1]+f[i])%mod; 53 for(int i=0;i<=k+4;i++) 54 { 55 w[i]=1; 56 for(int j=0;j<=k+4;j++) 57 if(i!=j)w[i]=(1ll*w[i]*(i-j)%mod+mod)%mod; 58 w[i]=power(w[i],mod-2); 59 } 60 ans=0; 61 for(int i=0;i<=k+4;i++) 62 { 63 sum=1ll*f[i]*w[i]%mod; 64 for(int j=0;j<=k+4;j++) 65 if(i!=j)sum=(1ll*sum*(n-j)%mod+mod)%mod; 66 ans=(1ll*ans+sum)%mod; 67 } 68 printf("%d ",ans); 69 } 70 int main() 71 { 72 T=read(); 73 while(T--)work(); 74 return 0; 75 }