??其实我不知道我是否真的理解了FFT,但是我会用FFT优化多项式乘法了QAQ。。
(以下大多摘自算导
前置知识
1. 多项式
??在一个代数域F上,关于变量x的多项式定义为形式和形式表示的函数
A(x)=∑j=0n?1ajxj,其中a0…an?1为多项式各项的系数
2. 多项式的次数界
??若多项式有非零系数的最高次项为xk,则称k为该多项式的次数,任何严格大于k的整数都是这个多项式的次数界。
3. 多项式的表示
(1)系数表示法
??对于一个次数界为n的多项式A(x)来说,其系数表示法可以看做是一个列向量a=(a0,a1,…,an?1)。系数表示法对于某些多项式的计算很方便,如对多项式A(x)在给定点x0的求值运算就是计算A(x0)的值。如果使用霍纳法则,则求值运算的运行时间为O(n),
A(x0)=a0+x0(a1+x0(a2+...+x0(an?2+x0(an?1))...))
??另外,加法运算的时间复杂度是O(n),暴力进行乘法运算的时间复杂度是O(n2)。
(2)点值表示法
?? 有n个点(x0,A(x0),),(x1,A(x1)),..,(xn?1,A(n?1)),当所有xk各不相同时,这n个点可以唯一表示一个次数界为n的多项式,但是一个次数界为n的多项式可以有多个点值表示。通俗一点说,已知一个多项式函数的n个函数值可以唯一确定这个函数,而知道这个函数可以知道不止n个函数值。
??已知点值表示求系数表示称为插值,用拉格朗日插值法可以做到O(n2)。
??若次数界为n的多项式A(x),B(x)的n个点的xk是对应相同的,点值表示法的加法操作时间复杂度是O(n),只要把对应A/B(xk)相加即可,若A(x),B(x)都已知2n个点,那把A/B(xk)相乘同样可以在O(n)时间内完成乘法运算。。
进入正题。。
1. 单位复根
??n次单位复根是满足wn=1的复数w,有n个,他们均匀的分布在以复
平面的原点为圆心的单位圆上,为e2πki/n,k=0,1,…,n?1,复数幂定义为eui=cos(u)+sin(u)i 。wn=e2πi/n称为主n次单位根,所有的其他n次单位根都是wn的幂。
??因为有wnn=w0n=1 ,所以有wjnwkn=w(j+k)modnn 。(性质1
重要性质
相消引理
?? 对任何整数n≥0,k≥0,d>0,有 wdkdn=wkn,利用定义不难证明
?? 推论:对任意偶数n>0,有wn/2n=w2=1 。
折半引理
?? 如果n>0为偶数,n个n次单位复根的平方等于n/2个n/2次单位复根。利用相消引理和性质1可证。
求和引理
??对于任意整数n≥1和不能被n整除的非零整数k,有 ∑n?1j=0wkjn=0。
证明:
?? 原式=(wkn)n?1wkn?1=(wnn)k?1wkn?1=1?1wkn?1=0 ,只要k不整除n就可以保证分母不为0。
2. DFT
??A(x)是一个次数界为n的多项式,不失一般性地假定n是2的幂,因为
次数界总是可以增大的。有列向量y=(y0,y1,…,yn?1),其中yk=A(wkn),则称y是系数向量a的离散傅里叶变换,也写作y=DFTn(a)。
3. 快速傅里叶变换(FFT)
??FFT是一种可以在O(nlogn)时间内计算出DFTn(a)的算法,主要思想
是分治。
??定义A0(x)=a0+a2x+a4x2+…+an?2xn/2?1,包含了A(x)
所有偶数下标的系数,A1(x)=a1+a3x+a5x2+…+an?1xn/2?1,包含了A(x)所有奇数下标的系数。易得A(x)=A0(x2)+xA1(x2),所以我们可以先求A0和A1的DFT,然后再组合起来。组合的时候有
A(wk+n/2n)=A0(w2k+nn)+wk+n/2nA1(w2k+nn)=A0(wkn)+wknw2A1(wkn)A0(wkn)?wknA1(wkn)
所以用n2个yk就可以推出全部。根据折半引理,问题的规模缩小一半,每次组合的时间复杂度是O(n),所以总时间复杂度是O(nlogn)。
?? 求出DFTn(a)后,要对单位复根进行插值,将点值表示转化为系数表示,有y=Vna,其中
Vn=????????????1111...11wnw2nw3n...wn?1n1w2nw4nw6n...w2n1w3nw6nw9n...w3n..................1wn?1nw2(n?1)nw3(n?1)n...w(n?1)(n?1)n????????????
是范德蒙特矩阵。
a=Vn?1?y,考虑Vn?1在(j,k)处的数Vj,k,根据V?1nVn=In(n阶单位矩阵)有
∑t=0n?1Vj,twktn={1(j=k)0(j!=k)
定理:Vj,k=w?kjn/n
证明:∑n?1t=0Vj,twktn=1n∑n?1t=0wt(k?j)n,当j=k时,上式=1,否则由求和引理可得上式为0,得证。
??那么有 aj=1n∑n?1k=0ykw?kjn,因此在求逆DFT的时候可以类似于求DFT,只要把y和a角色互换,然后让 wn=w?1n,再做FFT即可。
??写的时候发现非递归要比递归快很多。。
递归:
#include<cstdio>
#include<iostream>
#include<cmath>
#include<memory.h>
#define N 400010
using namespace std;
const double pi=acos(-1);
struct complex{
double x,i;
complex(){}
complex(double x,double i):x(x),i(i){}
complex operator+(complex a) {return complex(x+a.x,i+a.i);}
complex operator-(complex a) {return complex(x-a.x,i-a.i);}
complex operator*(complex a) {return complex(x*a.x-i*a.i,x*a.i+i*a.x);}
}a[N],b[N];
int n,m,i,nn;
void fft(complex *a,int n,int t)
{
if (n==1) return;
complex a0[n>>1],a1[n>>1];
for (int i=0;i<n;i+=2) a0[i>>1]=a[i],a1[i>>1]=a[i+1];
fft(a0,n>>1,t);fft(a1,n>>1,t);
complex wn(cos(2*pi/n),t*sin(2*pi/n)),w(1,0);
for (int i=0;i<(n>>1);i++,w=w*wn) a[i]=a0[i]+w*a1[i],a[i+(n>>1)]=a0[i]-w*a1[i];
}
int main()
{
scanf("%d%d",&n,&m);
memset(a,0,sizeof(a));memset(b,0,sizeof(b));
for (i=0;i<=n;i++) scanf("%lf",&a[i].x);
for (i=0;i<=m;i++) scanf("%lf",&b[i].x);
nn=1;while (nn<=n+m) nn<<=1;
fft(a,nn,1);fft(b,nn,1);
for (i=0;i<=nn;i++) a[i]=a[i]*b[i];
fft(a,nn,-1);
for (i=0;i<=n+m;i++) printf("%d ",(int)(a[i].x/nn+0.5));
}
非递归:
#include<cstdio>
#include<iostream>
#include<cmath>
#include<memory.h>
#define N 400010
using namespace std;
const double pi=acos(-1);
struct complex{
double x,i;
complex(){}
complex(double x,double i):x(x),i(i){}
complex operator+(complex a) {return complex(x+a.x,i+a.i);}
complex operator-(complex a) {return complex(x-a.x,i-a.i);}
complex operator*(complex a) {return complex(x*a.x-i*a.i,x*a.i+i*a.x);}
}a[N],b[N];
int n,m,i,nn,len,rev[N];
void fft(complex *a,int n,int t)
{
for (int i=0;i<n;i++) if (i<rev[i]) swap(a[i],a[rev[i]]);
for (int j=1;j<n;j<<=1)
{
complex wn(cos(2*pi/(j<<1)),t*sin(2*pi/(j<<1)));
for (int i=0;i<n;i+=(j<<1))
{
complex w(1,0),t0,t1;
for (int k=0;k<j;k++,w=w*wn) t0=a[i+k],t1=w*a[i+j+k],a[i+k]=t0+t1,a[i+j+k]=t0-t1;
}
}
}
int main()
{
freopen("FFT.in","r",stdin);
scanf("%d%d",&n,&m);
memset(a,0,sizeof(a));memset(b,0,sizeof(b));
for (i=0;i<=n;i++) scanf("%lf",&a[i].x);
for (i=0;i<=m;i++) scanf("%lf",&b[i].x);
nn=1;len=0;while (nn<=n+m) nn<<=1,len++;
rev[0]=0;
for (i=1;i<nn;i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<(len-1));
fft(a,nn,1);fft(b,nn,1);
for (i=0;i<=nn;i++) a[i]=a[i]*b[i];
fft(a,nn,-1);
for (i=0;i<=n+m;i++) printf("%d ",(int)(a[i].x/nn+0.5));
}
??非递归这里有一个翻转的函数,意在让所有数按合并时候的二叉树的叶子节点的顺序排列,不难发现翻转过来就是它的新位置。。