快速傅里叶变换（FFT）

引入：
多项式：
$A(x)=a_0+a_{1}x+a_2{x}^2+...+a_{n-1}{x}^{n-1}$
$B(x)=b_0+b_{1}x+b_2{x}^2+...+b_{m-1}{x}^{m-1}$
在计算多项式的乘法，如$A(x)\times B(x)$时，往往需要$O(n^2)$的复杂度来计算：

（平方的竖式计算）
这是因为我们用的系数表示法来存储的多项式（即把多项式的系数按顺序存在数组里）

然而，有一种更好的表示方法——点值表示法：给一个多项式带入不同的x值，得到不同的y值（如同一个函数图像，从上面取n个点的坐标）
即 $y_k=A(x_k)(所有x_k不相同，k=0,1,2,...,n-1)$
这样，我们可以固定一串$x_k$，用数组存储一串$y_k$来表示多项式。那么两个多项式相乘$A(x_k)\times B(x_k)=A_{yk}\times B_{yk}$，只需要将数组对应乘起来就可以了，$O(n)$！

然而又有问题了，好像没有人也没有题读多项式用点值表示法吧。。。我们需要把系数表示法转换为点值表示法，乘了过后，又转换回去。
如何转换？（数学老师：把点带进函数，解方程。。。好像又变回$O(n^2)$了，那刚刚那个$O(n)$乘法并没有卵用）

又有一个新的圣器——快速傅里叶变换
选用一个神奇的叫单位复根的东西来作为$x_k$，用分治的思想$O(nlogn)$完成转换。

单位复根
$n$次单位复根即n次方为1的数，${\omega }^n=1$的复数$\omega$，而n次单位复根恰好有n个，分别是$e^{2\pi ik/n}$，$(k=0,1,...,n-1)$。
又∵$e^{i\theta }=cos\theta +isin\theta$
∴n个n次单位复根就会像这样呈现↓

${\omega }_a^b$就是一个半径为1的圆，转$\frac{b}{a}$圈的位置。
这样，单位复根就有了一些神奇的数学性质：

• ${\omega }_{dn}^{dk}={\omega }_n^k$，证明很简单，转到$\frac{dk}{dn}$圈就等于$\frac{k}{n}$。
• $n>0$且n为偶数，则n个n次单位复根的平方等于n/2个n/2次单位复根。
  证明：
  设$k\le \frac{n}{2}$
  $({\omega }_n^k{)}^2={\omega }_n^{2k}={\omega }_{2n/2}^{2k}={\omega }_{n/2}^k$
  $({\omega }_n^{n/2+k}{)}^2={\omega }_n^{n+2k}={\omega }_n^n{\omega }_n^{2k}=({\omega }_n^k{)}^2$（同上） （∵${\omega }_n^n=1$）
• 对于任意整数$n\ge 1$和$k$不为$n$的倍数，即${\omega }_n^k̸=1$，则${\Sigma }_j^{n-1}({\omega }_n^k{)}^j=0$
  证明：
  这为等比数列求和，利用公式↓
  ${\Sigma }_j^{n-1}({\omega }_n^k{)}^j=\frac{({\omega }_n^k{)}^n-1}{{\omega }_n^k-1}=\frac{1^k-1}{{\omega }_n^k-1}=0$
• ${\omega }_n^{k+n/2}=-{\omega }_n^k$，相当于再$\frac{k}{n}$处再转半圈，正好转到对面

FFT之Cooley-Tukey算法：分治思想。
为保证分治成功，需要将n扩大到最近的2的幂。
我们要计算$A(x)=a_0+a_{1}x+a_2{x}^2+...+a_{n-1}{x}^{n-1}$
的点值表示法数组y
$y_i=A({\omega }_n^i)$
设
$A^{[0]}(x)=a_0+a_{2}x+a_4{x}^2+...+a_{n-2}{x}^{n/2-1}$
$A^{[1]}(x)=a_1+a_{3}x+a_5{x}^2+...+a_{n-1}{x}^{n/2-1}$
$\therefore A(x)=A^{[0]}(x^2)+x{A}^{[1]}(x^2)$
所以我们可以先分治求解数组$y^{[0]}$与$y^{[1]}$。
$y_k^{[0]}=A^{[0]}({\omega }_n^{2k})=A^{[0]}({\omega }_{n/2}^k)$
$y_k^{[1]}=A^{[1]}({\omega }_n^{2k})=A^{[1]}({\omega }_{n/2}^k)$
$(k=0,1,2...\frac{n}{2}-1)$
规模正好缩小一半。
然后我们需要用$y^{[0]}$与$y^{[1]}$来合成数组$y$。
因为$k<\frac{n}{2}$，所以$y$数组需要一半一半的求。
前半部分：
$y_k=A({\omega }_n^k)=A^{[0]}({\omega }_{n/2}^k)+{\omega }_n^k{A}^{[1]}({\omega }_{n/2}^k)=y_k^{[0]}+{\omega }_n^k{y}_k^{[1]}$
后半部分：
$y_{k+n/2}=A({\omega }_n^{k+n/2})=A^{[0]}({\omega }_{n/2}^{k+n/2})+{\omega }_n^{k+n/2}{A}^{[1]}({\omega }_{n/2}^{k+n/2})=A^{[0]}({\omega }_{n/2}^k)-{\omega }_n^k{A}^{[1]}({\omega }_{n/2}^k)=y_k^{[0]}-{\omega }_n^k{y}_k^{[1]}$
分治直到长度为1，此时$y=A(...)=a_0$
然后使用上面的方法合并。

逆FFT：将点值表示法转换回系数表示法。
只需将点值表示法数组y当作系数，将${\omega }_n^{-k}$代入计算，将结果除以n，即得到$a_k$。
$a_k=\frac{y_0+y_1{\omega }_n^{-k}+y_2{\omega }_n^{-2k}+...+y_{n-1}{\omega }_n^{-(n-1)k}}{n}$
又∵$y_k=A({\omega }_n^k)=a_0+a_1{\omega }_n^k+a_2{\omega }_n^{2k}+...+a_{n-1}{\omega }_n^{(n-1)k}$
∴在上面分子中所有y分离出$a_i$单独计算：
①$i=k$，∴可分离出：
$a_i+a_i{\omega }_n^i{\omega }_n^{-k}+a_i{\omega }_n^{2i}{\omega }_n^{-2k}+...+a_i{\omega }_n^{(n-1)i}{\omega }_n^{-(n-1)k}=a_i+a_i+a_i+...+a_i=n{a}_k$
②$i̸=k$，每项均可分离出：
$a_i+a_i{\omega }_n^i{\omega }_n^{-k}+a_i{\omega }_n^{2i}{\omega }_n^{-2k}+...+a_i{\omega }_n^{(n-1)i}{\omega }_n^{-(n-1)k}$
$=a_i{\omega }_n^{(i-k)0}+a_i{\omega }_n^{(i-k)}+a_i{\omega }_n^{(i-k)2}+...+a_i{\omega }_n^{(i-k)(n-1)}$

$=a_i{\Sigma }_j^{n-1}({\omega }_n^{(i-k)}{)}^j=0$

综上所诉，∴分子$y_0+y_1{\omega }_n^{-k}+y_2{\omega }_n^{-2k}+...+y_{n-1}{\omega }_n^{-(n-1)k}=n{a}_k$
∴$\frac{y_0+y_1{\omega }_n^{-k}+y_2{\omega }_n^{-2k}+...+y_{n-1}{\omega }_n^{-(n-1)k}}{n}=a_k$

逆FFT做法就是：将点值表示法数组y当作系数，将${\omega }_n^{-k}$代入计算，将结果除以n，即得到$a_k$。

一个利用FFT做多项式乘积的程序。
代码：（递归版本）

#include<cstdio>
#include<cmath>
#include<cstring>
#include<complex>
using namespace std;
const int MAXN=1005;
typedef complex<double> cpxd;
cpxd out[MAXN];
void FFT(cpxd *A,cpxd *out,int n,int flag,int step=1)
/*
A是原数列
out是返回的点值表示法
n是长度
flag=1 FFT；flag=-1 逆FFT
step表示相邻两个数的距离。
因为FFT不停的提出奇偶项递归，则此时递归的数列在原数列A中，相邻两个数的距离为step
比如1,2,3,4,5,6,7,8 -> 1,3,5,7 -> 1,5  (这两个数在原数列中距离为4)

*/
{
	if(n<1)return;
	if(n==1)
	{out[0]=A[0];return;}
	FFT(A,out,n/2,flag,step*2);
	FFT(A+step,out+n/2,n/2,flag,step*2);
	for(int i=0;i<n/2;i++)
	{
		cpxd temp=exp(cpxd(0,flag*acos(-1)*2.0*i/n))*out[n/2+i];
		//exp那一坨生成单位复根
		out[n/2+i]=out[i]-temp;
		out[i]=out[i]+temp;
	}
	if(flag==-1&&step==1)
		for(int i=0;i<n;i++)
			out[i]/=n;
}
void mul(cpxd *A,cpxd *B,cpxd *C,int n)
{
	int len=1;
	n<<=1;
	while(len<n)len<<=1;
	FFT(A,out,len,1);
	memcpy(A,out,sizeof(cpxd)*len);
	FFT(B,out,len,1);
	memcpy(B,out,sizeof(cpxd)*len);
	for(int i=0;i<len;i++)
		C[i]=A[i]*B[i];
	FFT(C,out,len,-1);
	memcpy(C,out,sizeof(cpxd)*len);
}
int main()
{
	int n;
	scanf("%d",&n);
	cpxd a[MAXN],b[MAXN],c[MAXN];
	double x;
	for(int i=1;i<=n;i++)
	{
		scanf("%lf",&x);
		a[i]=cpxd(x,0.0);
	}
	for(int i=1;i<=n;i++)
	{
		scanf("%lf",&x);
		b[i]=cpxd(x,0.0);
	}
	mul(a+1,b+1,c+1,n);
	for(int i=1;i<2*n-1;i++)
		printf("%lf ",c[i].real());
	printf("%lf\n",c[2*n-1].real());
	return 0;
}
/*
3
1 2 3
3 2 1
*/
/*
3.000000 8.000000 14.000000 8.000000 3.000000
*/

再来改为非递归版本，速度更快
看看递归版本如何迭代的：

如果我们把初始的数组排成最后的样子，就不需要递归了。
观察最后的数字的二进制
0 4 2 6 1 5 3 7
000 100 010 110 001 101 011 111
正好是把数字倒过来按顺序排列

代码：（非递归版本）

#include<cstdio>
#include<cmath>
#include<complex>
using namespace std;
const int MAXN=1005;
typedef complex<double> cpxd;
void FFT(cpxd *A,int n,int flag)
{
	for(int i=0,j=0;i<n;i++)
	{
		if(j>i)swap(A[i],A[j]);//排序
		int k=n;
		while(j&(k>>=1))
			j^=k;
		j|=k;
		//把二进制数字倒过来后，生成下一个数
	}
	double pi=flag*acos(0)*2.0;
	for(int i=1;i<n;i<<=1)//每次递归的长度的一半
	{
		double tmp=pi/i;//因为i为长度的一半，所以÷i就已经乘以了2
		for(int j=0;j<i;j++)
		{
			cpxd w=exp(cpxd(0,tmp*j));//生成单位复根
			for(int l=j,r;l<n;l+=(i<<1))//计算每一组长度为2i的数列，一次合并
			{
				r=l+i;
				cpxd temp=w*A[r];
				A[r]=A[l]-temp;
				A[l]=A[l]+temp;
			}
		}
	}
	if(flag==-1)
		for(int i=0;i<n;i++)
			A[i]/=n;
}
void mul(cpxd *A,cpxd *B,cpxd *C,int n)
{
	int len=1;
	n<<=1;
	while(len<n)len<<=1;
	FFT(A,len,1);
	FFT(B,len,1);
	for(int i=0;i<len;i++)
		C[i]=A[i]*B[i];
	FFT(C,len,-1);
}
int main()
{
	int n;
	scanf("%d",&n);
	cpxd a[MAXN],b[MAXN],c[MAXN];
	double x;
	for(int i=1;i<=n;i++)
	{
		scanf("%lf",&x);
		a[i]=cpxd(x,0.0);
	}
	for(int i=1;i<=n;i++)
	{
		scanf("%lf",&x);
		b[i]=cpxd(x,0.0);
	}
	mul(a+1,b+1,c+1,n);
	for(int i=1;i<2*n-1;i++)
		printf("%lf ",c[i].real());
	printf("%lf\n",c[2*n-1].real());
	return 0;
}
/*
3
1玩
2 3
3 2 1
*/
/*
3.000000 8.000000 14.000000 8.000000 3.000000
*/

完