快速傅里叶变换

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

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

正常求两个多项式乘积

$$\begin{gathered} A(x)=\sum _{i=0}^n{A}_i{x}^i,B(x)=\sum _{i=0}^n{B}_i{x}^i \\ C(x)=\sum _{i=0}^n\sum _{j=0}^n{A}_i{B}_j{x}^{i+j} \end{gathered}$$

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

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

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

$$\begin{gathered} 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))\} \end{gathered}$$

横坐标相同的点相乘

$$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)$分成奇函数部分与偶函数部分

$$\begin{gathered} P(x)=a_0+a_{1}x+a_2{x}^2+...+a_n{x}^n \\ P(x)=(a_0+a_2{x}^2+a_4{x}^4+...)+x(a_1+a_3{x}^2+a_5{x}^4+...) \\ P(x)=P_e(x^2)+x{P}_o(x^2) \\ P(x_i)=P_e(x_i^2)+x_i{P}_o(x_i^2) \\ P(-x_i)=P_e(x_i^2)-x_i{P}_o(x_i^2) \end{gathered}$$

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

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

那么我们取复数呢？

$$\begin{gathered} i^2=(-i{)}^2=-1 \\ {1}^2=(-1{)}^2=1 \\ \mathrm{因}\mathrm{此}\mathrm{有} \\ i,-i,1,-1:x_1,-x_1,x_2,-x_2 \\ -1,1:x_1^2,x_2^2 \\ 1:x_1^4 \\ \mathrm{可}\mathrm{解}\mathrm{决}n\le 4\mathrm{的}\mathrm{问}\mathrm{题} \end{gathered}$$

那么更大的呢？

我们令$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;
}