快速傅里叶变换(FFT)相关内容汇总
（原稿：https://paste.ubuntu.com/p/yJNsn3xPt8/）
快速傅里葉變換，是求兩個多項式卷積的算法，其時間復雜度為$O(n\log n)$，優于普通卷積求法，且根據有關證明，快速傅里葉變換是基于變換求卷積的理論最快算法。
關于FFT的介紹，最詳細易懂的是《算法導論》上的內容。
其大致介紹與代碼在這里：http://www.cnblogs.com/rvalue/p/7351400.html.
?
1.FFT&NTT模板
1 #include<cmath> 2 #include<cstdio> 3 #include<algorithm> 4 #define rep(i,l,r) for (int i=l; i<=r; i++) 5 using namespace std; 6 7 const int N=300010; 8 const double pi=acos(-1.); 9 int n,m,L,x,rev[N]; 10 11 struct C{ 12 double x,y; 13 C(double _x=0,double _y=0):x(_x),y(_y){} 14 C operator +(C &b){ return C(x+b.x,y+b.y); } 15 C operator -(C &b){ return C(x-b.x,y-b.y); } 16 C operator *(C &b){ return C(x*b.x-y*b.y,x*b.y+b.x*y); } 17 }a[N],b[N]; 18 19 void DFT(C a[],int n,bool f){ 20 for (int i=0; i<n; i++) if (i<rev[i]) swap(a[i],a[rev[i]]); 21 for (int i=1; i<n; i<<=1){ 22 C wn=C(cos(pi/i),(f?1:-1)*sin(pi/i)); 23 for (int p=i<<1,j=0; j<n; j+=p){ 24 C w(1,0); 25 for (int k=0; k<i; k++,w=w*wn){ 26 C x=a[j+k],y=w*a[i+j+k]; a[j+k]=x+y; a[i+j+k]=x-y; 27 } 28 } 29 } 30 if (!f) for (int i=0; i<n; i++) a[i].x/=n; 31 } 32 33 int main(){ 34 freopen("Dft.in","r",stdin); 35 freopen("Dft.out","w",stdout); 36 scanf("%d%d",&n,&m); 37 rep(i,0,n) scanf("%d",&x),a[i]=x; 38 rep(i,0,m) scanf("%d",&x),b[i]=x; 39 rep(i,0,max(n,m)) a[i]=C(a[i].x+b[i].x,a[i].x-b[i].x); 40 m=n+m; for (n=1; n<=m; n<<=1) L++; 41 rep(i,0,n-1) rev[i]=(rev[i>>1]>>1)|((i&1)<<(L-1)); 42 DFT(a,n,1); rep(i,0,n-1) a[i]=a[i]*a[i]; DFT(a,n,0); 43 rep(i,0,m) printf("%d ",(int)(a[i].x/4+0.5)); 44 return 0; 45 } FFT 1 void NTT(int a[],int n,bool f){ 2 for (int i=0; i<n; i++) if (i<rev[i]) swap(a[i],a[rev[i]]); 3 for (int i=1; i<n; i<<=1){ 4 int wn=ksm(3,f ? (mod-1)/(i<<1) : (mod-1)-(mod-1)/(i<<1)); 5 for (int p=i<<1,j=0; j<n; j+=p){ 6 int w=1; 7 for (int k=0; k<i; k++,w=1ll*w*wn%mod){ 8 int x=a[j+k],y=1ll*w*a[i+j+k]%mod; 9 a[j+k]=(x+y)%mod; a[i+j+k]=(x-y+mod)%mod; 10 } 11 } 12 } 13 if (f) return; 14 int inv=ksm(n,mod-2); 15 for (int i=0; i<n; i++) a[i]=1ll*a[i]*inv%mod; 16 } NTT 1 #include<cmath> 2 #include<cstdio> 3 #include<algorithm> 4 #include<cstring> 5 #define mem(a) memset(a,0,sizeof(a)) 6 #define rep(i,l,r) for (int i=(l); i<=(r); i++) 7 using namespace std; 8 9 const int N=500010,mod=998244353,inv2=(mod+1)/2; 10 int n,k,rev[N],inv[N],X[N],Y[N],A[N],B[N],C[N],D[N],E[N],F[N],G[N]; 11 12 void Print(int a[],int n=::n){ for (int i=0; i<n; i++) printf("%d ",a[i]); puts(""); } 13 14 int ksm(int a,int b){ 15 int res=1; 16 for (; b; a=1ll*a*a%mod,b>>=1) 17 if (b & 1) res=1ll*res*a%mod; 18 return res; 19 } 20 21 void NTT(int a[],int n,bool f){ 22 for (int i=0; i<n; i++) if (i<rev[i]) swap(a[i],a[rev[i]]); 23 for (int i=1; i<n; i<<=1){ 24 int wn=ksm(3,f ? (mod-1)/(i<<1) : (mod-1)-(mod-1)/(i<<1)); 25 for (int p=i<<1,j=0; j<n; j+=p){ 26 int w=1; 27 for (int k=0; k<i; k++,w=1ll*w*wn%mod){ 28 int x=a[j+k],y=1ll*w*a[i+j+k]%mod; 29 a[j+k]=(x+y)%mod; a[i+j+k]=(x-y+mod)%mod; 30 } 31 } 32 } 33 if (f) return; 34 int inv=ksm(n,mod-2); 35 for (int i=0; i<n; i++) a[i]=1ll*a[i]*inv%mod; 36 } 37 38 void mul(int a[],int b[],int l){ 39 int n=1,L=0; 40 for (; n<(l<<1); n<<=1) L++; 41 for (int i=0; i<n; i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<(L-1)); 42 NTT(a,n,1); NTT(b,n,1); 43 for (int i=0; i<n; i++) a[i]=1ll*a[i]*b[i]%mod; 44 NTT(a,n,0); NTT(b,n,0); 45 } 46 47 void Inv(int a[],int b[],int l){ 48 if (l==1){ b[0]=ksm(a[0],mod-2); return; } 49 Inv(a,b,l>>1); int n=1,L=0; 50 for (; n<(l<<1); n<<=1) L++; 51 for (int i=0; i<n; i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<(L-1)); 52 for (int i=0; i<l; i++) A[i]=a[i]; 53 NTT(A,n,1); NTT(b,n,1); 54 for (int i=0; i<n; i++) b[i]=1ll*b[i]*(2-1ll*A[i]*b[i]%mod+mod)%mod; 55 NTT(b,n,0); 56 for (int i=l; i<n; i++) b[i]=0; 57 for (int i=0; i<n; i++) A[i]=0; 58 } 59 60 void Sqrt(int a[],int b[],int l){ 61 if (l==1){ b[0]=sqrt(a[0]); return; } 62 Sqrt(a,b,l>>1); Inv(b,B,l); int n=1,L=0; 63 for (; n<(l<<1); n<<=1) L++; 64 for (int i=0; i<n; i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<(L-1)); 65 for (int i=0; i<l; i++) C[i]=a[i]; 66 NTT(C,n,1); NTT(B,n,1); NTT(b,n,1); 67 for (int i=0; i<n; i++) b[i]=1ll*inv2*(b[i]+1ll*C[i]*B[i]%mod)%mod; 68 NTT(b,n,0); 69 for (int i=l; i<n; i++) b[i]=0; 70 for (int i=0; i<n; i++) C[i]=B[i]=0; 71 } 72 73 void Deri(int a[],int b[],int l){ 74 for (int i=1; i<l; i++) b[i-1]=1ll*i*a[i]%mod; 75 } 76 77 void Inte(int a[],int b[],int l){ 78 for (int i=1; i<l; i++) b[i]=1ll*a[i-1]*inv[i]%mod; b[0]=0; 79 } 80 81 void Ln(int a[],int b[],int l){ 82 Deri(a,D,l); Inv(a,E,l); mul(D,E,l); Inte(D,b,l); 83 for (int i=0; i<(l<<1); i++) D[i]=E[i]=0; 84 } 85 86 void Exp(int a[],int b[],int l){ 87 if (l==1){ b[0]=1; return; } 88 Exp(a,b,l>>1); Ln(b,F,l); int n=1,L=0; 89 for (; n<(l<<1); n<<=1) L++; 90 for (int i=0; i<n; i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<(L-1)); 91 for (int i=0; i<l; i++) F[i]=(-F[i]+a[i]+mod)%mod; F[0]=(F[0]+1)%mod; 92 NTT(F,n,1); NTT(b,n,1); 93 for (int i=0; i<n; i++) b[i]=1ll*b[i]*F[i]%mod; 94 NTT(b,n,0); 95 for (int i=l; i<n; i++) b[i]=0; 96 for (int i=0; i<n; i++) F[i]=0; 97 } 98 99 void Ksm(int a[],int b[],int k,int l){ 100 Ln(a,G,l); 101 for (int i=0; i<l; i++) G[i]=1ll*k*G[i]%mod; 102 Exp(G,b,l); 103 } 104 105 int main(){ 106 freopen("polynomial.in","r",stdin); 107 freopen("polynomial.out","w",stdout); 108 scanf("%d%d",&n,&k); 109 for (int i=0; i<n; i++) scanf("%d",&X[i]); 110 int l=1; for (; l<=n; l<<=1); inv[0]=inv[1]=1; 111 rep(i,2,l) inv[i]=1ll*(mod-mod/i)*inv[mod%i]%mod; 112 Sqrt(X,Y,l); mem(X); //Print(Y); 113 Inv(Y,X,l); mem(Y); //Print(X); 114 Inte(X,Y,l); mem(X); //Print(Y); 115 Exp(Y,X,l); mem(Y); //Print(X); 116 Inv(X,Y,l); Y[0]=(Y[0]+1)%mod; mem(X); //Print(Y); 117 Ln(Y,X,l); X[0]=(X[0]+1)%mod; mem(Y); //Print(X); 118 Ksm(X,Y,k,l); mem(X); //Print(Y); 119 Deri(Y,X,n); mem(Y); Print(X); 120 return 0; 121 } 多項式全家桶(COGS2189)?
2.牛頓迭代法
現已知$G(F(x))\equiv 0\ (\text{mod}\ x^n)$，考慮如何倍增求出$F(x)$。
我們先將x擴展到2的次冪項（多余位補0），方便后面的倍增。
當n=1時，直接算出常數項即可。
否則，先遞歸下去求出$F_0(x)$使得$G(F_0(x))\equiv 0\ (\text{mod}\ x^\frac{n}{2})$。
接著將$G(F(x))$在$F_0(x)$處泰勒展開，有：$$G(F(x))=G(F_0(x))+G'(F_0(x))(F(x)-F_0(x))+\frac{G''(F_0(x))}{2}(F(x)-F_0(x))^2+...$$由于$F(x)$與$F_0(x)$前$\frac{n}{2}$項都相同，所以上式從第三項開始也都是0。
于是有：$G(F(x))\equiv G(F_0(x))+G'(F_0(x))(F(x)-F_0(x))\equiv 0(\text{mod}\ x^n)$
移項得：$F(x)=F_0(x)-\frac{G(F_0(x))}{G'(F_0(x))}\ (\text{mod}\ x^n)$
（順便一提，這里的$F_0(x)$只是一個多項式，可以當作普通泰勒展開中的$x$看待，不要混淆概念將$G'(F_0(x))$鏈式法則了）
這樣，我們遞歸下去求出$F_0(x)$，再多項式求逆與乘法，就得到了$F(x)$。
時間復雜度：$T(n)=T(n/2)+O(n\log n)=O(n\log n)$
?
3.多項式求逆
利用上面牛頓迭代法，考慮如何求出一個多項式在模意義下的逆元，即對于多項式$A(x)$，求出$B(x)$使得$A(x)B(x)\equiv\ 1(\text{mod}\ x^n)$（當然這里不用牛頓迭代法也能得出相同的結論）。
構造函數$G(B(x))=A(x)B(x)-1$，則有$G(B(x))\equiv 0\ (\text{mod}\ x^n)$
同樣將次數界擴展到2的次冪，然后先遞歸下去求出$B_0(x)$使得$A(x)B_0(x)\equiv 1\ (\text{mod}\ x^\frac{n}{2})$。
接下來暴力帶入牛頓迭代的式子：$B(x)=B_0(x)-\frac{G(B_0(x))}{G'(B_0(x))}=B_0(x)-\frac{A(x)B_0(x)-1}{A(x)}=B_0(x)-B_0(x)(A(x)B_0(x)-1)=B_0(x)(2-A(x)B_0(x))$
?
4.多項式開根
和求逆一樣套牛頓迭代的式子：設$G(B(x))=B^2(x)-A(x)$，則遞歸下去求出$B_0(x)$，再$B(x)\equiv B_0(x)-\frac{G(B_0(x))}{G'(B_0(x))}=B_0(x)-\frac{B_0^2(x)-A(x)}{2B_0(x)}=\frac{B_0^2(x)+A(x)}{2B_0(x)}\ (\text{mod}\ x^n)$
1 void Sqrt(int a[],int b[],int l){ 2 if (l==1){ b[0]=sqrt(a[0]); return; } 3 Sqrt(a,b,l>>1); Inv(b,B,l); int n=1,L=0; 4 for (; n<(l<<1); n<<=1) L++; 5 for (int i=0; i<n; i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<(L-1)); 6 for (int i=0; i<l; i++) C[i]=a[i]; 7 NTT(C,n,1); NTT(B,n,1); NTT(b,n,1); 8 for (int i=0; i<n; i++) b[i]=1ll*inv2*(b[i]+1ll*C[i]*B[i]%mod)%mod; 9 NTT(b,n,0); 10 for (int i=l; i<n; i++) b[i]=0; 11 for (int i=0; i<n; i++) C[i]=B[i]=0; 12 } Sqrt?
5.多項式求導與積分
這個很簡單，直接按定義來即可，復雜度線性。
1 void Deri(int a[],int b[],int l){ 2 for (int i=1; i<l; i++) b[i-1]=1ll*i*a[i]%mod; 3 } 4 5 void Inte(int a[],int b[],int l){ 6 for (int i=1; i<l; i++) b[i]=1ll*a[i-1]*inv[i]%mod; b[0]=0; 7 } Derivation&Integration?
6.多項式求ln
這個不需要牛頓迭代，因為設$\ln(A(x))=B(x)$，求導得$\int\frac{A'(x)}{A(x)}=B(x)$，于是做一遍多項式求導、一遍求逆、一遍乘法和一遍積分即可。要求常數項為1。
1 void Ln(int a[],int b[],int l){ 2 Deri(a,D,l); Inv(a,E,l); mul(D,E,l); Inte(D,b,l); 3 for (int i=0; i<(l<<1); i++) D[i]=E[i]=0; 4 } Ln?
7.多項式求exp
$exp(A(x))=B(x)$，兩邊求導得$A(x)=\ln(B(x))$，設$G(B(x))=\ln(B(x))-A(x)$，遞歸下去求出$B_0(x)$，再$B(x)=B_0(x)-\frac{G(B_0(x))}{G'(B_0(x))}=B_0(x)-\frac{\ln(B(x))-A(x)}{\frac{1}{B_0(x)}}$=$B_0(x)(1-\ln(B(x))+A(x))$。要求常數項為0。
1 void Exp(int a[],int b[],int l){ 2 if (l==1){ b[0]=1; return; } 3 Exp(a,b,l>>1); Ln(b,F,l); int n=1,L=0; 4 for (; n<(l<<1); n<<=1) L++; 5 for (int i=0; i<n; i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<(L-1)); 6 for (int i=0; i<l; i++) F[i]=(-F[i]+a[i]+mod)%mod; F[0]=(F[0]+1)%mod; 7 NTT(F,n,1); NTT(b,n,1); 8 for (int i=0; i<n; i++) b[i]=1ll*b[i]*F[i]%mod; 9 NTT(b,n,0); 10 for (int i=l; i<n; i++) b[i]=0; 11 for (int i=0; i<n; i++) F[i]=0; 12 } Exp?
8.多項式快速冪
$A^k(x)=exp(k*\ln(A(x)))$
不管怎么套都是一個log，但顯然常數巨大。
1 void Ksm(int a[],int b[],int k,int l){ 2 Ln(a,G,l); 3 for (int i=0; i<l; i++) G[i]=1ll*k*G[i]%mod; 4 Exp(G,b,l); 5 } Pow?
9.多項式除法
(引用這里)
要求的就是給定兩個多項式$A(x)$,$B(x)$，其項數為$n$,$m$
求解一個$n-m$項的多項式$C(x)$，以及一個小于$n-m$項的多項式$R(x)$。
滿足：$A(x)=B(x)*C(x)+R(x)$。
定義一個操作$R$,其中$R$就是$Reverse$，$A^R(x)=x^nA(\frac{1}{x})$。這個操作說白點就是$A(x)[x^i]$對應$A^R(x)[x^{n-i}]$，也就是把所有的系數給翻轉過來。
然后推式子：
$$\begin{aligned} A(x)&=B(x)*C(x)+R(x)\\ A(\frac{1}{x})&=B(\frac{1}{x})*C(\frac{1}{x})+R(\frac{1}{x})\\ x^nA(\frac{1}{x})&=(x^m*B(\frac{1}{x}))*(x^{n-m}*C(\frac{1}{x}))+x^nR(\frac{1}{x})\\ A^R(x)&=B^R(x)*C^R(x)+R^R(x)*x^{n-m+1} \end{aligned}$$
到了這里我們把等號換成同余，把整個式子在模$x^{n-m+1}$意義下進行。
$$\begin{aligned} A^R(x)&\equiv B^R(x)*C^R(x)+R^R(x)*x^{n-m+1}\\ &\equiv B^R(x)*C^R(x) \end{aligned}$$
所以在模意義下，我們可以利用多項式求逆求解$C^R(x)=\frac{A^R(x)}{B^R(x)}mod\ x^{n-m+1}$
那么就有$R(x)=A(x)-B(x)*C(x)$。
?
10.分治FFT
分治FFT是求這樣一類卷積：現已知$g[0],g[1],...,g[n-1]$和$f[0]$，$f[i]=\sum\limits_{j=1}^{i}f[i-j]g[j]$，求$f[1],...,f[n]$。
由于式子左右都有$f$，需要用類似$CDQ$分治的做法來實現。復雜度$O(n\log^2 n)$
假設我們現在需要求$f[l],...,f[r]$先遞歸到左邊處理出$f[l],...,f[mid]$，再計算左半邊對右半邊的影響。這里直接卷上$g$即可，主要要作適當移位。（當然如果g直接給出也是可以直接多項式求逆的）
?
11.MTT
任意模數FFT，在模數不為998244353，直接做FFT又可能炸精度的情況下使用。
有兩種寫法，一種是myy的三模數FFT，一種是下面的做7次FFT的做法。(引用這里)
求$F(x)$與$G(x)$在任意模數下的卷積。
為什么不能直接FFT乘然后再取模？因為直接乘結果會爆long long。
考慮拆系數。設一個常數M，把$F(x)$和$G(x)$拆成：$A(x)=\frac{F(x)}M$，$B(x)=F(x)\%M$，$C(x)=\frac{G(x)}M$，$D(x)=G(x)\%M$。
這樣答案就變成了：
$M^2A(x)C(x)+M(A(x)D(x)+B(x)C(x))+B(x)D(x)$
直接FFT就不會爆long long了。$M$取$\sqrt {mod}$（我取了32768）最優。
這種方法一共要做7次DFT。
?
12.應用
?1.加速卷積運算
卷積運算一般會用在容斥、DP等題目中，所有形如$c_{k}=\sum_{i=0}^{k}a_{i}b_{k-i}$的計算都可以用FFT加速。同時，形如$c_{k}=\sum_{i=k}^{n-1}a_{i}b_{i-k}$的問題也可以，注意到只需要將$b[]$翻轉一下做一次卷積，然后$a[n],...,a[2n-1]$的位置存的就是$c[0],...,c[n-1]$的答案了，同時分治FFT也能解決相同形式的問題。
另外，在快速求解第一、第二類斯特林數時也有應用，見這里。
例題：HDU4609,BZOJ5306,BZOJ4555
?
2.加速字符串匹配
比如要求對一個字符串的每一個位置，統計以它為軸的最長回文子序列長度。
對字符集里的每種字符做一遍，所有等于當前處理字符的位置賦為1，否則賦為0。以x軸回文意味著s[x-i]=s[x+i]，注意到這時s[x-i]*s[x+i]為1，且如果做卷積，這個值會被統計到一個固定位置x-i+x+i=2x處。所以做字符集大小次FFT即可。
同理，普通字符串匹配也可以類似地做，但卷積后每個位置記錄的就是固定的差的信息了。
關于帶通配符的字符串匹配問題也有FFT應用：BZOJ3160,BZOJ4503。
?
3.優化DP轉移
一個比較常見的套路是，需要求一個滿足題設條件的DP數組f[]，可以先放寬題設限制求一個數組g[]，然后講f和g都視為生成函數，根據二者的關系做FFT。這種方法往往出現于帶有組合數、容斥的DP中。例題：BZOJ3456,Luogu5162
?
4.分治+FFT
求解形如$\prod(x^{k_i}+a_i)$且$\sum k_i$在1e5左右范圍內時可以應用的一種做法。
分治，遞歸到兩邊分別求解，然后一次卷積合并。$O(n\log^2 n)$
第一類斯特林數的快速求解方法之一。注意維護好次數界，這是復雜度的保證。
應用：LOJ2541,CF960G
?
5.生成函數
很深的部分，一個比較簡單的應用是：在做一些生成函數模型較明顯的計數問題時，可以通過等比數列求和公式、泰勒展開式將無窮項式化為閉形式，再通過數學知識快速求出某次項系數。
比較基礎的應用：
OGF（普通型生成函數）的冪次表示生成序列，EGF（指數型生成函數）的冪次表示生成集合。
EGF就是在OGF的基礎上多了卷積是附帶的C(n,i)，所以可以做多個帶標號基本元素的有序排列。
如：輪換計數、樹的計數、環套樹的計數、無向聯通圖計數等。
應用：POJ3734
?
6.其它
我不會的黑科技們
(1)多項式多點求值與快速插值
(2)常系數齊次線性遞推
(3)拉格朗日反演
(4)樹、圖的EGF
?
轉載于:https://www.cnblogs.com/HocRiser/p/8207295.html
總結
以上是生活随笔為你收集整理的快速傅里叶变换(FFT)相关内容汇总的全部內容，希望文章能夠幫你解決所遇到的問題。
- 上一篇： 遍历map的方法
- 下一篇： javas的四种状态 无锁状态 偏向锁状