快速傅里叶变换,是求两个多项式卷积的算法,其时间复杂度为$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=;
const double pi=acos(-.);
int n,m,L,x,rev[N]; struct C{
double x,y;
C(double _x=,double _y=):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=; i<n; i++) if (i<rev[i]) swap(a[i],a[rev[i]]);
for (int i=; i<n; i<<=){
C wn=C(cos(pi/i),(f?:-)*sin(pi/i));
for (int p=i<<,j=; j<n; j+=p){
C w(,);
for (int k=; 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=; 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,,n) scanf("%d",&x),a[i]=x;
rep(i,,m) scanf("%d",&x),b[i]=x;
rep(i,,max(n,m)) a[i]=C(a[i].x+b[i].x,a[i].x-b[i].x);
m=n+m; for (n=; n<=m; n<<=) L++;
rep(i,,n-) rev[i]=(rev[i>>]>>)|((i&)<<(L-));
DFT(a,n,); rep(i,,n-) a[i]=a[i]*a[i]; DFT(a,n,);
rep(i,,m) printf("%d ",(int)(a[i].x/+0.5));
return ;
}
FFT
void NTT(int a[],int n,bool f){
for (int i=; i<n; i++) if (i<rev[i]) swap(a[i],a[rev[i]]);
for (int i=; i<n; i<<=){
int wn=ksm(,f ? (mod-)/(i<<) : (mod-)-(mod-)/(i<<));
for (int p=i<<,j=; j<n; j+=p){
int w=;
for (int k=; 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-);
for (int i=; i<n; i++) a[i]=1ll*a[i]*inv%mod;
}
NTT
#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=,mod=,inv2=(mod+)/;
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=; i<n; i++) printf("%d ",a[i]); puts(""); } int ksm(int a,int b){
int res=;
for (; b; a=1ll*a*a%mod,b>>=)
if (b & ) res=1ll*res*a%mod;
return res;
} void NTT(int a[],int n,bool f){
for (int i=; i<n; i++) if (i<rev[i]) swap(a[i],a[rev[i]]);
for (int i=; i<n; i<<=){
int wn=ksm(,f ? (mod-)/(i<<) : (mod-)-(mod-)/(i<<));
for (int p=i<<,j=; j<n; j+=p){
int w=;
for (int k=; 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-);
for (int i=; i<n; i++) a[i]=1ll*a[i]*inv%mod;
} void mul(int a[],int b[],int l){
int n=,L=;
for (; n<(l<<); n<<=) L++;
for (int i=; i<n; i++) rev[i]=(rev[i>>]>>)|((i&)<<(L-));
NTT(a,n,); NTT(b,n,);
for (int i=; i<n; i++) a[i]=1ll*a[i]*b[i]%mod;
NTT(a,n,); NTT(b,n,);
} void Inv(int a[],int b[],int l){
if (l==){ b[]=ksm(a[],mod-); return; }
Inv(a,b,l>>); int n=,L=;
for (; n<(l<<); n<<=) L++;
for (int i=; i<n; i++) rev[i]=(rev[i>>]>>)|((i&)<<(L-));
for (int i=; i<l; i++) A[i]=a[i];
NTT(A,n,); NTT(b,n,);
for (int i=; i<n; i++) b[i]=1ll*b[i]*(-1ll*A[i]*b[i]%mod+mod)%mod;
NTT(b,n,);
for (int i=l; i<n; i++) b[i]=;
for (int i=; i<n; i++) A[i]=;
} void Sqrt(int a[],int b[],int l){
if (l==){ b[]=sqrt(a[]); return; }
Sqrt(a,b,l>>); Inv(b,B,l); int n=,L=;
for (; n<(l<<); n<<=) L++;
for (int i=; i<n; i++) rev[i]=(rev[i>>]>>)|((i&)<<(L-));
for (int i=; i<l; i++) C[i]=a[i];
NTT(C,n,); NTT(B,n,); NTT(b,n,);
for (int i=; i<n; i++) b[i]=1ll*inv2*(b[i]+1ll*C[i]*B[i]%mod)%mod;
NTT(b,n,);
for (int i=l; i<n; i++) b[i]=;
for (int i=; i<n; i++) C[i]=B[i]=;
} void Deri(int a[],int b[],int l){
for (int i=; i<l; i++) b[i-]=1ll*i*a[i]%mod;
} void Inte(int a[],int b[],int l){
for (int i=; i<l; i++) b[i]=1ll*a[i-]*inv[i]%mod; b[]=;
} 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=; i<(l<<); i++) D[i]=E[i]=;
} void Exp(int a[],int b[],int l){
if (l==){ b[]=; return; }
Exp(a,b,l>>); Ln(b,F,l); int n=,L=;
for (; n<(l<<); n<<=) L++;
for (int i=; i<n; i++) rev[i]=(rev[i>>]>>)|((i&)<<(L-));
for (int i=; i<l; i++) F[i]=(-F[i]+a[i]+mod)%mod; F[]=(F[]+)%mod;
NTT(F,n,); NTT(b,n,);
for (int i=; i<n; i++) b[i]=1ll*b[i]*F[i]%mod;
NTT(b,n,);
for (int i=l; i<n; i++) b[i]=;
for (int i=; i<n; i++) F[i]=;
} void Ksm(int a[],int b[],int k,int l){
Ln(a,G,l);
for (int i=; 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=; i<n; i++) scanf("%d",&X[i]);
int l=; for (; l<=n; l<<=); inv[]=inv[]=;
rep(i,,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[]=(Y[]+)%mod; mem(X); //Print(Y);
Ln(Y,X,l); X[]=(X[]+)%mod; mem(Y); //Print(X);
Ksm(X,Y,k,l); mem(X); //Print(Y);
Deri(Y,X,n); mem(Y); Print(X);
return ;
}
多项式全家桶(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))$
void Inv(int a[],int b[],int l){
if (l==){ b[]=ksm(a[],mod-); return; }
Inv(a,b,l>>); int n=,L=;
for (; n<(l<<); n<<=) L++;
for (int i=; i<n; i++) rev[i]=(rev[i>>]>>)|((i&)<<(L-));
for (int i=; i<l; i++) A[i]=a[i];
NTT(A,n,); NTT(b,n,);
for (int i=; i<n; i++) b[i]=1ll*b[i]*(-1ll*A[i]*b[i]%mod+mod)%mod;
NTT(b,n,);
for (int i=l; i<n; i++) b[i]=;
for (int i=; i<n; i++) A[i]=;
}
Inv
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==){ b[]=sqrt(a[]); return; }
Sqrt(a,b,l>>); Inv(b,B,l); int n=,L=;
for (; n<(l<<); n<<=) L++;
for (int i=; i<n; i++) rev[i]=(rev[i>>]>>)|((i&)<<(L-));
for (int i=; i<l; i++) C[i]=a[i];
NTT(C,n,); NTT(B,n,); NTT(b,n,);
for (int i=; i<n; i++) b[i]=1ll*inv2*(b[i]+1ll*C[i]*B[i]%mod)%mod;
NTT(b,n,);
for (int i=l; i<n; i++) b[i]=;
for (int i=; i<n; i++) C[i]=B[i]=;
}
Sqrt
5.多项式求导与积分
这个很简单,直接按定义来即可,复杂度线性。
void Deri(int a[],int b[],int l){
for (int i=; i<l; i++) b[i-]=1ll*i*a[i]%mod;
} void Inte(int a[],int b[],int l){
for (int i=; i<l; i++) b[i]=1ll*a[i-]*inv[i]%mod; b[]=;
}
Derivation&Integration
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=; i<(l<<); i++) D[i]=E[i]=;
}
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。
void Exp(int a[],int b[],int l){
if (l==){ b[]=; return; }
Exp(a,b,l>>); Ln(b,F,l); int n=,L=;
for (; n<(l<<); n<<=) L++;
for (int i=; i<n; i++) rev[i]=(rev[i>>]>>)|((i&)<<(L-));
for (int i=; i<l; i++) F[i]=(-F[i]+a[i]+mod)%mod; F[]=(F[]+)%mod;
NTT(F,n,); NTT(b,n,);
for (int i=; i<n; i++) b[i]=1ll*b[i]*F[i]%mod;
NTT(b,n,);
for (int i=l; i<n; i++) b[i]=;
for (int i=; i<n; i++) F[i]=;
}
Exp
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=; i<l; i++) G[i]=1ll*k*G[i]%mod;
Exp(G,b,l);
}
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直接给出也是可以直接多项式求逆的)
void CDQ(int l,int r){
if (l==r) return;
int mid=(l+r)>>,lim=r-l+,n=,L=;
CDQ(l,mid);
while (n<lim) n<<=,L++;
for (int i=; i<n; i++) rev[i]=(rev[i>>]>>)|((i&)<<(L-));
for (int i=; i<n; i++) a[i]=b[i]=;
rep(i,l,mid) a[i-l]=f[i];
rep(i,,r-l) b[i]=g[i];
NTT(a,n,); NTT(b,n,);
for (int i=; i<n; i++) a[i]=1ll*a[i]*b[i]%mod;
NTT(a,n,);
rep(i,mid+,r) f[i]=(f[i]+a[i-l])%mod;
CDQ(mid+,r);
}
CDQ
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=,L=;
for (; n<l; n<<=,L++);
for (int i=; i<n; i++) rev[i]=(rev[i>>]>>)|((i&)<<(L-));
for (int i=; i<n; i++) A[i]=a1[i]>>,B[i]=a1[i]&0x7fff;
for (int i=; i<n; i++) C[i]=a2[i]>>,D[i]=a2[i]&0x7fff;
FFT(A,n,); FFT(B,n,); FFT(C,n,); FFT(D,n,);
for (int i=; 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,-); FFT(s2,n,-); FFT(s3,n,-);
for (int i=; i<n; i++)
a1[i]=((((ll)round(s1[i].x)%mod)<<)%mod+(((ll)round(s2[i].x)%mod)<<)%mod+(ll)round(s3[i].x)%mod)%mod;
}
MTT
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也能解决相同形式的问题。
另外,在快速求解第一、第二类斯特林数时也有应用,见这里。
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)$
第一类斯特林数的快速求解方法之一。注意维护好次数界,这是复杂度的保证。
5.生成函数
很深的部分,一个比较简单的应用是:在做一些生成函数模型较明显的计数问题时,可以通过等比数列求和公式、泰勒展开式将无穷项式化为闭形式,再通过数学知识快速求出某次项系数。
比较基础的应用:
OGF(普通型生成函数)的幂次表示生成序列,EGF(指数型生成函数)的幂次表示生成集合。
EGF就是在OGF的基础上多了卷积是附带的C(n,i),所以可以做多个带标号基本元素的有序排列。
如:轮换计数、树的计数、环套树的计数、无向联通图计数等。
应用:POJ3734
6.其它
我不会的黑科技们
(1)多项式多点求值与快速插值
(2)常系数齐次线性递推
(3)拉格朗日反演
(4)树、图的EGF