快速傅里葉變換(FFT)相關內容匯總

　　快速傅里葉變換，是求兩個多項式卷積的算法，其時間復雜度為$O(n\log n)$，優於普通卷積求法，且根據有關證明，快速傅里葉變換是基於變換求卷積的理論最快算法。

　　關於FFT的介紹，最詳細易懂的是《算法導論》上的內容。

　　其大致介紹與代碼在這里：http://www.cnblogs.com/rvalue/p/7351400.html.

 

1.FFT&NTT模板

#include<cmath>
#include<cstdio>
#include<algorithm>
#define rep(i,l,r) for (int i=l; i<=r; i++)
using namespace std;

const int N=300010;
const double pi=acos(-1.);
int n,m,L,x,rev[N];

struct C{
    double x,y;
    C(double _x=0,double _y=0):x(_x),y(_y){}
    C operator +(C &b){ return C(x+b.x,y+b.y); }
    C operator -(C &b){ return C(x-b.x,y-b.y); }
    C operator *(C &b){ return C(x*b.x-y*b.y,x*b.y+b.x*y); }
}a[N],b[N];

void DFT(C a[],int n,bool f){
    for (int i=0; i<n; i++) if (i<rev[i]) swap(a[i],a[rev[i]]);
    for (int i=1; i<n; i<<=1){
        C wn=C(cos(pi/i),(f?1:-1)*sin(pi/i));
        for (int p=i<<1,j=0; j<n; j+=p){
            C w(1,0);
            for (int k=0; k<i; k++,w=w*wn){
                C x=a[j+k],y=w*a[i+j+k]; a[j+k]=x+y; a[i+j+k]=x-y;
            }
        }
    }
    if (!f) for (int i=0; i<n; i++) a[i].x/=n;
}

int main(){
    freopen("Dft.in","r",stdin);
    freopen("Dft.out","w",stdout);
    scanf("%d%d",&n,&m);
    rep(i,0,n) scanf("%d",&x),a[i]=x;
    rep(i,0,m) scanf("%d",&x),b[i]=x;
    rep(i,0,max(n,m)) a[i]=C(a[i].x+b[i].x,a[i].x-b[i].x);
    m=n+m; for (n=1; n<=m; n<<=1) L++;
    rep(i,0,n-1) rev[i]=(rev[i>>1]>>1)|((i&1)<<(L-1));
    DFT(a,n,1); rep(i,0,n-1) a[i]=a[i]*a[i]; DFT(a,n,0);
    rep(i,0,m) printf("%d ",(int)(a[i].x/4+0.5));
    return 0;
}

void NTT(int a[],int n,bool f){
    for (int i=0; i<n; i++) if (i<rev[i]) swap(a[i],a[rev[i]]);
    for (int i=1; i<n; i<<=1){
        int wn=ksm(3,f ? (mod-1)/(i<<1) : (mod-1)-(mod-1)/(i<<1));
        for (int p=i<<1,j=0; j<n; j+=p){
            int w=1;
            for (int k=0; k<i; k++,w=1ll*w*wn%mod){
                int x=a[j+k],y=1ll*w*a[i+j+k]%mod;
                a[j+k]=(x+y)%mod; a[i+j+k]=(x-y+mod)%mod;
            }
        }
    }
    if (f) return;
    int inv=ksm(n,mod-2);
    for (int i=0; i<n; i++) a[i]=1ll*a[i]*inv%mod;
}

#include<cmath>
#include<cstdio>
#include<algorithm>
#include<cstring>
#define mem(a) memset(a,0,sizeof(a))
#define rep(i,l,r) for (int i=(l); i<=(r); i++)
using namespace std;

const int N=500010,mod=998244353,inv2=(mod+1)/2;
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];

void Print(int a[],int n=::n){ for (int i=0; i<n; i++) printf("%d ",a[i]); puts(""); }

int ksm(int a,int b){
    int res=1;
    for (; b; a=1ll*a*a%mod,b>>=1)
        if (b & 1) res=1ll*res*a%mod;
    return res;
}

void NTT(int a[],int n,bool f){
    for (int i=0; i<n; i++) if (i<rev[i]) swap(a[i],a[rev[i]]);
    for (int i=1; i<n; i<<=1){
        int wn=ksm(3,f ? (mod-1)/(i<<1) : (mod-1)-(mod-1)/(i<<1));
        for (int p=i<<1,j=0; j<n; j+=p){
            int w=1;
            for (int k=0; k<i; k++,w=1ll*w*wn%mod){
                int x=a[j+k],y=1ll*w*a[i+j+k]%mod;
                a[j+k]=(x+y)%mod; a[i+j+k]=(x-y+mod)%mod;
            }
        }
    }
    if (f) return;
    int inv=ksm(n,mod-2);
    for (int i=0; i<n; i++) a[i]=1ll*a[i]*inv%mod;
}

void mul(int a[],int b[],int l){
    int n=1,L=0;
    for (; n<(l<<1); n<<=1) L++;
    for (int i=0; i<n; i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<(L-1));
    NTT(a,n,1); NTT(b,n,1);
    for (int i=0; i<n; i++) a[i]=1ll*a[i]*b[i]%mod;
    NTT(a,n,0); NTT(b,n,0);
}

void Inv(int a[],int b[],int l){
    if (l==1){ b[0]=ksm(a[0],mod-2); return; }
    Inv(a,b,l>>1); int n=1,L=0;
    for (; n<(l<<1); n<<=1) L++;
    for (int i=0; i<n; i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<(L-1));
    for (int i=0; i<l; i++) A[i]=a[i];
    NTT(A,n,1); NTT(b,n,1);
    for (int i=0; i<n; i++) b[i]=1ll*b[i]*(2-1ll*A[i]*b[i]%mod+mod)%mod;
    NTT(b,n,0);
    for (int i=l; i<n; i++) b[i]=0;
    for (int i=0; i<n; i++) A[i]=0;
}

void Sqrt(int a[],int b[],int l){
    if (l==1){ b[0]=sqrt(a[0]); return; }
    Sqrt(a,b,l>>1); Inv(b,B,l); int n=1,L=0;
    for (; n<(l<<1); n<<=1) L++;
    for (int i=0; i<n; i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<(L-1));
    for (int i=0; i<l; i++) C[i]=a[i];
    NTT(C,n,1); NTT(B,n,1); NTT(b,n,1);
    for (int i=0; i<n; i++) b[i]=1ll*inv2*(b[i]+1ll*C[i]*B[i]%mod)%mod;
    NTT(b,n,0);
    for (int i=l; i<n; i++) b[i]=0;
    for (int i=0; i<n; i++) C[i]=B[i]=0;
}

void Deri(int a[],int b[],int l){
    for (int i=1; i<l; i++) b[i-1]=1ll*i*a[i]%mod;
}

void Inte(int a[],int b[],int l){
    for (int i=1; i<l; i++) b[i]=1ll*a[i-1]*inv[i]%mod; b[0]=0;
}

void Ln(int a[],int b[],int l){
    Deri(a,D,l); Inv(a,E,l); mul(D,E,l); Inte(D,b,l);
    for (int i=0; i<(l<<1); i++) D[i]=E[i]=0;
}

void Exp(int a[],int b[],int l){
    if (l==1){ b[0]=1; return; }
    Exp(a,b,l>>1); Ln(b,F,l); int n=1,L=0;
    for (; n<(l<<1); n<<=1) L++;
    for (int i=0; i<n; i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<(L-1));
    for (int i=0; i<l; i++) F[i]=(-F[i]+a[i]+mod)%mod; F[0]=(F[0]+1)%mod;
    NTT(F,n,1); NTT(b,n,1);
    for (int i=0; i<n; i++) b[i]=1ll*b[i]*F[i]%mod;
    NTT(b,n,0);
    for (int i=l; i<n; i++) b[i]=0;
    for (int i=0; i<n; i++) F[i]=0;
}

void Ksm(int a[],int b[],int k,int l){
    Ln(a,G,l);
    for (int i=0; i<l; i++) G[i]=1ll*k*G[i]%mod;
    Exp(G,b,l);
}

int main(){
    freopen("polynomial.in","r",stdin);
    freopen("polynomial.out","w",stdout);
    scanf("%d%d",&n,&k);
    for (int i=0; i<n; i++) scanf("%d",&X[i]);
    int l=1; for (; l<=n; l<<=1); inv[0]=inv[1]=1;
    rep(i,2,l) inv[i]=1ll*(mod-mod/i)*inv[mod%i]%mod;
    Sqrt(X,Y,l); mem(X); //Print(Y);
    Inv(Y,X,l); mem(Y); //Print(X);
    Inte(X,Y,l); mem(X); //Print(Y);
    Exp(Y,X,l); mem(Y); //Print(X);
    Inv(X,Y,l); Y[0]=(Y[0]+1)%mod; mem(X); //Print(Y);
    Ln(Y,X,l); X[0]=(X[0]+1)%mod; mem(Y); //Print(X);
    Ksm(X,Y,k,l); mem(X); //Print(Y);
    Deri(Y,X,n); mem(Y); Print(X);
    return 0;
}

 

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))$

void Inv(int a[],int b[],int l){
    if (l==1){ b[0]=ksm(a[0],mod-2); return; }
    Inv(a,b,l>>1); int n=1,L=0;
    for (; n<(l<<1); n<<=1) L++;
    for (int i=0; i<n; i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<(L-1));
    for (int i=0; i<l; i++) A[i]=a[i];
    NTT(A,n,1); NTT(b,n,1);
    for (int i=0; i<n; i++) b[i]=1ll*b[i]*(2-1ll*A[i]*b[i]%mod+mod)%mod;
    NTT(b,n,0);
    for (int i=l; i<n; i++) b[i]=0;
    for (int i=0; i<n; i++) A[i]=0;
}

 

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)$

void Sqrt(int a[],int b[],int l){
    if (l==1){ b[0]=sqrt(a[0]); return; }
    Sqrt(a,b,l>>1); Inv(b,B,l); int n=1,L=0;
    for (; n<(l<<1); n<<=1) L++;
    for (int i=0; i<n; i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<(L-1));
    for (int i=0; i<l; i++) C[i]=a[i];
    NTT(C,n,1); NTT(B,n,1); NTT(b,n,1);
    for (int i=0; i<n; i++) b[i]=1ll*inv2*(b[i]+1ll*C[i]*B[i]%mod)%mod;
    NTT(b,n,0);
    for (int i=l; i<n; i++) b[i]=0;
    for (int i=0; i<n; i++) C[i]=B[i]=0;
}

 

5.多項式求導與積分

這個很簡單，直接按定義來即可，復雜度線性。

void Deri(int a[],int b[],int l){
    for (int i=1; i<l; i++) b[i-1]=1ll*i*a[i]%mod;
}

void Inte(int a[],int b[],int l){
    for (int i=1; i<l; i++) b[i]=1ll*a[i-1]*inv[i]%mod; b[0]=0;
}

 

6.多項式求ln

這個不需要牛頓迭代，因為設$\ln(A(x))=B(x)$，求導得$\int\frac{A'(x)}{A(x)}=B(x)$，於是做一遍多項式求導、一遍求逆、一遍乘法和一遍積分即可。要求常數項為1。

void Ln(int a[],int b[],int l){
    Deri(a,D,l); Inv(a,E,l); mul(D,E,l); Inte(D,b,l);
    for (int i=0; i<(l<<1); i++) D[i]=E[i]=0;
}

 

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。

void Exp(int a[],int b[],int l){
    if (l==1){ b[0]=1; return; }
    Exp(a,b,l>>1); Ln(b,F,l); int n=1,L=0;
    for (; n<(l<<1); n<<=1) L++;
    for (int i=0; i<n; i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<(L-1));
    for (int i=0; i<l; i++) F[i]=(-F[i]+a[i]+mod)%mod; F[0]=(F[0]+1)%mod;
    NTT(F,n,1); NTT(b,n,1);
    for (int i=0; i<n; i++) b[i]=1ll*b[i]*F[i]%mod;
    NTT(b,n,0);
    for (int i=l; i<n; i++) b[i]=0;
    for (int i=0; i<n; i++) F[i]=0;
}

 

8.多項式快速冪

$A^k(x)=exp(k*\ln(A(x)))$

不管怎么套都是一個log，但顯然常數巨大。

void Ksm(int a[],int b[],int k,int l){
    Ln(a,G,l);
    for (int i=0; i<l; i++) G[i]=1ll*k*G[i]%mod;
    Exp(G,b,l);
}

 

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直接給出也是可以直接多項式求逆的）

void CDQ(int l,int r){
    if (l==r) return;
    int mid=(l+r)>>1,lim=r-l+1,n=1,L=0;
    CDQ(l,mid);
    while (n<lim) n<<=1,L++;
    for (int i=0; i<n; i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<(L-1));
    for (int i=0; i<n; i++) a[i]=b[i]=0;
    rep(i,l,mid) a[i-l]=f[i];
    rep(i,0,r-l) b[i]=g[i];
    NTT(a,n,1); NTT(b,n,1);
    for (int i=0; i<n; i++) a[i]=1ll*a[i]*b[i]%mod;
    NTT(a,n,0);
    rep(i,mid+1,r) f[i]=(f[i]+a[i-l])%mod;
    CDQ(mid+1,r);
}

 

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。

void calc(int a1[],int a2[],int l){
    int n=1,L=0;
    for (; n<l; n<<=1,L++);
    for (int i=0; i<n; i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<(L-1));
    for (int i=0; i<n; i++) A[i]=a1[i]>>15,B[i]=a1[i]&0x7fff;
    for (int i=0; i<n; i++) C[i]=a2[i]>>15,D[i]=a2[i]&0x7fff;
    FFT(A,n,1); FFT(B,n,1); FFT(C,n,1); FFT(D,n,1);
    for (int i=0; i<n; i++) s1[i]=A[i]*C[i],s2[i]=A[i]*D[i]+B[i]*C[i],s3[i]=B[i]*D[i];
    FFT(s1,n,-1); FFT(s2,n,-1); FFT(s3,n,-1);
    for (int i=0; i<n; i++)
        a1[i]=((((ll)round(s1[i].x)%mod)<<30)%mod+(((ll)round(s2[i].x)%mod)<<15)%mod+(ll)round(s3[i].x)%mod)%mod;
}

 

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