快速傅里叶变换(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