快速傅里叶变换

快速计算多项式相乘系数 FFT快速傅里叶变换

快速傅里叶变换(FFT)——有史以来最巧妙的算法？

正常求两个多项式乘积

\[A(x)=\sum_{i=0}^{n}{A_ix^i},B(x)=\sum_{i=0}^{n}{B_ix^i} \\ C(x)=\sum_{i=0}^{n}\sum_{j=0}^nA_iB_jx^{i+j} \]

复杂度为\(O(n^2)\)

分别设n+1个在A与B上的点

根据高斯消元法，n+1个不同的点一定能确定一个唯一的n次多项式

\[A(x)=\{(x_0,A(x_0)),(x_1,A(x_1)),(x_2,A(x_2))...(x_n,A(x_n))\}\\ B(x)=\{(x_0,B(x_0)),(x_1,B(x_1)),(x_2,B(x_2))...(x_n,B(x_n))\}\\ \]

横坐标相同的点相乘

\[C(x)=\{(x_0,A(x_0)B(x_0)),(x_1,A(x_1)B(x_1))...(x_n,A(x_n)B(x_n))\} \]

那么就有两个问题：

1.怎么快速得到n+1个点

2.怎么将点值表达法变回系数

如果\(A(x)\)是一个偶函数，那么\(A(x)=A(-x)\)，因此带入一个值也就得到了两个点

同理如果\(A(x)\)是一个奇函数，那么有\(A(x)=-A(-x)\)，也可以带一个点得到两个点

因此可以将\(A(x)\)分成奇函数部分与偶函数部分

\[P(x)=a_0+a_1x+a_2x^2+...+a_nx^n\\ P(x)=(a_0+a_2x^2+a_4x^4+...)+x(a_1+a_3x^2+a_5x^4+...)\\ P(x)=P_e(x^2)+xP_o(x^2)\\ P(x_i)=P_e(x_i^2)+x_iP_o(x_i^2)\\ P(-x_i)=P_e(x_i^2)-x_iP_o(x_i^2) \]

然后如果将\(P_e(x),P_o(x)\)继续分解

但我们的\(P_e(x),P_o(x)\)是通过\(P(x)=P_e(x^2)+xP_o(x^2)\)得来的，因此取相反数也无法让它们一一符号相反

那么我们取复数呢？

\[i^2=(-i)^2=-1\\ 1^2=(-1)^2=1\\ 因此有\\ i,-i,1,-1:\ x_1,-x_1,x_2,-x_2\\ -1,1:x_1^2,x_2^2\\ 1:x_1^4 \\可解决n\leq4的问题 \]

那么更大的呢？

我们令\(w=e^{\frac{2\pi i}{n}}=cos(\frac{2\pi}{n})+isin(\frac{2\pi}{n})\)，我们可以得到它的性质正好可以解决该问题

反过来呢？

将\(w=e^{-\frac{2\pi i}{n}}\)，然后给每个数都乘上\(\frac{1}{n}\)即可

//初始化每个位置最终到达的位置 
void init(int k)
{
    int len=1<<k;
	for(int i=0;i<len;i++)
	rev[i]=(rev[i>>1]>>1)|((i&1)<<(k-1));
}
//a表示要操作的系数，n表示序列长度
//若flag为1，则表示FFT，为-1则为IFFT(需要求倒数）
void fft(cp *a,int n,int flag){ 
    for(int i=0;i<n;i++)
	{
	 //i小于rev[i]时才交换，防止同一个元素交换两次，回到它原来的位置。 
	  if(i<rev[i])swap(a[i],a[rev[i]]);
	}
	for(int h=1;h<n;h*=2)//h是准备合并序列的长度的二分之一
	{
		cp wn=exp(cp(0,flag*pie/h));//求单位根w_n^1 
		for(int j=0;j<n;j+=h*2)//j表示合并到了哪一位
		{
			cp w(1,0);
			for(int k=j;k<j+h;k++)//只扫左半部分，得到右半部分的答案
			{
			    cp x=a[k];
			    cp y=w*a[k+h];
		        a[k]=x+y;  //这两步是蝴蝶变换 
		        a[k+h]=x-y;
		        w*=wn; //求w_n^k 
			}
		 }
	 }
	 //判断是否是FFT还是IFFT 
	 if(flag==-1)
	 for(int i=0;i<n;i++)
     a[i]/=n;
}