12/20/2018 重写。
概述
快速傅里叶变换（Fast Fourier Transform，FFT）是一种可在 $\Theta(n \log n)$ 时间内完成的离散傅里叶变换（Discrete Fourier transform，DFT）算法。
这篇文章将简单地讲解它的工作原理，以及使用 C++ 的代码实现。
基本概念
多项式
在数学中，由若干个单项式相加组成的代数式叫做多项式。多项式中的每个单项式叫做多项式的项，这些单项式中的最高项次数，就是这个多项式的次数。
系数表示
设 $A(x)$ 表示一个 $n-1$ 次多项式，所有项系数组成的 $n$ 维向量 $(a_0,a_1,\cdots,a_{n-1})$ 唯一确定了这个多项式。
$$\begin{aligned} A(x)=& \sum_{i=1}^{n-1}a_ix^i \\ =& a_0+a_1x+a_2x^2\cdots+a_{n-1}x^{n-1}\end{aligned}$$
点值表示
将一组互不相同的 $n$ 个 $x$ 带入多项式，得到的 $n$ 个值。设它们组成的 n 维向量 $(x_0,x_1,\cdots,x_{n-1})$，$(y_0,y_1,\cdots,y_{n-1})$ 唯一确定了这个多项式。
$$\begin{aligned} y_i=& A(x_i) \\ =&\sum_{j=1}^{n-1}a_j x_i^j\end{aligned}$$
多项式乘法
令 $A(x)=\sum_{i=0}^{n-1}a_i x^i$，$B(x)=\sum_{i=0}^{n-1}b_i x^i$，有
$$C(x)=A(x)\times B(x) = \sum_{k=0}^{2n-2}(\sum_{k=i+j}a_i b_j)x^k$$
两个 $n-1$ 次多项式相乘，得到的是一个 $2n-2$ 次多项式，朴素算法的时间复杂度为 $\Theta(n^2)$。
如果用两个多项式在 $2n-1$ 个点取得的点值表示，有
$$C_{y_i}=(\sum_{j=0}^{2n-1}a_j x_i^j)\times(\sum_{j=0}^{2n-1}b_j x_i^j)=A_{y_i}\times B_{y_i}$$
已知 $A$ 和 $B$ 点值表示，求出 $C$ 的点值表示，时间复杂度为 $\Theta(n)$。
复数
复数，为实数的延伸，它使任一多项式方程都有根。复数当中有个“虚数单位” $i$，它是 $-1$ 的一个平方根，即 $i^{2}=-1$。任一复数都可表达为 $x+y_i$，其中 $x$ 及 $y$ 皆为实数，分别称为复数之“实部”和“虚部”。
复平面
数学中，复平面（complex plane）是用水平的实轴与垂直的虚轴建立起来的复数的几何表示。它可视为一个具有特定代数结构笛卡儿平面（实平面），一个复数的实部用沿着 x-轴的位移表示，虚部用沿着 y-轴的位移表示。
复平面的想法提供了一个复数的几何解释。在加法下，它们像向量一样相加；两个复数的乘法在极坐标下的表示最简单——乘积的长度或模长是两个绝对值或模长的乘积，乘积的角度或辐角是两个角度或辐角的和。特别地，用一个模长为 1 的复数相乘即为一个旋转。
单位根
在复平面上，以原点为圆心，1 为半径作得的圆称为单位圆。
数学上，$n$ 次单位根是 $n$ 次幂为 $1$ 的复数。它们位于复平面的单位圆上，构成正 $n$ 边形的顶点，其中一个顶点是 $1$。设所得的幅角为正且最小的单位根为 $\omega_n$，其余的记为 $\omega_n^k$，有
$$\omega_n^k = \cos k \dfrac{2\pi}{n}+i \sin k \dfrac{2\pi}{n}$$
其中 $\omega_n^0 = \omega_n^n = 1$。
性质
$$\begin{aligned} \omega_{2n}^{2k}=&\omega_n^k \\ \omega_n^{k+\frac{n}{2}}=& -\omega_n^k \end{aligned}$$
欧拉公式
$$e^{ix}=\cos x+i\sin x$$
单位根的代码实现
由欧拉公式，求 $n$ 次单位根 $\omega_n$，有
const double PI = acos(-1.0);
complex<double> wn=exp(complex<double>(0,PI/n));其中，$\text{exp}$ 用于计算 $e$ 的幂，$\text{complex}$ 是复数类。
离散傅里叶变换 (Discrete Fourier Transform)
考虑多项式 $A(x)$ 的表示。将 $n$ 次单位根的 $0$ 到 $n – 1$ 次幂带入多项式的系数表示，所得点值向量 $(A(\omega_n^0),A(\omega_n^1),\cdots,A(\omega_n^{n-1}))$ 称为其系数向量 $(a_0,a_1,\cdots,a_{n-1})$ 的离散傅里叶变换。
快速傅里叶变换 (Fast Fourier Transform)
考虑将多项式按照系数下标的奇偶分成两部分
$$A(x)=(a_0 + a_2 x^2 + a_4 x^4 + \cdots) + (a_1 x + a_3 x^3 + a_5 x^5 + \cdots) $$
令
$$\begin{aligned}A_1(x)=&a_0+a_2x+a_4x^2+\cdots \\ A_2(x)=& a_1 + a_3 x + a_5 x^2 + \cdots \end{aligned}$$
则有
$$A(x)=A_1(x^2)+x A_2(x^2)$$
设 $k< \dfrac{n}{2}$，有
$$\begin{aligned} A(\omega_n ^ k) &= A_1(\omega_n ^ {2k}) + \omega_n ^ {k} A_2(\omega_n ^ {2k}) \\ &= A_1(\omega_{\frac{n}{2}} ^ {k}) + \omega_n ^ {k} A_2(\omega_{\frac{n}{2}} ^ {k}) \ \end{aligned} $$
$$\begin{aligned} A(\omega_n ^ {k + \frac{n}{2}}) &= A_1(\omega_n ^ {2k + n}) + \omega_n ^ {k + \frac{n}{2}} A_2(\omega_n ^ {2k + n}) \\ &= A_1(\omega_n ^ {2k} \times \omega_n ^ n)-\omega_n ^ {k} A_2(\omega_n ^ {2k} \times \omega_n ^ n) \\ &= A_1(\omega_n ^ {2k})-\omega_n ^ {k} A_2(\omega_n ^ {2k}) \\ &= A_1(\omega_{\frac{n}{2}} ^ {k})-\omega_n ^ {k} A_2(\omega_{\frac{n}{2}} ^ {k}) \ \end{aligned}$$
也就是说，已知 $A_1(x)$ 和 $A_2(x)$ 在 $\omega_{ \frac{n}{2} } ^ 0,\ \omega_{ \frac{n}{2} } ^ 1,\ \dots,\ \omega_{ \frac{n}{2} } ^ { \frac{n}{2} – 1 }$ 的点值表示，就可以在 $\Theta(n)$ 的时间内求得 $A(x)$ 在 $ \omega_n ^ 0,\ \omega_n ^ 1,\ \dots,\ \omega_n ^ {n – 1} $ 处的点值表示。而关于 $A_1(x)$ 和 $A_2(x)$ 的问题都是相对于 $A(x)$ 规模缩小至 $\dfrac{2}{1}$ 的子问题，分治的边界为常数 $a_i$。这个分治算法在 $\Theta(n\log n)$ 的时间内完成了 DDF。
傅里叶逆变换 (Inverse Fourier Transform)
将点值表示的多项式转化为系数表示，同样可以使用快速傅里叶变换，这个过程叫做傅里叶逆变换。
设 $ (y_0,\ y_1,\ y_2,\ \dots,\ y_{n – 1}) $ 为 $ (a_0,\ a_1,\ a_2,\ \dots,\ a_{n – 1}) $ 的傅里叶变换。考虑另一个向量 $ (c_0,\ c_1,\ c_2,\ \dots,\ c_{n – 1}) $ 满足
$$ c_k = \sum\limits_{i = 0} ^ {n – 1} y_i (\omega_n ^ {-k}) ^ i $$
即多项式 $ B(x) = y_0 + y_1 x + y_2 x ^ 2 + \dots + y_{n – 1} x ^ {n – 1} $ 在 $\omega_n ^ 0,\ \omega_n ^ {-1},\ \omega_n ^ {-2},\ \dots,\ \omega_n ^ {-(n – 1)} $ 处的点值表示。
将上式展开，得
$$\begin{aligned} c_k &= \sum\limits_{i = 0} ^ {n – 1} y_i (\omega_n ^ {-k}) ^ i \\ &= \sum\limits_{i = 0} ^ {n – 1} (\sum\limits_{j = 0} ^ {n – 1} a_j (\omega_n ^ i) ^ j) (\omega_n ^ {-k}) ^ i \\ &= \sum\limits_{i = 0} ^ {n – 1} (\sum\limits_{j = 0} ^ {n – 1} a_j (\omega_n ^ j) ^ i) (\omega_n ^ {-k}) ^ i \\ &= \sum\limits_{i = 0} ^ {n – 1} (\sum\limits_{j = 0} ^ {n – 1} a_j (\omega_n ^ j) ^ i (\omega_n ^ {-k}) ^ i) \\ &= \sum\limits_{i = 0} ^ {n – 1} \sum\limits_{j = 0} ^ {n – 1} a_j (\omega_n ^ {j – k}) ^ i \\ &= \sum\limits_{j = 0} ^ {n – 1} a_j (\sum\limits_{i = 0} ^ {n – 1} (\omega_n ^ {j – k}) ^ i) \end{aligned} $$
对于 $ \sum\limits_{i = 0} ^ {n – 1} (\omega_n ^ {j – k}) ^ i $，当 $k \ne j$ 时，由单位根性质可知（或者用等比数列求和），它的值为 $0$；当 $k = j$ 时，显然为 $n$。
所以，当且仅当 $j=k$ 时，$a_j$ 对 $c_k$ 有 $na_j$ 的贡献，即
$$ \begin{aligned} c_i &= n a_i \\ a_i &= \frac{1}{n} c_i \end{aligned} $$
所以，使用单位根的倒数代替单位根，做一次类似快速傅里叶变换的过程，再将结果每个数除以 $n$，即为傅里叶逆变换的结果。
代码实现
预处理下标
直接分治计算的效率较低，我们使用迭代完成。
具体地，我们以 $n=7$ 为例，寻求分治过程中下标的变化
000 001 010 011 100 101 110
0 1 2 3 4 5 6
0 2 4 6 1 3 5
0 4 2 6 1 5 3
0 4 2 6 1 5 3
000 100 010 110 001 101 011抵达分治边界时的下标，恰好是原下标的二进制翻转。我们可以预处理出下标的对应关系。
蝴蝶操作
考虑合并两个子问题的过程，假设 $A_1(\omega_{ \frac{n}{2} } ^k)$ 和 $A_2(\omega_{ \frac{n}{2} } ^k)$ 分别存在 $a(k)$ 和 $a(\frac{n}{2} + k)$ 中，$A(\omega_n ^ {k})$ 和 $A(\omega_n ^ {k + \frac{n}{2} })$ 将要被存放在 $b(k)$ 和 $b(\frac{n}{2} + k)$ 中，合并的操作可表示为
$$\begin{aligned} b(k) & \leftarrow a(k) + \omega_n ^ k \times a(\frac{n}{2} + k) \\ b(\frac{n}{2} + k) & \leftarrow a(k)-\omega_n ^ k \times a(\frac{n}{2} + k) \ \end{aligned} $$
蝴蝶操作的意义是，引入临时变量，使得合并在原地进行，无需新增一个数组
$$ \begin{aligned} t & \leftarrow \omega_n ^ k \times a(\frac{n}{2} + k) \\ a(\frac{n}{2}+k) & \leftarrow a(k)-t \\ a(k) & \leftarrow a(k)+t \end{aligned} $$
代码
这份代码里，我们利用 $FFT$ 求出了 $A(x)$ 和 $B(x)$ 在 $s$ 次单位根处的点值表示，$\Theta(n)$ 求出卷积的点值表示，然后再次 $FFT$ 求得卷积的系数表示。
#include <bits/stdc++.h>
#define maxl 4294153
using namespace std;
typedef complex<double> cd;
const double PI=acos(-1.0);
cd a[maxl],b[maxl];
int n,m,s,bit,rev[maxl];
int read(){
int x=0,flag=1;char ch=getchar();
while(!isdigit(ch)&&ch!='-') ch=getchar();
if (ch=='-') flag=-1,ch=getchar();
while(isdigit(ch)) x=(x<<3)+(x<<1)+(ch-'0'),ch=getchar();
return x*flag;
}
void init(){
for (bit=1,s=2;(1<<bit)<n+m-1;++bit)
s<<=1;
for (int i=0;i<s;++i)
rev[i]=((rev[i>>1]>>1)|((i&1)<<(bit-1)));
}
// 预处理下标，二进制翻转
void fft(cd *a,int n,int dft){
for (int i=0;i<n;++i)
if (i<rev[i]) swap(a[i],a[rev[i]]);
for (int step=1;step<n;step<<=1){
cd wn=exp(cd(0,dft*PI/step));
// wn 提供了特定的幅角，原理见上文 复数：欧拉公式 和 单位根的代码实现
for (int j=0;j<n;j+=step<<1){
cd wnk(1,0);
for (int k=j;k<j+step;++k,wnk*=wn){
cd t=wnk*a[k+step];
a[k+step]=a[k]-t;
a[k]+=t;
}
}
}
if (dft==-1)
for (int i=0;i<n;++i)
a[i]/=n;
}
int main(){
n=read(),m=read();
n++,m++;init();
for (int i=0;i<n;++i)
a[i].real(read());
for (int i=0;i<m;++i)
b[i].real(read());
fft(a,s,1);
fft(b,s,1);
for (int i=0;i<s;++i)
a[i]*=b[i];
fft(a,s,-1);
for (int i=0;i<n+m-1;++i)
printf("%d ",(int)(a[i].real()+0.5));
putchar('\n');
return 0;
}