快速傅里叶变换(Fast Fourier Transform, FFT)

介绍

计算两个多项式的乘积，按照定义直接计算需要 $O(n^2)$ 的时间，但通过快速傅里叶变换，就可以在 $O(n \log n)$ 的时间计算。

多项式

形如 $\displaystyle \sum_{i=0}^{n-1}a_ix^i$ 的代数表达式叫做多项式，其中 $x$ 是多项式的未知数， $a_{0..n-1}$ 是多项式的系数， $n-1$ 是多项式的次数(实际上次数要求 $a_{n-1}$ 非零，下面为了方便允许 $a_{n-1}=0$ )。

对于一个多项式 $A(x)$ ，一般的表示方式是系数表示，就是用多项式的每一项的系数来表示这个多项式；多项式 $A(x)=\displaystyle \sum_{i=0}^{n-1} a_ix^i$ 的系数表示就是向量 $a=(a_0,a_1,\cdots,a_{n-1})$ 。而多项式 $A(x)$ 的点值表示就是选 $n$ 个不同的数对多项式进行求值，得到的点值对组成的集合，即 $\{(x_i,y_i):0 \le i <n \}$ ，使得对于 $0 \le i <n$ ， $x_i$ 互不相同， $y_i=A(x_i)$ 。

一个多项式可以有很多种不同的点值表示，但每一种点值表示都唯一对应着一个多项式。这是因为根据 $n$ 个点值对我们可以以 $n$ 个系数为变量列出 $n$ 条方程，然后把系数解出来。至于说方程为什么一定有解，可以利用范德蒙矩阵的性质来证明，可参考算法导论。

由多项式的系数表示得到点值表示，称为求值；由多项式的点值表示得到系数表示，称为插值。

多项式的乘法

对于两个 $n$ 次多项式 $\displaystyle A(x)=\sum_{i=0}^{n}a_ix^i,B(x)=\sum_{i=0}^{n}b_ix^i$ ，他们的乘积 $C(x)=\displaystyle \sum_{i=0}^{2n}c_ix^i$ ，其中 $c_i=\displaystyle \sum_{j}a_jb_{i-j}(0 \le j \le n,0 \le i-j \le n)$ 。

在系数表示下，多项式乘法的直接计算需要 $O(n^2)$ 的时间；但在点值表示下，只需要 $O(n)$ 的时间。因为如果 $C(x)=A(x)B(x)$ ，则对于 $x_k$ ， $C(x_k)=A(x_k)B(x_k)$ 。只需要把 $A$ 和 $B$ 的点值表达逐点相乘即可。不过，需要注意到一个问题，如果 $A$ 和 $B$ 的次数为 $n$ ，那么 $C$ 的次数为 $2n$ ，因此我们需要 $2n+1$ 个点值对才能确定 $C$ 。只需要对 $A$ 和 $B$ 都求出 $2n+1$ 个点值对即可。

考虑用点值表示下的 $O(n)$ 乘法来加速系数表示下的乘法。那么问题的关键就在于如何快速的把一个多项式从系数表示转化成点值表示(求值)，以及转回来(插值)。这时，求值点的选择就非常重要。如果选择"复数单位根"作为求值点，就可以在 $O(n\log n)$ 时间内完成求值和插值。

复数单位根

复数的 $n$ 次单位根是指满足 $w^n=1$ 的复数 $w$ 。根据复数乘法相当于模长相乘，辐角相加可知， $n$ 次单位根的模长一定是 $1$ ，辐角的 $n$ 倍一定是 $0$ 。因此， $n$ 次单位根也就是 $e^{{2 \pi k i}/{n} },k=0,1,\cdots,n-1$ 。 (根据欧拉公式 $e^{\theta i}=\cos \theta +i \sin \theta$ ， $e^{\theta i}$ 就相当于单位圆上辐角为 $\theta$ 的点) 记 $w_n=e^{2 \pi i /n}$ ，那么 $n$ 次单位根就是 $w_n^{0},w_n^{1},\cdots,w_n^{n-1}$ 。

求值(DFT)

现在有一个 $n-1$ 次多项式 $A(x)=\displaystyle \sum_{i=0}^{n-1}a_ix^i$ ，我们希望计算它在 $n$ 个 $n$ 次复数单位根处的点值，即 $y_k=A(w_n^k)=\displaystyle\sum_{j=0}^{n-1}w_n^{kj},k=0,1,\cdots,n-1$ 。这里，向量 $y=(y_0,y_1,\cdots,n-1)$ 称为向量 $a=(a_0,a_1,\cdots,a_{n-1})$ 的 离散傅里叶变换 (DFT)，也记作 $y=DFT_n(a)$ 。

快速傅里叶变换 $(FFT)$ 是一种计算 $DFT_n(a)$ 的算法。这里要求 $n$ 是 $2$ 的次幂。在多项式乘法中，我们可以在高次项补 $0$ 做到这一点。考虑将 $A(x)$ 分成偶次项和奇次项。定义 $A^{[0]}(x)=a_0+a_2x+a_4x^2+\cdots+a_{n-2}x^{n/2-1}$ $A^{[1]}(x)=a_1+a_3x+a_5x^2+\cdots+a_{n-1}x^{n/2-1}$ 于是有 $A(x)=A^{[0]}(x^2)+xA^{[1]}(x^2)$ 。所以，要求 $A(x)$ 在 $w_n^0,w_n^1,\cdots,w_n^{n-1}$ 处的值，就可以分别求出 $A^{[0]}(x)$ 和 $A^{[1]}(x)$ 在 $(w_n^0)^2,(w_n^1)^2,\cdots,(w_n^{n-1})^2$ 处的值然后进行合并。由于 $(w_n^k)^2=w_n^{2k}=w_{n/2}^{k}$ ， $(w_n^{k+n/2})^2=w_{n/2}^{k+n/2}=w_{n/2}^{k}$ ，因此 $n$ 个 $n$ 次单位根的平方恰好就是所有 $\displaystyle \frac{n}{2}$ 次单位根、每个出现 $2$ 次。至此，问题转化为求两个次数为 $\displaystyle \frac{n}{2}-1$ 的多项式 $A^{[0]}(x)$ 和 $A^{[1]}(x)$ 在所有 $\displaystyle \frac{n}{2}$ 次单位根处的值，即两个 $DFT_{n/2}$ 的问题，可以递归处理。

递归 $\log n$ 层，每层 $O(n)$ ，所以时间复杂度 $O(n \log n)$ 。

插值(IDFT)

现在我们希望将一个多项式从点值表示转化回系数表示，这个变换称作 $IDFT$ ，也就是 $DFT$ 的逆。

根据 $n$ 个点值对我们可以以 $n$ 个系数为变量列出 $n$ 条方程，问题相当于是解这个方程组：

写成矩阵的形式就是：

记 $(1)$ 中的系数矩阵(就是最左边那个)为 $V$ 。考虑求出 $V$ 的逆矩阵 $V^{-1}$ ，即满足 $V^{-1} \times V=I_n$ ，其中 $I_n$ 为 $n \times n$ 的单位矩阵，然后在 $(1)$ 两边同时左乘 $V^{-1}$ 来解出答案。

考虑下面这个矩阵 $d_{i,j}=w_n^{-ij}$

设他们相乘的结果为 $E=D \times V$ $e_{i,j}=\displaystyle \sum_{k=0}^{n-1}d_{i,k}v_{k,j} =\sum_{k=0}^{n-1}w_n^{-ik}w_n^{kj} = \sum_{k=0}^{n-1}w_n^{k(j-i)}$

当 $i=j$ 时， $e_{i,j}=n$ ；

当 $i \not= j$ 时， $e_{i,j}=\displaystyle\sum_{k=0}^{n-1}(w_n^{j-i})^{k}=\frac{1-(w_n^{j-i})^n}{1-w_n^{j-i}}=0$ 。

因此 $E=nI_n$ ，所以 $V^{-1}=\displaystyle\frac{1}{n}D$ ，将 $(1)$ 两边同时左乘 $V^{-1}$ ，得：

因此， $IDFT$ 就相当于把 $w_n$ 换成 $w_n^{-1}$ ，求 $DFT$ ，再把结果除掉 $n$ 。由于 $w_n^{-1}$ 和 $w_n$ 具有类似的性质，所以这里求 $DFT$ 的方法和上面是一样的。

实现

模板题：uoj.ac/problem/34

直接实现是递归的，下面直接给出代码：uoj.ac/submission/…

但非递归实现代码更短，常数更小，是一般的写法。

考虑算法的过程，以 $n=8$ 为例，

初始系数序列为：

分治一次后：

分治两次后：

分治三次后：

考虑最终那个序列，相当于将所有下标以二进制从低到高的第 $1$ 位为第一关键字，第 $2$ 位为第二关键字， $\cdots$ ，第 $k$ 位为第 $k$ 关键字排序( $2^{k}=n$ )。那么我们将下标二进制的 $k$ 位翻转一下，就是第 $1$ 位变成第 $k$ 位，第 $2$ 位变成第 $k-1$ 位， $\cdots$ ，第 $k$ 位变成第 $1$ 位，那么原来的排序就相当于按下标从小到大排序，即不变。因此，只需要将 $x_i$ 跟 $x_{rev(i)}$ 交换一下，就能得到最终的序列，其中 $rev(i)$ 表示将 $i$ 二进制的 $k$ 位翻转得到的数。