FFT相關詳解：

http://blog.miskcoo.com/2015/04/polynomial-multiplication-and-fast-fourier-transform

https://blog.csdn.net/ggn_2015/article/details/68922404

練習題：

bzoj2179 FFT快速傅立葉:

模板題

#include<cstdio> #include<algorithm> #include<cstring> #include<cmath> #define inf 0x7f7f7f7f using namespace std; const int maxn=2e5; const double pi=acos(-1); int n,m,l; int sum[maxn*2+10],rev[maxn+10];int read() {char ch;int x=0,f=1;for (;ch<'0'||ch>'9';ch=getchar()) if (ch=='-') f=-1;for (;ch>='0'&&ch<='9';ch=getchar()) x=x*10+ch-'0';return x*f; }struct complex {double x,y;complex(){};complex(double _x,double _y){x=_x,y=_y;}complex operator +(const complex &a){return complex(x+a.x,y+a.y);}complex operator -(const complex &a){return complex(x-a.x,y-a.y);}complex operator *(const complex &a){return complex(x*a.x-y*a.y,x*a.y+y*a.x);} }A[maxn+10],B[maxn+10],Ans[maxn*2+10];int get() {int ch=getchar();for (;ch<'0'||ch>'9';ch=getchar());return ch-'0'; }void fft(complex *a,int len,int flag) {for (int i=0;i<len;i++) if (rev[i]<i) swap(a[i],a[rev[i]]);for (int i=2;i<=len;i<<=1){complex wn(cos(2*pi/i),sin(2*pi/i*flag));for (int j=0;j<len;j+=i){complex nw(1,0);for (int k=0;k<(i>>1);k++,nw=nw*wn){complex x=a[j+k],y=nw*a[j+k+(i>>1)];a[j+k]=x+y,a[j+k+(i>>1)]=x-y;}}}if (flag==-1) for (int i=0;i<len;i++) a[i].x/=len; }int main() {n=read();for (int i=n-1;~i;i--) A[i].x=get();for (int i=n-1;~i;i--) B[i].x=get();for (n*=2,m=1;m<=n;m<<=1)l++; for (int i=0;i<m;i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<(l-1));fft(A,m,1),fft(B,m,1);for (int i=0;i<m;i++) Ans[i]=A[i]*B[i];fft(Ans,m,-1);for (int i=0;i<m;i++) sum[i]=(int)(Ans[i].x+0.5);for (int i=0;i<m;i++){sum[i+1]+=sum[i]/10;sum[i]%=10;}while((!sum[n-1])&&n) n--;while(sum[n]) n++;for (int i=n-1;~i;i--) putchar(sum[i]+'0');printf("\n");return 0; }

bzoj2194 快速傅立葉之二:

模板題+1。計算\(C_k=\sum_{i=k}^{n-1}{a_i \times b_{i-k}}\),將\(b[]\)翻轉，則\(C_k=\sum_{i=k}^{n-1}{a_i \times b_{n+1+k-i}}\)，即\(C_k=\sum\limits_{i+j=n+1+k,0<=i,j<=2 \times n}{a_i \times b_j}\)。直接FFT即可。

#include<cstdio> #include<algorithm> #include<cstring> #include<cmath> #define inf 0x7f7f7f7f using namespace std; const int maxn=4e5; const double pi=acos(-1); int n,m,l; int rev[(maxn<<2)+5];int read() {char ch;int x=0,f=1;for (;ch<'0'||ch>'9';ch=getchar()) if (ch=='-') f=-1;for (;ch>='0'&&ch<='9';ch=getchar()) x=x*10+ch-'0';return x*f; }struct complex {double x,y;complex(){};complex(double _x,double _y){x=_x,y=_y;}complex operator +(const complex &a){return complex(x+a.x,y+a.y);};complex operator -(const complex &a){return complex(x-a.x,y-a.y);};complex operator *(const complex &a){return complex(x*a.x-y*a.y,x*a.y+y*a.x);}; }a[maxn+10],b[maxn+10],ans[(maxn<<2)+5];void fft(complex *a,int len,int flag) {for (int i=0;i<len;i++) if (rev[i]<i) swap(a[i],a[rev[i]]);for (int i=2;i<=len;i<<=1){complex wn(cos(2*pi/i),sin(2*pi/i*flag));for (int j=0;j<len;j+=i){complex nw(1,0);for (int k=0;k<(i>>1);k++,nw=nw*wn){complex x=a[j+k],y=nw*a[j+k+(i>>1)];a[j+k]=x+y,a[j+k+(i>>1)]=x-y;}}}if (flag==-1) for (int i=0;i<len;i++) a[i].x/=len; }int main() {n=read();for (int i=0;i<n;i++) a[i].x=read(),b[n+1-i].x=read();for (m=1;m<=n*2;m<<=1)l++;for (int i=0;i<m;i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<(l-1));fft(a,m,1),fft(b,m,1);for (int i=0;i<m;i++) ans[i]=a[i]*b[i];fft(ans,m,-1);for (int i=n+1;i<=n*2;i++) printf("%d\n",(int)(ans[i].x+0.5));return 0; }

bzoj4827 [Hnoi2017]禮物:

假設移動a位，整體增加c，則答案為\(\sum_{i=1}^{n}{(x_{((i+a-1) \% n+1)}-y_{i}+c)^{2}}\)。把此式展開\(\sum _{i=1} ^{n} {x_{( (i+a-1) \% n+1 )} ^{2} +y_{i} ^{2} +c ^{2} +2 \times c \times x_{( (i+a-1) \% n+1 )} -2 \times y_{i} \times c -2 \times x_{( (i+a-1) \% n+1)} \times y_{i}}\)。

\(\sum _{i=1} ^{n} {x_{( (i+a-1) \% n+1 )} \times y_{i}}\)可以用FFT求，之后枚舉a和c就行了。

#include<cstdio> #include<algorithm> #include<cstring> #include<cmath> #define ll long long #define inf 0x7f7f7f7f using namespace std; const int maxn=5e4; const double pi=acos(-1); int n,m,k,l,s,ans=1e9; int a[maxn+10],b[maxn+10],rev[maxn*4+10];int read() {char ch;int x=0,f=1;for (;ch<'0'||ch>'9';ch=getchar()) if (ch=='-') f=-1;for (;ch>='0'&&ch<='9';ch=getchar()) x=x*10+ch-'0';return x*f; }struct complex {double x,y;complex(){};complex (double _x,double _y){x=_x,y=_y;};complex operator +(const complex &a){return complex(x+a.x,y+a.y);}complex operator -(const complex &a){return complex(x-a.x,y-a.y);}complex operator *(const complex &a){return complex(x*a.x-y*a.y,x*a.y+y*a.x);} }x1[maxn*4+10],x2[maxn*4+10];void fft(complex *a,int len,int flag) {for (int i=0;i<len;i++) if (rev[i]<i) swap(a[i],a[rev[i]]);for (int i=2;i<=len;i<<=1){complex wn(cos(2*pi/i),sin(2*pi/i*flag));for (int j=0;j<len;j+=i){complex nw(1,0);for (int k=0;k<(i>>1);k++,nw=nw*wn){complex x=a[j+k],y=nw*a[j+k+(i>>1)];a[j+k]=x+y,a[j+k+(i>>1)]=x-y;}}}if (flag==-1) for (int i=0;i<len;i++) a[i].x/=len; }int main() {n=read(),k=read();for (int i=0;i<n;i++) a[i]=read(),x1[i].x=a[i];for (int i=0;i<n;i++) b[i]=read(),x2[n-i-1].x=b[i],s+=b[i];for (m=1;m<=n*2;m<<=1)l++;for (int i=0;i<m;i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<(l-1));fft(x1,m,1),fft(x2,m,1);for (int i=0;i<m;i++) x1[i]=x1[i]*x2[i];fft(x1,m,-1);for (int c=0;c<=k;c++){int sum=0;for (int i=0;i<n;i++) sum+=(a[i]+c)*(a[i]+c)+b[i]*b[i];ans=min(ans,sum-2*((int)(x1[n-1].x+0.5)+s*c));for (int i=1;i<n;i++){int tmp=(int)(x1[n+i-1].x+0.5)+(int)(x1[i-1].x+0.5)+s*c;ans=min(ans,sum-2*tmp);}}printf("%d\n",ans);return 0; }

bzoj4555 [Tjoi2016&Heoi2016]求和

求\(f_n=\sum_{i=0}^{n}\sum_{j=0}^{i}S_{i,j} \times 2^j \times j!\),對答案\(\%98244353\)。其中\(S_{i,j}\)為第二類斯特林數，意義為將n個不同的元素拆分成m個集合的方案數,所以當m>n時\(S_{n,m}\)為0。而\(S_{n,m} \times m!\)就是將n個不同的元素拆分成m個不同的集合的方案數。用容斥可得到\(S_{n,m} \times m!=\sum_{i=0}^{m}{(-1)^i \times \binom{m}{i} \times (m-i)^n}\)。

好，那我們開始推式子。
\(f_n=\sum_{i=0}^{n}\sum_{j=0}^{i}S_{i,j} \times 2^j \times j!\)
\(f_n=\sum_{i=0}^{n}\sum_{j=0}^{n}{2^j \times\sum_{k=0}^{j}}(-1)^k \times \binom{j}{k} \times (j-k)^i\)
\(f_n=\sum_{j=0}^{n}2^j \times \sum_{k=0}^{j}(-1)^k \times \frac{j!}{k!(j-k)!} \times \sum_{i=0}^{n} (j-k)^i?\)
\(f_n=\sum_{j=0}^{n} 2^j \times j! \times \sum_{k=0}^{j} \frac{(-1)^k}{k!} \times \frac{\sum_{i=0}^{n}(j-k)^i}{(j-k)!}\)

后面的\(\sum_{k=0}^{j} \frac{(-1)^k}{k!} \times \frac{\sum_{i=0}^{n}(j-k)^i}{(j-k)!}\)就是一個卷積的形式，因為要膜，所以直接上NTT即可。

#include<cstdio> #include<algorithm> #include<cstring> #include<cmath> #define ll long long #define inf 0x7f7f7f7f using namespace std; const int mod=998244353; const int maxn=1e5+5; int ans,n,m,l; int rev[maxn<<2],inv[maxn<<2],x[maxn<<2],y[maxn<<2];int read() {char ch;int x=0,f=1;for (;ch<'0'||ch>'9';ch=getchar()) if (ch=='-') f=-1;for (;ch>='0'&&ch<='9';ch=getchar()) x=x*10+ch-'0';return x*f; } int power(int a,int k) {int sum=1;for (;k;k>>=1,a=1ll*a*a%mod)if (k&1)sum=1ll*sum*a%mod;return sum; }void ntt(int *a,int len,int flag) {for (int i=0;i<len;i++) if (rev[i]<i) swap(a[rev[i]],a[i]);for (int i=2;i<=len;i<<=1){int gn=power(3,((mod-1)+flag*(mod-1)/i)%(mod-1));for (int j=0;j<len;j+=i){int ng=1;for (int k=0;k<(i>>1);k++,ng=1ll*ng*gn%mod){int x=a[j+k],y=1ll*ng*a[j+k+(i>>1)]%mod;a[j+k]=(x+y)%mod,a[j+k+(i>>1)]=1ll*((x-y)%mod+mod)%mod;}}}if (flag==-1){int invlen=power(len,mod-2);for (int i=0;i<len;i++) x[i]=1ll*x[i]*invlen%mod;} }int main() {n=read();for (m=1;m<n*2;m<<=1)l++;for (int i=0;i<m;i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<(l-1));inv[0]=1;for (int i=1;i<=n;i++) inv[i]=1ll*inv[i-1]*i%mod;inv[n]=power(inv[n],mod-2);for (int i=n-1;i;i--) inv[i]=1ll*inv[i+1]*(i+1)%mod;for (int i=0;i<=n;i++) x[i]=(inv[i]*(i&1?-1:1)+mod)%mod;y[0]=1,y[1]=n+1;for (int i=2;i<=n;i++) y[i]=1ll*(power(i,n+1)-1)*power(i-1,mod-2)%mod*inv[i]%mod;ntt(x,m,1),ntt(y,m,1);for (int i=0;i<m;i++) x[i]=1ll*x[i]*y[i]%mod;ntt(x,m,-1);int two_power=1,fact=1;ans=x[0];for (int i=1;i<=n;i++){two_power=1ll*two_power*2%mod,fact=1ll*fact*i%mod;ans=(ans+1ll*two_power*fact%mod*x[i]%mod)%mod;}printf("%d\n",ans);return 0; }

轉載于:https://www.cnblogs.com/Alseo_Roplyer/p/9449755.html

總結

以上是生活随笔為你收集整理的日常学习——FFT的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。