FFT入门

GNN’s blog，原帖

前言：

下了好大决心，花了一个晚上的时间，终于看懂了FFT的理论与流程，然后又花了一个晚上实现代码，做了几道模板题。
FFT理论很深，却又很实用，一个很常用的用法就是加速多项式乘法，使得原来O(n^2)的复杂度减小到O(nlogn)。
下面我来大概讲述一下FFT的基本理论与算法流程，帮助初学者了解FFT，同时也是自己的复习。

---

<1> 介绍：

FFT，全称快速傅里叶变换(fast Fourier transform)，是用来计算离散傅里叶变换(DFT)及其逆变换(IDFT)的快速算法。对于DFT，迷一点讲就是说它把时域信号转化为频域信号（不明白也没关系，并不影响学习算法）。在算法竞赛中，我们经常用FFT来加速卷积，做一些多项式乘法或是高精度乘法之类的。

我主要基于多项式乘法来介绍FFT，下面也都是用多项式乘法来讲的~

---

<2>一些定义：

我们将一个以x为变量的多项式A(x)表示为:

$$A(x) = \sum _{j=0}^{n-1}a_jx^j\\$$
若多项式 $A(x)$的最高次的非零系数为 $a^k$，则称这是一个k次多项式，而它的次数界为>k的任意一个整数（一般取k+1）。
若有两个以x为变量的n次多项式 $A(x)$， $B(x)$，则它们的积可表示为：
$$C(x) = \sum_{j=0}^{2n-2}c_jx^j\\$$
其中
$$c_j = \sum_{k=0}^ja_kb_j-k\\$$
$C(x)$是一个次数界为2n-1的多项式。
另外，下面还大量提到了复数，并用专用符号 $i$ 来表示虚数单位，$i^2=-1$。

---

<3>系数表达与点值表达：

对一个次数界为n的多项式$A(x)$而言，其系数表达是一个由系数组成的向量$a=(a_0,a_1,...,a_{n-1})$。
而它的点值表达是一个由n个点值对组成的集合

$$(x_0,y_0),(x_1,y_1),...,(x_{n-1},y_{n-1})\\$$
使得对$k = 0, 1, ..., n - 1$，所有的$x_k$各不相同，
$$y_k = A(x_k)\\$$
对于点值表达，我们可以看成是把多项式$A(x)$当做一个以$x$为自变量的函数，并把$y$作为因变量，从其图像中取n个点，取其坐标，就得到点值表示。
我们可以看到，求出这个多项式的n个点值复杂度为$O(n^2)$，后面可以看到，如果我们巧妙地选取点$x_k$，就可加速这一过程，使复杂度降至$O(nlogn)$。
我们将求值运算的逆称为插值（就是通过点值表达来求出系数表达），有一个定理（差值多项式的唯一性）：对于任意n个点值对组成的集合，其中所有的$x_k$都不同，那么存在唯一的多项式$A(x)$,满足$y_k = A(x_k), k = 0, 1, ..., n - 1$。证明不再赘述，我们大概可以从两点确定一条直线，三点确定一条抛物线来大概认定这个定理吧。
对于多项式乘法，点值表达是十分方便的。若$C(x) = A(x)B(x)$，则对于任意$x_k$，$C(x_k) = A(x_k)B(x_k)$，即：
若$A(x)$有$(x, y_a)$， $B(x)$有$(x, y_b)$， 则$C(x)$有$(x, y_ay_b)$。

---

<4>FFT的大致流程:

现在我们已经能在线性时间内将两个多项式进行乘法操作了，那么剩下的问题就是如何快速实现多项式系数表达与点值表达的快速转化。
上面提到如果巧妙地选取$x_k$，就可以加速这一过程，而我们要选取的点，就是“单位复数根”，具体介绍将在下节展开。
现在给出FFT算法的大致流程：
*1.加倍次数界：将两个多项式$A(x)$，$B(x)$先补为2n次（高次系数为0）；若不足2的整数次幂位，则再补至2的整数次幂位（以便二分分治处理，后面讲）。
*2.求值：将系数表达式转化为点值表达式（用2n阶的FFT(DFT)，求出多项式在2n次单位复根处的值）。
*3.逐点相乘：将两个多项式点值表达中的$y$逐个相乘，得道积的点值表达。
*4.差值：利用积的点值表达再用FFT(IDFT)求出系数表达。

如此，多项式乘法就做完了，我们从两个多项式的系数表达，得到了积的系数表达。

---

<5>单位复数根:

学FFT，这是极其重要的一点，希望能完全理解其中的变换与结论，最好能手动推算一下，并不难推，大部分是基础的指数运算，但推过与没推过是有很大区别的。

n次单位复根是满足$w^n=1$的复数$w$。n次单位复根恰好有n个，为$e^{2πik/n}$，k=0,1,…n-1。
复数的指数形式的定义：

$$e^{iu}=cos(u)+i *sin(u)\\$$
并且值
$$w_n=e^{2πi/n}\\$$
称为主n次单位根，所有其它n次单位复根都是$w_n$的幂次。
可以得到$w_n^n=w_n^0=1$，意味着$w_n^k=w_n^{kmodn}$，（$w_n^{-1}=w_n^{n-1}$）。
消去引理：对任何整数 n>=0，k>=0，以及d>0，
$$w_{dn}^{dk}=w_n^k\\$$
推论：对任意偶数，有
$$w_n^{n/2}=w_2=-1\\$$
折半引理：如果n>0为偶数，那么n个n次单位复数根的平方的集合就是n/2个n/2次单位复数根的集合。

证明：
根据消去引理，有$(w_n^k)^2=w_{n/2}^k=(w_n^{k+n/2})^2$。（$w_n^{n/2}=-1$）

求和引理：对任意整数n>=1和不能被n整除的非负整数k，有

$$\sum_{j=0}^{n-1}(w_n^k)^j=0\\$$
注：用等差数列的求和公式推导即可。

---

<6>DFT:

我们希望计算次数界为n的多项式$A(x)=\sum_{j=0}^{n-1}a_jx^j$在$w_n^0，w_n^1，w_n^2，...，w_n^{n-1}$处的值。对k=0，1，…，n-1，定义结果$y_k:$

$$y_k=A(w_n^k)=\sum_{j=0}^{n-1}a_jw_n^{kj}$$
我们也记$y=DFT_n(a)$。

---

<7>FFT:

通过FFT，利用单位复数根的特殊性质，我们就可以在$O(nlogn)$的时间内求出$DFT_n(a)$，此时n恰好是2的整数次幂。
FFT采用分治，将原多项式偶数下标与奇数下标的系数分开，得到两个新的n/2次多项式$A^{[0]}(x)$和$A^{[1]}(x):$

$$A^{[0]}(x)=a_0+a_2x+a_4x^2+...+a_{n-2}x^{n/2-1}$$
$$A^{[1]}(x)=a_1+a_3x+a_5x^2+...+a_{n-1}x^{n/2-1}$$
其中，$A^{[0]}(x)包含偶数下标，A^{[1]}(x)包含奇数下标。$
于是有：
$$A(x)=A^{[0]}(x^2)+xA^{[1]}(x^2)$$
$$（注：这样刚好就把奇数项与偶数项的次数插空，还原原式了。）$$
于是，就可以递归分治处理了。
下面给出FFT的伪代码：
$void\ \ FFT(a)${
$\ \ \ \ \ \ n=a.length$
$\ \ \ \ \ \ if(n==1)\ \ return\ a$
$\ \ \ \ \ \ w_n=e^{2πi/n}$
$\ \ \ \ \ \ w=1$
$\ \ \ \ \ \ a^{[0]}=(a_0，a_2，...，a_{n-2})$
$\ \ \ \ \ \ a^{[1]}=(a_1，a_3，...，a_{n-1})$
$\ \ \ \ \ \ y^{[0]}=FFT(a^{[0]})$
$\ \ \ \ \ \ y^{[1]}=FFT(a^{[1]})$
$\ \ \ \ \ \ for\ k=0\ \ to\ \ n/2-1\ $
$\ \ \ \ \ \ \ \ \ \ \ \ y_k=y_k^{[0]}+wy_k^{[1]}$
$\ \ \ \ \ \ \ \ \ \ \ \ y_{k+(n/2)}=y_k^{[0]}-wy_k^{[1]}$
$\ \ \ \ \ \ \ \ \ \ \ \ w=ww_n$
$\ \ \ \ \ \ return\ \ y$
}

注意：$y_k$与$y_k^{[0]}$和$y_k^{[1]}$求点值所用的单位复根是不同的，前者是$w_n^k$，而后者是$w_{n/2}^k$。后者等于前者的平方，对应前面的$"x^2"$。于是，$for$循环中的第一个式子就很容易理解了，对于第二个式子，可以自己推导一下（提示：要用到1:$w_n^{k+(n/2)}=-w_n^k$，2：$w_n^{2k+n}=w_n^{2k}$）。

---

<8>插值:

现在我们已经知道了如何快速的通过系数表达，在单位复根处快速求出点值，那么，接下来就只剩下最后一项任务，那就是在快速时间内完成DFT的逆过程IDFT，在单位复根处插值，得到结果的系数表达。

我们可以把DFT写成矩阵乘积$y=V_na$，其中$V_n$是一个由$w_n$适当幂次填充成的范德蒙德矩阵，矩阵第i行，第J列的元素为$w_n^{ij}$（i，j从0到n-1）。
定理：$V_n^{-1}的(j,k)处元素为w_n^{-kj}/n。$
至于证明，只需把矩阵乘法各点的式子写出来，结合求和引理，就可以发现$V_nV_n^{-1}$主对角线上的元素都为1，其余则都是0（也就是单位矩阵$I$）。
于是可以推出$DFT_n^{-1}(y):$

$$a_j=1/n*\sum_{k=0}^{n-1}y_kw_n^{-kj}$$
其中$j=0，1，...，n-1$。
我们比较第6节中的式子，发现只要做出以下修改：
*1：把a与y互换
*2：用$w_n^{-1}$替换$w_n$
*3：将计算结果每个除以n

最后再来一个定理
卷积定理：对任意两个长度为n的向量a和b，其中n是2的幂，

$$a⊗b=DFT_{2n}^{-1}(DFT_{2n}(a)*DFT_{2n}(b))$$

---

最后，感谢大家的阅读！
如有不足之处，尽可在评论中指出。

本文参考文献：
《算法竞赛入门经典训练指南》——刘汝佳
《算法导论》
（注：本文大部分内容均摘自算法导论，自己加以筛选、注释，第7节中提到的推导过程书中有介绍，若没想出来可以查阅（P534））