快速傅里叶变换FFT
参考文献:《算法导论》
多项式
定义：
多项式加法：
多项式乘法：
多项式的系数表达：
多项式值的点值表达：
插值（将点值还原成多项式）：
已知的将点值转换为多项式系数的方法有高斯消元（$Theta(N^3)$）、拉格朗日插值法（$Theta(N^2)$）等，但通过类似于求点值的方法，可以将复杂度降到$Theta(N*logN)$
点值的作用：
快速乘法：（其中单位复数根接下来会介绍，此处不必纠结）
（此处先抛出，大家不用纠结，以后我还会说）FFT的灵魂在于，
它利用单位复数根的性质将多项式的系数表达和点值表达快速互化，利用点值乘法的特点快速进行乘法运算，最终高效实现了多项式乘法。
（题外话：它将一维（向量）的问题升为二维（矩阵），再降回一维（启发：我们要想得到更好的方法，必须要有更高层的技术））
复数：
在介绍复数之前，首先介绍一些可能会用到的东西
向量
同时具有大小和方向的量
在几何中通常用带有箭头的线段表示
圆的弧度制
等于半径长的圆弧所对的圆心角叫做1弧度的角，用符号rad表示，读作弧度。用弧度作单位来度量角的制度叫做弧度制
公式:
$1^{circ }=frac{pi}{180^{circ}}rad$
$180^{circ }=pi rad$
平行四边形定则
平行四边形定则：$AB+AD=AC$
复数
定义
$a$ , $b$为实数，$i^2=-1$，形如$a+bi$的数叫负数，其中$i$被称为虚数单位，复数域是目前已知最大的域
在复平面中，$x$代表实数，$y$ 轴（除原点外的点）代表虚数，从原点 $(0,0)$ 到 $(a,b)$ 的向量表示复数 $a+bi$
模长：从原点$(0,0)$到点$(a,b)$的距离，即$sqrt{a^2+b^2}$
幅角：假设以逆时针为正方向，从$x$轴正半轴到已知向量的转角的有向角叫做幅角
运算法则
加法：
因为在复平面中，复数可以被表示为向量，因此复数的加法与向量的加法相同，都满足平行四边形定则（就是上面那个）
乘法：
几何定义：
复数相乘，模长相乘，幅角相加
代数定义：
$(a+bi)*(c+di)$
$=ac+adi+bci+bdi^2$
$=ac+adi+bci-bd$
$=(ac-bd)+(bc+ad)i$
代码实现：
struct COM{/*Complex*/ double real,imag; COM(){ real=0; imag=0; } inline COM(double x,double y){ real=x; imag=y; } friend inline COM operator + (COM x,COM y){ return COM(x.real+y.real,x.imag+y.imag); } friend inline COM operator - (COM x,COM y){ return COM(x.real-y.real,x.imag-y.imag); } friend inline COM operator * (COM x,COM y){ return COM(x.real*y.real-x.imag*y.imag,x.real*y.imag+x.imag*y.real); } inline void output(){ printf("(%lf,%lf)",real,imag); } };
单位复数根：
单位复数根是FFT的关键，FFT的关键就是运用单位根的性质加速运算
单位复数根的求法：
因为单位复数根均匀分布在单位圆上，所以可以根据定义直接用三角函数求解:
$∵e^{ui}=cos(u)+isin(u) , w_n=e^{2{pi}i/n}$
$∴w_n=cos(u)+isin(u)$
求$w^x_{y}$代码
inline COM unit(int x,int y){ double ii=2.0*Pi*x/y; return COM(cos(ii),sin(ii)); }
单位复数根的性质
证明：
根据消去引理：
$omega^{frac{n}{2}}_{n}=omega^{frac{n}{2}*frac{2}{n}}_{n*frac{2}{n}}=omega_{2}$
根据单位复数根的定义：
$omega^{n}_{2n}=-1$
∴推论30.4显然成立
式A.5为几何级数（小奥错次相减）：
终于开始讲核心内容了：
DFT
FFT
此处先给出递归FFT的伪代码，随后证明其正确性：
前面部分是分治，for循环中是核心代码，以下是数学证明
插值（从点值到系数）
（此处用了数学中的常用思想（设而不求），定义范德蒙德矩阵但不求解，利用其性质得出等式，帮助解题）
实现
递归实现
这样，我们就基本了解了FFT的思想，此时递归实现就很简单了（可以参考上文中给出的伪代码）：
1 #include<bits/stdc++.h> 2 using namespace std; 3 typedef long long LL; 4 const int INF=1e9+7,MAXN=2e6+50; 5 const double Pi=3.141592653589793238462; 6 struct COM{ 7 double real,imag; 8 COM(){ 9 real=0; 10 imag=0; 11 } 12 inline COM(double x,double y){ 13 real=x; 14 imag=y; 15 } 16 friend inline COM operator + (COM x,COM y){ 17 return COM(x.real+y.real,x.imag+y.imag); 18 } 19 friend inline COM operator - (COM x,COM y){ 20 return COM(x.real-y.real,x.imag-y.imag); 21 } 22 friend inline COM operator * (COM x,COM y){ 23 return COM(x.real*y.real-x.imag*y.imag,x.real*y.imag+x.imag*y.real); 24 } 25 }; 26 void FFT(int lim,COM *a,int opt){ 27 if(lim==1) 28 return; 29 COM a0[lim>>1],a1[lim>>1]; 30 for(int i=0;i<lim;i+=2){ 31 a0[i>>1]=a[i]; 32 a1[i>>1]=a[i+1]; 33 } 34 FFT(lim>>1,a0,opt); 35 FFT(lim>>1,a1,opt); 36 COM unit=COM(cos(2.0*Pi/lim),opt*(sin(2.0*Pi/lim))),w=COM(1,0); 37 for(int i=0;i<(lim>>1);i++){ 38 a[i]=a0[i]+w*a1[i]; 39 a[i+lim/2]=a0[i]-w*a1[i]; 40 w=w*unit; 41 } 42 } 43 int N,M; 44 COM a[MAXN],b[MAXN],c[MAXN]; 45 int limit=1; 46 int main(){ 47 scanf("%d%d",&N,&M); 48 for(int i=0;i<=N;i++) 49 scanf("%lf",&a[i].real); 50 for(int i=0;i<=M;i++) 51 scanf("%lf",&b[i].real); 52 while(limit<=N+M){ 53 limit<<=1; 54 } 55 FFT(limit,a,1); 56 FFT(limit,b,1); 57 for(int i=0;i<=limit;i++){ 58 c[i]=a[i]*b[i]; 59 } 60 FFT(limit,c,-1); 61 for(int i=0;i<=N+M;i++){ 62 printf("%d ",(int)(c[i].real/limit+0.5)); 63 } 64 return 0; 65 }
迭代实现
递归实现之所以不优，有两方面原因（空间，时间）
空间上，如果要简单递归实现，对栈的要求非常高
时间上又有两方面，一是递归本身的效率低下，二是复数的运算非常慢，故我们要选择迭代，以合并一些复数运算
蝴蝶操作
观察之前的伪代码
其中$omega y^1_k$出现了两次，则该值为公用子表达式，可以用临时变量优化成
光是这样，只能很小程度地优化，但如果我们将递归转成递推，同一个$omega$的值就可以用于多次蝴蝶操作
递归转递推
刚才讲到“如果能将递归转成递推”，但这转化并不像普通分治那么简单，因为我们在$FFT$中是按照奇偶分治的，所以要先对递推数组进行处理
如图，我们要将初始数组从${a_0,a_1,a_2,a_3,a_4,a_5,a_6,a_7}$转成${a_0,a_4,a_2,a_6,a_1,a_5,a_3,a_7}$
显然，我们可以在$Theta(NlogN)$内完成，但这样做显然增大了常数，我们需要$Theta(N)$的做法
此处先给出做法：
- 先预处理出一个存储每个数二进制翻转后的数组$rev$（例如$[0,2^3)$中$001$变成$100$）
for(int i=0;i<limit;i++) r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
- 然后对于原数组$A$，如果下标小于$rev$中的值($i<rev_i$)则交换$A_i,A_{rev_i}$
for(int i=0;i<limit;i++) if(i<r[i]) swap(A[i],A[r[i]]);
效果：
$rev$数组：
$000
ightarrow000$
$001
ightarrow100$
$010
ightarrow010$
$011
ightarrow110$
$100
ightarrow001$
$101
ightarrow101$
$110
ightarrow011$
$111
ightarrow111$
过程
$ 0 , 1 , 2 , 3 , 4 , 5 , 6 , 7 leftarrow {000
ightarrow000}$
$ 0 , 4 , 2 , 3 , 1 , 5 , 6 , 7 leftarrow {001
ightarrow100}$
$ 0 , 4 , 2 , 3 , 1 , 5 , 6 , 7 leftarrow {010
ightarrow010}$
$ 0 , 4 , 2 , 6 , 1 , 5 , 3 , 7 leftarrow {011
ightarrow110}$
$ 0 , 4 , 2 , 6 , 1 , 5 , 3 , 7 leftarrow {100
ightarrow001}$
$ 0 , 4 , 2 , 6 , 1 , 5 , 3 , 7 leftarrow {101
ightarrow101}$
$ 0 , 4 , 2 , 6 , 1 , 5 , 3 , 7 leftarrow {110
ightarrow011}$
$ 0 , 4 , 2 , 6 , 1 , 5 , 3 , 7 leftarrow {111
ightarrow111}$
证明：
FFT的递归树就是按照二进制反序分配的，即每个数的位置就是它的原下标的二进制反序。故预处理出二进制反序，每次当前下标小于rev数组则交换（防止重复交换）
至于代码实现，就很简单了：
1 #include<bits/stdc++.h> 2 using namespace std; 3 typedef long long LL; 4 const int INF=1e9+7,MAXN=4e6+50; 5 const double Pi=3.141592653589793238462; 6 struct COM{ 7 double real,imag; 8 COM(){ real=0;imag=0; } 9 inline COM(double x,double y){ real=x;imag=y; } 10 friend inline COM operator + (COM x,COM y){ return COM(x.real+y.real,x.imag+y.imag); } 11 friend inline COM operator - (COM x,COM y){ return COM(x.real-y.real,x.imag-y.imag); } 12 friend inline COM operator * (COM x,COM y){ return COM(x.real*y.real-x.imag*y.imag,x.real*y.imag+x.imag*y.real); } 13 inline void output(){ printf("(%lf,%lf)",real,imag); } 14 }; 15 int rev[MAXN]; 16 void FFT(int lim,COM *a,int opt){ 17 int lg=log2(lim); 18 for(int i=0;i<lim;i++) 19 if(i<rev[i]) 20 swap(a[i],a[rev[i]]); 21 for(int i=0;i<=lg;i++){ 22 int len=1<<i; 23 COM unit=COM(cos(2.0*Pi/len),sin(opt*2.0*Pi/len)); 24 for(int j=0;j<lim;j+=len){ 25 COM w=COM(1,0); 26 for(int k=0;k<(len>>1);k++){ 27 COM t1=a[j+k],t2=w*a[j+k+(len>>1)]; 28 a[j+k]=t1+t2; 29 a[j+k+(len>>1)]=t1-t2; 30 w=w*unit; 31 } 32 } 33 } 34 } 35 int N,M; 36 COM a[MAXN],b[MAXN],c[MAXN]; 37 int limit=1,lgn; 38 int main(){ 39 scanf("%d%d",&N,&M); 40 for(int i=0;i<=N;i++) 41 scanf("%lf",&a[i].real); 42 for(int i=0;i<=M;i++) 43 scanf("%lf",&b[i].real); 44 while(limit<=N+M){ 45 limit<<=1; 46 lgn++; 47 } 48 for(int i=0;i<limit;i++) 49 rev[i]=(rev[i>>1]>>1)|((i&1)<<(lgn-1)); 50 FFT(limit,a,1); 51 FFT(limit,b,1); 52 for(int i=0;i<limit;i++){ 53 c[i]=a[i]*b[i]; 54 } 55 FFT(limit,c,-1); 56 for(int i=0;i<=N+M;i++){ 57 printf("%d ",(int)(c[i].real/limit+0.5)); 58 } 59 return 0; 60 }