FFT
全称为快速傅里叶变换,利用复数运算实现DFT与IDFT的快速算法,常用于多项式系数的卷积加速。
多项式表示方法一般有两种:系数表示和点值表示。
DFT实际上对应系数到点值的变换,IDFT对应点值到系数的变换。FFT采用分治的思想,结合蝴蝶变换化递归为迭代,通过取特殊点——复数域的单位根 ω n k \omega_n^k ωnk (注意因果关系,分治的可行性是基于单位根的选取的)来实现时间复杂度为 O ( n l o g n ) O(nlogn) O(nlogn) 的DFT ,分治决定了 n n n 必须是二次幂。而这种特殊点下的点值表示,可以通过取 ω n − k \omega_n^{-k} ωn−k 来实现IDFT(最后要除以n)。(具体证明就不写啦)
代码:
//n一定是2的幂,flag代表是DFT(1)还是IDFT(-1)
void FFT(cp *a,int n,int flag)
{
for(int i=0;i<n;i++) if(i<r[i]) swap(a[i],a[r[i]]);
for(int k=1;k<n;k<<=1)
{
cp w(1,0),wn(cos(pi/k),sin(pi*flag/k)),x,y;
for(int i=0;i<n;i+=(k<<1),w=cp(1,0))
for(int j=0;j<k;j++,w=w*wn)
{
x=a[i+j]; y=w*a[i+j+k];
a[i+j]=x+y; a[i+j+k]=x-y;
}
}
if(flag==-1) for(int i=0;i<n;i++) a[i].a/=n,a[i].b/=n;
}
顺便来一个复数模板:
struct cp
{
double a,b;
cp(double aa=0,double bb=0) {this->a=aa; this->b=bb;}
cp operator + (const cp x) {return cp(a+x.a,b+x.b);}
cp operator - (const cp x) {return cp(a-x.a,b-x.b);}
cp operator * (const cp x) {return cp(a*x.a-b*x.b,b*x.a+a*x.b);}
cp operator / (const cp x) {return cp((a*x.a+b*x.b)/(x.a*x.a+x.b*x.b),(b*x.a-a*x.b)/(x.a*x.a+x.b*x.b));}
};
我们知道FFT是用来加速多项式乘法的,而多项式乘法本质上是一个卷积的过程:定义两个多项式 A A A 和 B B B,其中 a i a_i ai 表示多项式 A A A 中指数为 i 的项的系数, b i b_i bi 同理,设多项式 Y = A × B Y=A\times B Y=A×B,那么 y k = ∑ i = 0 k a i b k − i y_k=\sum_{i=0}^k a_ib_{k-i} yk=∑i=0kaibk−i ,可以看出这是一个卷积的形式。
要注意多项式乘法不是每一项都符合卷积的规律(比如超过某个多项式最高次幂的项),但是反过来,形如 y k = ∑ i = 0 k a i b k − i y_k=\sum_{i=0}^k a_ib_{k-i} yk=∑i=0kaibk−i 这样的卷积运算就可以通过FFT来加速。贴一个卷积(多项式乘法)的代码:
cp A[maxn],B[maxn];
/*
na,nb分别表示有多少项
返回ans的项数
*/
int conv(double *a,int na,double *b,int nb,double *ans)
{
int l=0;
mem(A,0); mem(B,0);
REP(i,0,na-1) A[i]=cp(a[i],0);
REP(i,0,nb-1) B[i]=cp(b[i],0);
for(nb+=na-1,na=1;na<nb;na<<=1) l++;
REP(i,0,na) r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
FFT(A,na,1); FFT(B,na,1);
REP(i,0,na-1) A[i]=A[i]*B[i];
FFT(A,na,-1);
REP(i,0,nb-1) ans[i]=A[i].a;
return nb;
}
MTT
对于特殊的模数,可以用三模数NTT来求卷积,但是这个不普适,考虑任意模数的FFT:
设 a = k a M + b a , b = k b M + b b a=k_aM+b_a, \ b=k_bM+b_b a=kaM+ba, b=kbM+bb ,那么 a b = k a k b M 2 + ( k a b b + k b b a ) M + b a b b ab=k_ak_bM^2+(k_ab_b+k_bb_a)M+b_ab_b ab=kakbM2+(kabb+kbba)M+babb ,于是可以4次DFT加上3次IDFT求解出 a a a和 b b b的系数卷积,总共的复杂度是7次FFT:
int main()
{
//freopen("input.txt","r",stdin);
int n=read()+1,m=read()+1,p=read(),l=0,M=(1<<15);
REP(i,0,n-1) a[i]=read();
REP(i,0,m-1) b[i]=read();
for(m+=n-1,n=1;n<m;n<<=1) l++;
for(int i=0;i<n;i++)
r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
REP(i,0,n-1) ka[i]=cp(a[i]/M,0),ba[i]=cp(a[i]%M,0);
REP(i,0,n-1) kb[i]=cp(b[i]/M,0),bb[i]=cp(b[i]%M,0);
FFT(ka,n,1);FFT(kb,n,1);FFT(ba,n,1);FFT(bb,n,1);
REP(i,0,n-1) A[i]=ka[i]*kb[i],B[i]=ka[i]*bb[i]+kb[i]*ba[i],C[i]=ba[i]*bb[i];
FFT(A,n,-1);FFT(B,n,-1);FFT(C,n,-1);
REP(i,0,m-1)
{
int x=(LL)floor(A[i].a+0.5)%p;
int y=(LL)floor(B[i].a+0.5)%p;
int z=(LL)floor(C[i].a+0.5)%p;
printf("%d ",(1ll*x*M*M%p+1ll*y*M%p+z)%p);
}
return 0;
}
还存在一种优化方法,可以优化到4次FFT,这个以后再补吧。