多项式乘法是$FFT$的一个很常见的应用，所以就从它开始讲起。
注意：因为要用到复数，这篇文章所有的$i$都表示$\sqrt{-1}$，求和的循环变量使用$j,k$。
对于一个函数$A\left(x\right)=\sum_\limits{j=0}^\limits{n-1}a_jx^j$，我们称它为一个多项式，其中称$a_j$为它的系数。
如果一个多项式$A\left(x\right)$的最高次项的非零系数为$a_k$，则称$k$为它的次数，记作$degree\left(A\right)$。
所有严格大于一个多项式的次数的整数都是它的次数界。
现在定义它们的基本运算。
对于次数界为$n$的两个多项式$$A\left(x\right)=\sum_\limits{j=0}^\limits{n-1}a_jx^j,B\left(x\right)=\sum_\limits{j=0}^\limits{n-1}b_jx^j$$有$$A\left(x\right)+B\left(x\right)=\sum_\limits{j=0}^\limits{n-1}\left(a_j+b_j\right)x^j$$
$$A\left(x\right)B\left(x\right)=\sum_\limits{j=0}^\limits{2n-2}\sum_\limits{k=0}^\limits{j}a_kb_{j-k}x^j$$
多项式$A\left(x\right)B\left(x\right)$的系数向量称为向量$a=\left(a_0,a_1,...,a_{n-1}\right)$和向量$b=\left(b_0,b_1,...,b_{n-1}\right)$的卷积。
显然多项式乘法的朴素算法是$\mathcal{O}\left(n^2\right)$的，而且似乎并不能想出什么更快的算法。
现在换一个角度来看这个问题，考虑如何表示一个多项式。之前我们用合式来表示多项式，这样如果已知它的系数，就确定了这个多项式，所以多项式的系数向量被称为系数表达。使用系数表达可以在$\mathcal{O}\left(n\right)$的时间对于给定的$x$求出对应的值，但是似乎并不能高效地进行乘法运算。
多项式还有另一种表示方法——点值表达。点值表达是一个由2元组组成的集合$$\left\{\left(x_0,y_0\right),\left(x_1,y_1\right),...,\left(x_{n-1},y_{n-1}\right)\right\}$$因为通过点值表达可以列出一个关于$n$个多项式系数的$n$元一次方程，所以只要$x_j$互不相同点值表达就可以唯一确定一个多项式。使用点值表达可以高效地计算多项式乘法，只需对两个多项式用同一组$x_j$求出点值表达，再把对应的$y_j$相乘。
现在就可以很容易地想出一个利用点值表达进行多项式乘法的算法了。
- 加倍多项式的系数界以容纳乘法的计算结果；
- 根据系数表达求出对应的点值表达；
- 用点值表达计算乘法；
- 通过点值表达求出对应的系数表达（这个过程称为插值）。
直接随便找$n$个值带入求值的复杂度是$\mathcal{O}\left(n^2\right)$的，而且插值还很麻烦。
解决方法是带入一些具有特殊性质的值来加速这个过程，单位复数根就是我们要使用的特殊值。
$n$次单位复数根是指满足$\omega^n=1$的复数$\omega$。$n$次单位复数根恰好有$n$个，它们是$\mathrm{e}^{2\pi i/k}\left(k\in\left[0,n\right)\cup\mathbb{Z}\right)$。
$\omega_n=\mathrm{e}^{2\pi i/n}$称为主$n$次单位根，其他的$n$次单位根都是它的幂次。
下面是单位复数根的一些性质。
- 根据定义，$\omega_n^k=\omega_n^{k\bmod n}$；
- 消去引理：对于非负整数$n,k$和正数$d$，有$\omega_{dn}^{dk}=\omega_n^k$；
- 对于正偶数$n$，有$\omega_n^{n/2}=\omega_2=-1$；
- 折半引理：对于正偶数$n$，$n$个$n$次单位根的平方的集合就是$n/2$个$n/2$次单位根的集合。
这些性质比较显然，就不证明了。
现在回到正题，计算一个次数界为$n$的多项式$A$在$n$个$n$次单位根处的值。（实际上是$2n$，因为要在计算前加倍次数界）假设$A$以系数向量$a=\left(a_0,a_1,...,a_{n-1}\right)$的形式给出，则结果为向量$y=\left(y_0,y_1,...,y_{n-1}\right)$，其中$y_j=A\left(\omega_n^j\right)=\sum_\limits{k=0}^\limits{n-1}a_k\omega_n^{jk}$。向量$y$称为向量$a$的离散傅里叶变换（$DFT$），记作$y=DFT_n\left(a\right)$
重点来了，可以使用一种叫做快速傅里叶变换算法，用$\mathcal{O}\left(n\log_2 n\right)$的时间计算$DFT_n\left(a\right)$。以下内容均假设$n$为$2$的整数幂。
首先定义两个新的多项式$$A_0\left(x\right)=\sum_{j=0}^{n/2-1}a_{2j}x^j$$$$A_1\left(x\right)=\sum_{j=0}^{n/2-1}a_{2j+1}x^j$$。
显然有$A\left(x\right)=A_0\left(x^2\right)+xA_1\left(x^2\right)$。
于是问题转化为：
- 求$A_0$和$A_1$在$\left(\omega_n^0\right)^2,\left(\omega_n^1\right)^2,...,\left(\omega_n^{n-1}\right)^2$这$n$个点处的值；
- 由$A_0,A_1$求出$A$。
虽然上面说是$n$个点，但根据折半引理，实际只需计算$n/2$个点即可。所以第一步可以递归求解，现在考虑第二步。
对于$j\in\left[0,n/2\right)\cup\mathbb{Z}$，有
$$y_j=A_0\left(\left(\omega_n^j\right)^2\right)+\omega_n^jA_1\left(\left(\omega_n^j\right)^2\right)$$$$=A_0\left(\omega_{n/2}^j\right)+\omega_n^jA_1\left(\omega_{n/2}^j\right)$$
$$y_{j+n/2}=A_0\left(\left(\omega_n^{j+n/2}\right)^2\right)+\omega_n^{j+n/2}A_1\left(\left(\omega_n^{j+n/2}\right)^2\right)$$$$=A_0\left(\omega_{n/2}^{j+n/2}\right)-\omega_n^jA_1\left(\omega_{n/2}^{j+n/2}\right)$$$$=A_0\left(\omega_{n/2}^j\right)-\omega_n^jA_1\left(\omega_{n/2}^j\right)$$
最后考虑如何进行$DFT$的逆运算，从而由在$n$个单位复数根处的点值表达$y$求出系数表达$a$。（记作$a=DFT^{-1}_{n}\left(y\right)$）
这个记住结论即可。（因为我不会证明
）
$$a_j=\frac{1}{n}\sum_{k=0}^{n-1}y_k\omega_n^{-jk}$$
下面给出代码。
注意到网上很多实现都手写了复数类，实际上C++库本身就有complex类，在<complex>头文件中。
看起来完了对吧，但实际上上面给出的代码基本上没有人用，因为还有一个常数更小、辅助空间更少的非递归（迭代）实现。
考虑如何用迭代的方法来实现$FFT$。
来看一个例子，对于$n=8$的一次$FFT$调用，观察$a$数组的递归结构。
$$\{a_0,a_1,a_2,a_3,a_4,a_5,a_6,a_7\}$$
$$\{a_0,a_2,a_4,a_6\}\{a_1,a_3,a_5,a_7\}$$
$$\{a_0,a_4\}\{a_2,a_6\}\{a_1,a_5\}\{a_3,a_7\}$$
$$\{a_0\}\{a_4\}\{a_2\}\{a_6\}\{a_1\}\{a_5\}\{a_3\}\{a_7\}$$
注意到如果能够计算出最底层递归的顺序，就可以逐层往上进行计算了。
计算方法很简单，就是把下标看作$\log_2 n$位2进制数，对它的每一位逆序排列的结果。例如对于$n=8$，$4=\left(100\right)_2$，那么它在递归底层的顺序就是$\left(001\right)_2=1$。证明很简单，这里就不证了。
下面给出最终的代码。
评论 (0)
