0.前言
本文对信息学竞赛中的 FFT 算法及其推导作出一点粗浅的探讨，希望能够帮到大家一点点。
另：作者默认大家已经有了一些高中数学基础，对于课本内基础知识部分可能会涉及较少。建议初中生先学完有关三角函数、复数、向量的部分之后再来学习 FFT 。
1.前置知识
多项式
形如(A(x)=sum_{i=0}^n a_ix^i)的式子称作多项式（Polynomial），多项式可以进行加减乘（除法在此不涉及）等运算。
这种表示方法称作多项式的系数表示（Coefficient representation），多项式还有一种点值表示（Point-value representation）法：
众所周知，平面上的(n)个点((x_1,A(x_1))cdots(x_n,A(x_n)))可以唯一确定一个(n-1)次多项式，于是我们就用这些点来表示一个多项式。
点值表示的优点：将多项式乘法的复杂度由系数表示的(O(n^2))减少到了(O(n))（只需要将点值相乘）。
缺点：系数表示与点值表示的互化还是需要(O(n^2))，不过它是可以优化的（这也正是 FFT 做的事情）。
复数及其运算
首先我们有(sqrt{-1}=i)，并将形如(z=a+bi (a,binmathbb{R}))的数称为复数，并定义其共轭(ar{z}=a-bi)。
其中，(Re z=a)称为实部，(Im z=b)称为虚部。
以上是复数的代数表示，其运算同多项式运算，运算时要注意(i^2=-1)。（例如((a+bi)(c+di)=(ac-bd)+(ad+bc)i)）
复数还可以与平面上的向量一一对应，我们把这种表示叫做复数的几何表示，该平面叫做复平面。
对于一个复数，它的几何表示由模长（向量的长度）和幅角（从(x)轴正半轴逆时针转到该向量的角度）组成。
复数几何表示的加减法运算同向量运算，而乘法的规则是模长相乘、幅角相加。
欧拉公式：(e^{i heta}=cos heta+isin heta)。
单位根及其性质
接下来车速加快，请仔细阅读！
定义方程(z^n=1)的复数根为(n)次单位根，记作(omega_n^k=e^{frac{2pi ki}{n}} (k=0,1,cdots,n-1))。容易发现，这(n)个单位根恰好将单位圆(n)等分。
如图，这是方程(z^3=1)的三个复数根(1,-dfrac 12+dfrac{sqrt{3}}2i,-dfrac 12-dfrac{sqrt{3}}2i)：
如果一个单位根的(0,cdots,n-1)次方恰好构成互不相同的所有单位根，我们就将其称为本原单位根。
显然(omega_n=e^{frac{2pi i}{n}})是一个(n)次本原单位根（为方便，下文中“本原单位根”即特指(omega_n)）。
单位根有一些美妙的性质，我们来看一下。（假设以下的(n)都是正偶数）
( ext{Theorem 1: }omega_n^{2k}=omega_{frac{n}2}^k)
( ext{Proof: })利用复数几何表示的相乘法则证明。
( ext{Theorem 2: }omega_n^{frac{n}2+k}=-omega_n^k)
( ext{Proof: })可以理解为多转了半圈之后恰与原来相反。
习题一：写出方程(z^6=1)的所有复数根。
2.Cooley-Tukey算法
该算法分为两个过程： DFT（Discrete Fourier Transform） 和 IDFT（Inverse Discrete Fourier Transform） ，主要思想是分治。
铺垫了这么多，我们要想，如何利用单位根的一些性质来减少运算量呢？
两位巨佬 Cooley 和 Tukey 想出了一种方法：将每一项按照指数奇偶分类。
我们要计算一个多项式的 DFT ，也就是计算其点值表示，需要将每一个单位根代入计算。
即(A(omega_n^k)=sum_{j=0}^{n-1}a_iomega_n^{kj}=sum_{j=0}^{frac{n}2-1}a_{2j}omega_n^{2kj}+omega_n^ksum_{j=0}^{frac{n}2-1}a_{2j+1}omega_n^{2kj})。
根据刚刚提到的( ext{Theorem 1})，我们可以将其改写成(sum_{j=0}^{frac{n}2-1}a_{2j}omega_{frac{n}2}^{kj}+omega_n^ksum_{j=0}^{frac{n}2-1}a_{2j+1}omega_{frac{n}2}^{kj})。
此时再根据( ext{Theorem 2})，我们可以写出(A(omega_n^k))和(A(omega_n^{k+frac{n}2}))的表示：
我们成功了！现在只需要分别求奇数项和偶数项的 DFT 就能求出来两个值了！运算量减少了一半！（可以证明这个过程的确是(O(nlog n))的）
但是先别高兴得太早，我们还有一步：将运算完的点值表示转化回系数表示（毕竟大多数时候我们要的是乘积多项式的系数）。这时候应该怎么办呢？
Notice:由于证明较复杂且超出了高中的知识范围，建议从此处开始跳过直接看结论。
注意到我们求值的过程实际上就是计算了一个矩阵乘法：
而我们的任务就是求系数矩阵（假设叫做(V)）的逆矩阵。
注意到这个矩阵是一个 Vandermonde 矩阵，对此，我们可以很容易地求出它的逆矩阵。
两个比较繁复的证明：
利用行列式及伴随矩阵
利用拉格朗日插值
然而最终的结果非常优美：设(V={v_{ij}})，则(V^{-1}=frac 1n{v_{ij}^{-1}})。
Notice:现在你可以在此停止了，结论就在下方。
所以，我们只需要将 DFT 中(omega_n^k)换成(omega_n^{-k})（体现为虚部取相反数），再将答案乘以(frac 1n)就行了！
由此，我们就可以用一个函数来实现 DFT 和 IDFT 的过程了。
3.递归实现 FFT
一层层向下递归，求出 FFT ：
void FFT(C *f,int len,int opt){//f是要进行DFT的数组，len是长度，opt是要进行的操作DFT/IDFT
if(len==1)return;
C *f0=f,*f1=f+(len>>1);
for(rg int i=0;i<len;i++)tmp[i]=f[i];
for(rg int i=0;i<(len>>1);i++){//按奇偶性分成两个多项式
f0[i]=tmp[i<<1],f1[i]=tmp[i<<1|1];
}
FFT(f0,len>>1,opt),FFT(f1,len>>1,opt);//分别进行FFT
C wn=C(cos(2*Pi/len),opt*sin(2*Pi/len)),w=C(1,0);//wn是本原单位根，w是wn的幂
for(rg int i=0;i<(len>>1);i++){
tmp[i]=f0[i]+w*f1[i];//按照上面推出来的式子进行计算
tmp[i+(len>>1)]=f0[i]-w*f1[i];
w=w*wn;
}
for(rg int i=0;i<len;i++)f[i]=tmp[i];
}
C++ STL 为我们提供了复数类，但是据说常数巨大，而手写复数类也不难写，所以还是建议大家手写：
struct C {
double re,im;
C(double _re=0.0,double _im=0.0){
re=_re,im=_im;
}
}f[N],g[N],tmp[N];
il C operator + (C a,C b){
return C(a.re+b.re,a.im+b.im);
}
il C operator - (C a,C b){
return C(a.re-b.re,a.im-b.im);
}
il C operator * (C a,C b){
return C(a.re*b.re-a.im*b.im,a.re*b.im+a.im*b.re);
}
而最终实现可以这么写：
FFT(f,len,1),FFT(g,len,1);
...//进行计算
FFT(f,len,-1);
4.迭代实现 FFT
注意到上面的实现方式需要进行大量递归，常数巨大（据说在洛谷会只有77分，但我实测能很宽松地通过），我们考虑如何优化常数。
这时候就需要拿出这张著名的图片了（注意下面最后两个框顺序反了）：
这是我们进行递归的顺序，大家观察一下二进制有什么发现？
可以发现，我们其实就是按照反转二进制排序之后依次分组的！如(0000,1000,0100,1100)倒过来就是(0000,0001,0010,0011)。
所以，我们可以将要进行 FFT 的数组按照反转二进制排序，然后直接迭代即可！
for(rg int i=1;i<len;i++){
r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));//求出i的反转二进制（其中l是二进制位数）
}
而我们不需要递归的 FFT 函数变成了这样：
void FFT(C *f,int len,int opt){
for(rg int i=0;i<len;i++){
if(i<r[i])swap(f[i],f[r[i]]);//排序
}
for(rg int i=2;i<=len;i<<=1){
int mid=i>>1;
C wn=C(cos(2*Pi/i),opt*sin(2*Pi/i));
for(rg int j=0;j<len;j+=i){
C w=C(1,0);
for(rg int k=0;k<mid;k++){
C z=f[j+k+mid]*w;
f[j+k+mid]=f[j+k]-z,f[j+k]=f[j+k]+z;
w=w*wn;
}
}
}
}
用时和内存都有了优化。
另外，关于 FFT 还有一个优化：三次变两次。实现方法是将多项式(g)放到多项式(f)的虚部，这时(f^2)的虚部除以二就是我们要的答案。这个方法也可以省去不少常数。
5.应用
FFT 的最大应用是优化多项式卷积的计算。所谓卷积其实就是多项式乘法，它的定义是：对于两个函数(f(n),g(n))，它们的卷积((f*g)(n)=sum_{i+j=n}f(i)g(j))。（就是个听起来挺nb的名字）
所以，当我们碰到一个恶心的式子时，可以尝试将其化为卷积的形式。
例题：洛谷P3338 [ZJOI2014]力
题目让我们求(E_j=sum_{i=1}^{j-1}dfrac{q_i}{(i-j)^2}-sum_{i=j+1}^{n}dfrac{q_i}{(i-j)^2})的值，我们想办法将它化一化。
为了好看，我们设(f_i=q_i,g_i=dfrac{1}{i^2})。这时候原式变成了(sum_{i=1}^{j-1}f_ig_{j-i}-sum_{i=j+1}^nf_ig_{i-j})。
我们优化一下上下界，令(f_0=g_0=0)，有(E_j=sum_{i=0}^jf_ig_{j-i}-sum_{i=j}^nf_ig_{i-j})。
发现第一项已经是卷积形式，我们来着重推一推第二项(sum_{i=j}^nf_ig_{i-j})。
运用和式中变量代换的技巧，这个式子就等于(sum_{i=0}^{n-j}f_{i+j}g_i)。
这里还有一个常用的技巧，就是引入反转函数(f'_i=f_{n-i})。此时原式变成了(sum_{i=0}^{n-j}f'_{n-j-i}g_i)，再代换一次就是(sum_{i=0}^mf'_{m-i}g_i)。
我们把它化成了卷积形式，现在可以用 FFT 做了！
核心代码（代码里h[i]就是我们说的(f'_i)）：
int main(){
Read(n);
for(rg int i=1;i<=n;i++){
cin>>f[i].re;h[n-i].re=f[i].re;
g[i].re=1.0/(1.0*i*i);
}
while(len<n+n+1)len<<=1,l++;
for(rg int i=1;i<len;i++){
r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
}
FFT(f,len,1),FFT(g,len,1),FFT(h,len,1);
for(rg int i=0;i<=len;i++)f[i]=f[i]*g[i],h[i]=h[i]*g[i];
FFT(f,len,-1),FFT(h,len,-1);
for(rg int i=1;i<=n;i++){
printf("%.3lf
",(f[i].re-h[n-i].re)/len);
}
return 0;
}
6.总结
FFT 是多项式的基础算法，也是许多新手学多项式要过的第一关。理解了此处的知识，才能为之后的学习打下良好基础。
习题二：洛谷P3803 【模板】多项式乘法（FFT）
习题三：洛谷P1919 【模板】A*B Problem升级版（FFT快速傅里叶）