引入
可能有不少OIer都知道FFT这个神奇的算法, 通过一系列玄学的变化就可以在 $O(nlog(n))$ 的总时间复杂度内计算出两个向量的卷积(或者多项式乘法/高精度乘法), 而代码量却非常小. 博主一年半前曾经因COGS的一道叫做"神秘的常数 $\pi$"的题目而去学习过FFT, 但是基本就是照着板子打打完并不知道自己在写些什么鬼畜的东西OwO
不过...博主这几天突然照着算法导论自己看了一遍发现自己似乎突然意识到了什么OwO然后就打了一道板子题还1A了OwO再加上午考试差点AK以及日更频率即将不保于是就有了这篇博文233
博主在写这篇文字的过程中也发现了不少自己之前忽略的FFT细节, 感觉对FFT似乎有了更深刻的理解, 希望想学习FFT的读者能认真看完这篇文章并仔细体味FFT的优雅性质.
本文可能并不会介绍关于使用FFT进行信号处理的相关信息, 主要介绍在OI中的应用即快速卷积.
Rush了好久终于写出来了QAQ
那么现在, 就让我们变玄学为科学了解 $FFT$ 背后的工作原理与它的优雅性质
前置知识: 多项式的表示
多项式的系数表达
相信大家都知道对于一个次数界为 $n$ 的多项式 $A(x)$ 可以表示为如下形式:
\[A(x) = \sum _{i=0} ^{n-1} a_i x^i\]
这时, 多项式 $A(x)$ 的系数所组成的向量 $\vec{a}=(a_0,a_1,...a_{n-1})$ 是这个多项式的系数表达. (实际上是列向量不是行向量OwO)
使用系数表达来表示多项式可以进行一些方便的运算, 比如对其进行求值, 这时我们可以采用秦九昭算法来在 $O(n)$ 时间复杂度内对多项式进行求值.
\[A(x_0) = a_0+x_0(a_1+x_0(a_2+...+x_0(a_{n-2}+x_0 a_{n-1})...))\]
但是当我们想把两个采用系数表达的多项式乘起来的时候, 恭喜你, 问题来了: 一个多项式中的每个系数都要与另一个多项式中的每个系数相乘. 然后时间复杂度就变成了鬼畜的 $\Theta(n^2)$ ...
也就是对于两个多项式的系数表示 $\vec{a}$ 和 $\vec{b}$ 来说, 两个多项式乘积的系数表达 $\vec{c}$ 中的值的求值公式如下:
\[c_i=\sum_{j=0}^i a_jb_{i-j}\]
然而多项式乘法和卷积在很多地方都需要快速的实现, 而对于多项式乘法来说, 另一种表示方法似乎更为合适:
多项式的点值表达
首先我们对点值表达进行定义:
一个次数界为 $n$ 的多项式 $A(x)$ 的点值表达为一个由 $n$ 个点值对所组成的集合:
\[\{(x_0,y_0),(x_1,y_1),(x_2,y_2),...,(x_{n-1},y_{n-1})\}\]
使得对于 $k=0,1,2,...,n-1$ , 所有的 $x_k$ 各不相同
(我们可以先暂且把 $A(x)$ 看作一个函数)
其中 $y_k=A(x_k)$ , 也就是说选取 $n$ 个不同的值分别代入求值之后产生 $n$ 个值, 就会得到 $n$ 个未知数的值与多项式值的点对. 这 $n$ 个点对就是多项式的一个点值表达.
不难看出点值表达并不像多项式表达一样具有唯一性, 因为 $x_k$ 是没有限制条件的.
对于一个系数表达的多项式, 显然求出它的一个点值表达非常方便, 只需取 $n$ 个不同的 $x$ 值代入并求出对应的点值即可. 计算出这些点值需要 $O(n^2)$ 的时间.
....
等等...? 怎么还是 $O(n^2)$ ? 别急, 很快我们就会发现, 如果我们像某些丧病出题人一样精心选择数据, 我们就可以在 $O(nlog(n))$ 的时间复杂度内完成这个运算.
而从多项式的点值表达转换为唯一的系数表达就没有那么显然了. 首先我们介绍两个概念:
求值与插值
从一个多项式的系数表达确定其点值表达的过程称为求值(似乎很显然的样子? 毕竟我们求点值表达的过程就是取了 $n$ 个 $x$ 的值然后扔进了多项式求了 $n$ 个值出来233), 而求值运算的逆运算(也就是从一个多项式的点值表达确定多项式的系数表达)被称为插值. 而插值多项式的唯一性定理告诉我们只有多项式的次数界等于已知的点值对的个数, 插值过程才是明确的(能得到一个确定的多项式表达) , 联想之前的矩阵与线性方程组和增广矩阵可以得到如下证明:
定理(插值多项式的唯一性): 对于任意 $n$ 个点值对组成的集合 $\{(x_0,y_0),(x_1,y_1),...,(x_{n-1},y_{n-1}\}$ ,且其中所有 $x_k$ 不同, 则存在唯一的次数界为 $n$ 的多项式 $A(x)$ , 满足 $y_k=A(x_k),k=0,1,...,n-1$ .
证明 证明主要是根据某个矩阵存在逆矩阵. 多项式系数表达等价于下面的矩阵方程:
\[
\begin{bmatrix}
1 & x_0 & x_0^2 & \dots & x_0^{n-1} \\
1 & x_1 & x_1^2 & \dots & x_1^{n-1} \\
\vdots & \vdots & \vdots & \ddots & \vdots \\
1 & x_{n-1} & x_{n-1}^2 & \dots & x_{n-1}^{n-1}
\end{bmatrix}
\begin{bmatrix} a_0 \\ a_1 \\ \vdots \\ a_{n-1} \end{bmatrix}
=
\begin{bmatrix} y_0 \\ y_1 \\ \vdots \\ y_{n-1} \end{bmatrix}
\]左边的矩阵表示为 $V(x_0,x_1,...,x_{n-1})$ , 称为范德蒙矩阵. 而这个矩阵的行列式为:
\[\prod _{0\leq j < k \leq n-1} (x_k-x_j)\]
所以根据定理 "$n\times n$ 矩阵是奇异的, 当且仅当该矩阵的行列式为 $0$", 如果 $x_k$ 皆不相同, 则这个矩阵是可逆的. 因此给定点值表达 $\vec y$ , 则能确定唯一的系数表达 $\vec{a}$ 使:
\[\vec{a}=V(x_0,x_1,...,x_{n-1})^{-1}\vec{y}\]
$\blacksquare$
上面的证明过程实际上也给出了插值的一种算法, 即求出范德蒙矩阵的逆矩阵. 我们可以在 $O(n^3)$ 时间复杂度内求出逆矩阵所以相应地可以在 $O(n^3)$ 时间复杂度内计算出点值表达所对应的系数表达.
然而这特么比求值还慢...=_=
所幸我们还有一种基于 $n$ 个点的插值算分, 基于拉格朗日公式:
\[A(x)=\sum_{k=0}^{n-1} y_k \frac{\prod _{j\neq k} (x-x_j)} {\prod_{j\neq k} (x_k-x_j)}\]
基于这个公式我们可以在 $O(n^2)$ 时间复杂度内计算出点值表达所对应的系数表达.
然后插值勉强追上了求值的时间复杂度.
那么问题来了, 计算这个点值表达有啥用?
点值表达的运算
点值表达在很多多项式相关操作上比系数表达要方便很多, 比如对于加法, 如果 $C(x)=A(x)+B(x)$ ,则对于任意点值 $x_k$ , 满足 $C(x_k)=A(x_k)+b(x_k)$ . 而正确性是显而易见的.所以对于两个点值表达的多项式进行加法操作, 时间复杂度为 $O(n)$.
与加法类似, 多项式乘法在点值表达的情况下也变得更加方便: 如果 $C(x)=A(x)B(x)$, 则对于任意点 $x_k$ , $C(x_k)=A(x_k)B(x_k)$ . 也就是说点值表达常常可以大大简化多项式间的运算!
然而问题又来了: 多项式 $C$ 的次数界显然是 $A$ 的次数界与 $B$ 的次数界的和. 所以如果我们要插值并获得唯一的系数表达式就需要更多的点值对. 这时我们就有必要对 $A$ 与 $B$ 的点值表达进行"扩展", 至少扩展到 $C$ 的次数界.
但是我们似乎发现了一个严重的问题: 虽然点值表达中计算乘法速度比系数表达快很多, 但是到目前为止我们在两种表达方式之间转换的算法都是 $O(n^2)$ 的, 也就是说转换之后时间复杂度相较于朴素的多项式乘法实现不但没有提升, 甚至还增加了不少运算量.
对于这个问题, 我们可以采取一些方法: 因为我们在计算点值表达时可以采用任何点作为求值点, 通过精心挑选求值点就可以做到 $O(nlog(n))$ 的时间复杂度进行两种表示法之间的转换! 也就是说我们可以采取如下策略来加速卷积:
1.加倍次数界: 在多项式 $A$ 和多项式 $B$ 中加入若干值为 $0$ 的高阶系数使其达到多项式 $C$ 所需的次数界, 这一过程为 $O(n)$ 时间复杂度.
2.求值: 在 $O(nlog(n))$ 时间复杂度内计算出多项式 $A$ 和多项式 $B$ 的点值表达.
3.逐点相乘: 在 $O(n)$ 时间复杂度内计算出多项式 $C$ 的点值表达.
4.插值: 在$O(nlog(n))$ 时间复杂度内计算出多项式 $C$ 的系数表达
也就是说我们可以在 $O(nlog(n))$ 的总时间复杂度内计算出多项式 $A$ 与多项式 $B$ 的乘积.
而这个关键的求值与插值算法又是什么呢? 没错, 就是快速傅里叶变换(Fast Fourier Transform, FFT)
DFT与FFT
我们刚刚说到, 当我们仔细选择求值点的话就可以加速求值过程, 而我们所选择的值就是一种具有特殊性质的数:n次单位复数根.
n次单位复数根
n次单位复数根是满足 $\omega ^n=1$的复数 $\omega$ . 根据复数乘法运算的几何意义(幅角相加, 模相乘)和同余, 我们可以得出 $n$ 次单位复数根恰好有 $n$ 个且均匀分布在以原点为圆心, 一单位长度为半径的圆周上. 如下两图所示:
所以 . 根据复数在指数上的定义 $e^{iu}=cos(u) + i\:sin(u)$ , 我们可以得到对于 $k=0,1,...,n-1$ , $n$ 次单位根是 $e^{2\pi i k/n}$
其中, $\omega _n=e^{2\pi i /n}$ 为主 $n$ 次单位根. 因为所有其它单位根实际上都是它的整次幂.
实际上 $n$ 次单位根形成了一个乘法群(感兴趣请自行查阅群伦相关资料, 在这里写的话估计又要扯很久).
我们注意到对于一个 $n$ 次单位根 $\omega_n^k$, 它在极坐标系中的幅角就是 $\frac{k}{n}$ 个周角. 因为它们均匀分布在单位圆上. 而对于分数, 我们是可以进行约分的, 这也就引出了对后续FFT十分重要的结论:
消去引理, 折半引理与求和引理
直接引用算法导论中给出的证明:
引理(消去引理): 对任意整数 $n \geq 0, k\geq 0 $, 以及 $d\geq 0$, 有:
\[\omega_{dn}^{dk}=\omega_n^k\]
证明: 由单位复数根的定义式可直接推出引理, 因为
\[\omega_{dn}^{dk}=(e^{2\pi i / dn})^{dk}=(e^{2\pi i /}n)^k= \omega_n^k\]
$\blacksquare$
十分优雅而精简的证明, 正是运用了单位复数根的特殊性质. 证明过程中使用了幂的幂中指数相乘的运算法则, 将 $d$ 乘进括号内, 与括号内在指数分母上的另一个 $d$ 相抵消.
根据消去引理, 我们还可以推出另一个引理:
引理(折半引理): 如果 $n>0$ 为偶数, 那么 $n$ 个 $n$ 次单位复数根的平方的集合就是 $n/2$ 个$n/2$ 次单位复数根的集合.
证明: 根据消去引理, 对任意非负整数 $k$ , 我们有 $(\omega_n^k)^2=\omega_{n/2}^k$ . 注意, 如果对所有 $n$ 次单位复数根进行平方, 那么获得每个 $n/2$ 次单位根正好两次. 因为:
\[(\omega_n^{k+n/2})^2=\omega_n^{2k+n}=\omega_n^{2k}\omega_n^n=\omega_n^{2k}=(\omega_n^k)^2\]
因此, \(\omega_n^{k+n/2}\) 与 \(\omega_n^k\) 平方相同.
$\blacksquare$
从消去引理推出的折半引理是后面FFT中使用的分治策略的正确性的关键, 因为折半引理保证了递归子问题的规模只是递归调用前的一半, 从而保证 $O(nlog(n))$ 时间复杂度.
另一个重要的引理是求和引理, 是后文中使用FFT快速插值的基础.
引理(求和引理): 对于任意整数 $n\geq 1$ 与不能被 $n$ 整除的非负整数 $k$ , 有:
\[\sum_{j=0}^{n-1}(\omega_n^k)^j=0\]
证明: 根据几何级数的求和公式:
\[\sum_{k=0}^nx^k=\frac{x^{n+1}-1}{x-1}\]
我们可以得到如下推导:
\[\sum_{j=0}^{n-1}(\omega_n^k)^j=\frac{(\omega_n^k)^n-1}{\omega_n^k-1}=\frac{(\omega_n^n)^k-1}{\omega_n^k-1}=\frac{(1)^k-1}{\omega_n^k-1}=0\]
定理中要求 $k$ 不可被 $n$ 整除, 而仅当 $k$ 被 $n$ 整除时有 $\omega_n^k=1$, 所以保证分母不为 $0$ .
$\blacksquare$
介绍了这么多前置知识与理论基础, 我们该开始今天的重头戏了:
离散傅里叶变换 DFT
首先我们还是给出一个DFT的定义:
还记得我们需要FFT去做什么吗? 我们想要快速从多项式的系数表达计算出它的点值表达, 而上一节我们说过, 我们要用 $n$ 个 $n$ 次单位复数根处进行这一求值过程.
对于 $A$ 的系数表达 $\vec{a}= (a_0,a_1,...,a_{n-1})$ 和 $k=0,1,...,n-1$ , 定义结果 $y_k$:
\[y_k=A(\omega_n^k)=\sum_{i=0}^{n-1}a_i\omega^{ki}_n\]
则 $\vec{y}=(y_0,y_1,...,y_{n-1})$ 就是系数向量 $\vec{a}$ 的离散傅里叶变换(DFT). 简记为 $\vec{y}=DFT_n(\vec{a})$
然而问题一个接一个: 从这个定义式来看, 我们计算出一个多项式 $A(x)$ 的离散傅里叶变换依然需要 $O(n^2)$ 的时间复杂度, 依然没有任何提升的样子?
好了, 真正的主角该登场了:
快速傅里叶变换 FFT
如果我们使用FFT利用单位复数根的特殊性质,我们可以在 $O(nlog(n))$ 时间复杂度内计算出 $DFT_n(\vec a)$ . 如果直接按定义式计算的话依然是和求值相同的时间复杂度.
在下文中为了分治方便, 我们将假设 $n$ 是 $2$ 的整次幂, 虽然对于非整次幂的 $n$ 的算法已经出现, 但是平常 OI 中使用的时候一般采取补一大坨值为 $0$ 的高次项系数强行补到 $2$ 的整次幂23333
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$ 中下标为偶数的系数, 而 $A^{[1]}(x)$ 则包含了 $A$ 中下标为奇数的系数. 所以实际上我们根据下标可以得到如下结论:
\[A(x)=A^{[0]}(x^2)+A^{[1]}(x^2)\]
因为次数界折半了所以代入的是 $x^2$ , 然后对于奇数下标系数还要乘上一个 $x$ 作为补偿.
于是我们就把求 $A(x)$ 在 $\omega_n^0,\omega_n^1,...,\omega_n^{n-1}$ 处的值转化成了这样:
1.求多项式 $A^{[0]}(x)$ 和 多项式 $A^{[1]}(x)$ 在下列点上的取值:
\[(\omega_n^0)^2,(\omega_n^1)^2,...(\omega_n^{n-1})^2\]
2.根据刚刚推导出的将两个按奇偶拆开的多项式的值合并的公式合并结果.
但是根据折半引理, 这 $n$ 个要代入子多项式中求值的点显然并不是 $n$ 个不同值, 而是 $n/2$ 个 $n/2$ 次单位复数根组成. 所以我们递归来对次数界为 $n/2$ 的多项式 $A^{[0]}(x)$ $A^{[1]}(x)$ 在 $n/2$ 个 $n/2$ 次单位复数根处求值. 然我们成功地缩小了了子问题的规模. 也就是说我们将 $DFT_n$ 的计算划分为两个 $DFT_{n/2}$ 来计算. 这样我们就可以计算出 $n$ 个元素组成的向量 $\vec a$ 的DFT了.
然后让我们结合FFT的C++实现来理解一下:
1 void FFT(Complex* a,int len){
2 if(len==1)
3 return;
4 Complex* a0=new Complex[len/2];
5 Complex* a1=new Complex[len/2];
6 for(int i=0;i<len;i+=2){
7 a0[i/2]=a[i];
8 a1[i/2]=a[i+1];
9 }
10 FFT(a0,len/2);
11 FFT(a1,len/2);
12 Complex wn(cos(2*Pi/len),sin(2*Pi/len));
13 Complex w(1,0);
14 for(int i=0;i<(len/2);i++){
15 a[i]=a0[i]+w*a1[i];
16 a[i+len/2]=a0[i]-w*a1[i];
17 w=w*wn;
18 }
19 delete[] a0;
20 delete[] a1;
21 }
让我们逐行解读一下这个板子:
第2,3行是递归边界, 显然一个元素的 $DFT$ 就是它本身, 因为 $\omega_1^0=1$, 所以
\[y_0=a_0\omega_1^0=a_0\]
第4~9行将多项式的系数向量按奇偶拆成两部分.
第12,13和第17行结合起来用来更新 $\omega$ 的值, 12行计算主 $n$ 次单位复数根, 13行将变量 $w$ 初始化为 $\omega_n^0$ 也就是 $1$ , 17行使在组合两个子多项式的答案时可以避免重新计算 $\omega_n$ 的幂, 而是采用递推方式最大程度利用计算中间结果(OI中的常用技巧, 对于主 $n$ 次单位复数根还可以打个表来节约计算三角函数时的巨大耗时).
第10,11行将拆出的两个多项式递归处理. 我们设多项式 $A^{[0]}(x)$ 的 $DFT$ 结果为 $\vec{y^{[0]}}$ , $A^{[1]}(x)$ 的 $DFT$ 结果为 $\vec{y^{[1]}}$, 则根据定义, 对于 $k=0,1,...,n/2-1$, 有:
\[y_k^{[0]}=A^{[0]}(\omega_{n/2}^k)\]
\[y_k^{[1]}=A^{[1]}(\omega_{n/2}^k)\]
然后根据消去引理, 我们可以得到 $\omega_{n/2}^k=\omega_n^{2k}$, 所以
\[y_k^{[0]}=A^{[0]}(\omega_n^{2k})\]
\[y_k^{[1]}=A^{[1]}(\omega_n^{2k})\]
设整个多项式的 $DFT$ 结果为 $\vec y$, 则第15,16行用如下推导来将递归计算两个多项式的计算结果综合为一个多项式的结果:
第15行, 对于 $y_0,y_1,...,y_{n/2-1}$, 我们设 $k=0,1,...,n/2-1$ :
\begin{aligned} y_k&=y^{[0]}_k+\omega_n^ky_k^{[1]} \\ &=A^{[0]}(\omega_n^{2k})+\omega_n^kA^{[1]}(\omega_n^{2k}) \\ &= A(\omega_n^k) \end{aligned}
上面的推导中, 第一行是我们在代码中进行的运算, 第二行将 $y_k^{[0]}$ 和 $y_k^{[1]}$ 按定义代入, 第三行基于我们刚刚推导出的组合两个子多项式 $DFT$ 结果的式子.
第16行, 对于 $y_{n/2}, y_{n/2+1}, ..., y_{n-1}$ , 我们同样设 $k=0,1,...,n/2-1$ :
\begin{aligned} y_{k+(n/2)} &= y^{[0]}_{k+(n/2)} - \omega_n^ky_{k+(n/2)}^{[1]} \\ &= y_k^{[0]} + \omega_n^{k+(n/2)} y_k^{[1]}y_k^{[1]} \\ &= A^{[0]}(\omega_n^{2k}) +\omega_n^{k+(n/2)} A^{[1]}(\omega_n^{2k}) \\ &= A^{[0]}(\omega_n^{2k+n}) + \omega_n^{k+(n/2)} A^{[1]}(\omega_n^{2k+n}) \\ &= A(\omega_n^{k+(n/2)}) \end{aligned}
上面的推导中, 第一行同样是我们在代码中进行的运算, 第二行从 $\omega_n^{k+(n/2)}=-\omega_n^k$ 推导出, 第三行将 $y_k^{[0]}$ 和 $y_k^{[1]}$ 按定义代入, 第四行从 $\omega_n^n=1$ 推导出, 最后一行基于推导出的组合递归结果用的等式.
根据这两段推导, 我们在代码中所作的计算能够得到正确的$\vec y$.
第19,20行用于释放我们开出的临时内存空间.
上面的几段推导都印证了 $n$ 次单位复数根与DFT的优美对称性, 所以我们得以减少重复计算来得到更优的时间复杂度.
其中我们把 $\omega_n^k$ 称为旋转因子.
上述的算法实际上对应于下面的递归式:
\[T(n)=2T(n/2)+\Theta(n)=\Theta(nlog(n))\]
然而还有一个问题
我们到现在为止一直在讨论快速求值, 并且成功地通过FFT将时间复杂度降到了 $O(nlog(n))$, 然而我们还需要插值来将点值表达转化为系数表达. 这时我们就得...
通过 FFT 在单位复数根处插值
根据我们的四步快速卷积计算方案(倍次/求值/点积/插值), 我们还剩下最后一步就可以完成这一切了.
首先, $DFT$ 的计算过程按照定义可以表示成一个矩阵, 根据上文关于范德蒙德矩阵的等式 , 我们可以将 $DFT$ 写成矩阵乘积 $\vec y=V_n \vec a$, 其中 $V_n$ 是用 $\omega_n$ 的一定幂次填充的矩阵, 这个矩阵大概长这样:
\[
\begin{bmatrix} y_0 \\ y_1 \\ y_2 \\ y_3 \\ \vdots \\ y_{n-1} \end{bmatrix}=
\begin{bmatrix}
1 & 1 & 1 & 1 & \dots & 1 \\
1 & \omega_n & \omega_n^2 & \omega_n^ 3 & \dots & \omega_n^{n-1} \\
1 & \omega_n^2 & \omega_n^4 & \omega_n^6 & \dots & \omega_n^{2(n-1)} \\
1 & \omega_n^3 & \omega_n^6 & \omega_n^9 & \dots & \omega_n^{3(n-1)} \\
\vdots & \vdots & \vdots & \vdots & \ddots & \vdots \\
1 & \omega_n^{n-1} & \omega_n^{2(n-1)} & \omega_n^{3(n-1)} & \dots & \omega_n^{(n-1)(n-1)}
\end{bmatrix}
\begin{bmatrix} a_0 \\ a_1 \\ a_2 \\ a_3 \\ \vdots \\ a_{n-1} \end{bmatrix}
\]
观察矩阵中 $\omega_n$ 的幂次, 我们似乎发现 $\omega_n$ 的指数似乎组成了一张....乘法表???
因为我们将 $DFT$ 表示为 $\vec y=V_n \vec a$ , 则对于它的逆运算 $\vec a=DFT^{-1}(\vec y)$ , 我们可以选择直接求出 $V_n$ 的逆矩阵 $V^{-1}_n$ 并乘上 $\vec y$ 来计算.
然后就会有一个玄学的定理来告诉我们逆矩阵的形式, 算导上的证明如下:
定理: 对于 $j,k=0,1,...,n-1$, $V_n^{-1}$ 的 $(j,k)$ 处元素为 $\omega_n^{-kj}/n$.
证明: 我们证明 $V_n^{-1}V_n=I_n$, 其中 $I_n$ 为一个 $n\times n$ 的单位矩阵. 考虑 $V_n^{-1}V_n$ 中 $(j,j‘)$ 处的元素:
\[[V_n^{-1}V_n]_{jj‘}=\sum_{k=0}^{n-1}(\omega_n^{-kj}/n)(\omega_n^{kj‘})=\sum_{k=0}^{n-1}\omega_n^{k(j‘-j)}/n \]
当 $j‘=j$ 时, 此值为 $1$ , 否则根据求和引理, 此和为 $0$.注意, 只有 $-(n-1)\leq j‘-j \leq n-1$ , 从而使得 $j‘-j$ 不能被 $n$ 整除,才能应用求和引理.
$\blacksquare$
上面证明的思路是, 求出 $V_n^{-1}V_n$ 的各元素的值, 并和单位矩阵 $I_n$ 比较, 发现求出的值与单位矩阵相符, 原命题得证.
得到这个定理之后我们就可以推导出 $DFT^{-1}(\vec y)$ 了, 而我们发现它是长这样的:
\[a_j=\frac{1}{n} \sum_{k=0}^{n-1} y_k\omega_n^{-kj}\]
我们发现它和 $DFT$ 的定义极其相似! 唯一的区别仅仅在于前面的多出的 $\frac{1}{n}$ 和 $\omega_n$ 指数上多出来的负号. 也就是说, 我们可以通过稍微修改一下 $FFT$ 就可以用来计算 $DFT^{-1}$ 了!
最终我们可以得到像这样的代码, 通过第三个参数指定计算 $DFT$ 还是 $DFT^{-1}$ . 参数为 $1$ 则计算 $DFT$ , 参数为 $-1$ 则计算 $DFT^{-1}$. 除以 $n$ 的过程在函数执行后在函数外进行即可.
1 void FFT(Complex* a,int len,int opt){
2 if(len==1)
3 return;
4 Complex* a0=new Complex[len/2];
5 Complex* a1=new Complex[len/2];
6 for(int i=0;i<len;i+=2){
7 a0[i/2]=a[i];
8 a1[i/2]=a[i+1];
9 }
10 FFT(a0,len/2);
11 FFT(a1,len/2);
12 Complex wn(cos(2*Pi/len),opt*sin(2*Pi/len));
13 Complex w(1,0);
14 for(int i=0;i<(len/2);i++){
15 a[i]=a0[i]+w*a1[i];
16 a[i+len/2]=a0[i]-w*a1[i];
17 w=w*wn;
18 }
19 delete[] a0;
20 delete[] a1;
21 }
FFT 的速度优化与迭代实现方式
看到上面的实现, 根据OIer的常识, 我们显然可以意识到: 递归常数比较大. 而且上述实现中拆分多项式的时候还要重新分配一段空间, 造成了时间与空间上的多余开销. 那么, 我们是否可以通过模拟递归时对系数向量的重排与操作来取得更好的执行效率呢?
答案是肯定的.
首先我们可以注意到, 在上面的实现代码中计算了两次 $\omega_n^k y^{[1]}_k$, 我们可以将它只计算一次并将结果放在一个临时变量中. 然后我们的循环大概会变成这样:
1 for(int i=0;i<(len/2);i++){
2 int t=w*a1[i];
3 a[i]=a0[i]+t;
4 a[i+len/2]=a0[i]-t;
5 w=w*wn;
6 }
我们定义将旋转因子与 $y^{[1]}_k$ 相乘并存入临时变量 $t$ 中, 并从 $y^{[0]}_k$ 中增加与减去 $t$ 的操作为一个蝴蝶操作.
接着我们考虑如何将递归调用转化为迭代. 首先我们观察递归时的调用树:
我们可以发现如果我们自底向上观察这棵树, 如果将初始向量按照叶子的位置预先重排好的话, 完全可以自底向上一步一步合并结果. 具体的过程大概像这样:
首先成对取出元素, 对于每对元素进行 $1$ 次蝴蝶操作计算出它们的 $DFT$ 并用它们的 $DFT$ 替换原来的两个元素, 这样 $\vec a$ 中就会存储有 $n/2$ 个二元素 $DFT$ ;
接着继续成对取出这 $n/2$ 个 $DFT$ , 对于每对 $DFT$ 进行 $2$ 次蝴蝶操作计算出包含 $4$ 个元素的 $DFT$ , 并用计算出的 $DFT$ 替换掉原来的值.
重复上面的步骤, 直到计算出 $2$ 个长度为 $n/2$ 的 $DFT$ , 最后使用 $n/2$ 次蝴蝶操作即可计算出整个向量的 $DFT$.
实际实现时, 我们发现计算该轮要进行蝴蝶操作的两个元素的下标是比较方便的, 最外层循环维护当前已经计算的 $DFT$ 长度, 每次循环完成后翻倍; 次级循环维护当前所在的 $DFT$ 的开始位置的下标, 每次循环完成后加上 $DFT$ 长度值; 最内层循环枚举 $DFT$ 内的偏移量. 如果三个循环的循环变量分别为 $i,j,k$, 则 $j+k$ 指向当前要操作的第一个值, $j+k+i$ 指向下一个 $DFT$ 中需要计算的值.
然后关键的问题就在于如何快速重排得到递归计算时的叶子顺序.
我们以长度为 $8$ 的 $DFT$ 为例, 让我们打表找规律看看能不能发现什么:
| 原位置 | 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 |
| 原位置的二进制表示 | 000 | 001 | 010 | 011 | 100 | 101 | 110 | 111 |
| 重排后下标的二进制表示 | 000 | 100 | 010 | 110 | 001 | 101 | 011 | 111 |
| 重排结果 | 0 | 4 | 2 | 6 | 1 | 5 | 3 | 7 |
看起来似乎...是把二进制表示给翻过来了?
这个过程在算导上似乎叫做位逆序置换, 这个过程可以直接遍历一遍并交换二进制表示互为逆序的各个下标所对应的值, 也可以预处理出来保存在一个数组中来在多次定长 $FFT$ 中减少重复计算.
代码实现类似这样(我在这里将位逆序置换预处理掉存在 $rev$ 数组里了)
1 void FFT(Complex* a,int len,int opt){
2 for(int i=0;i<len;i++)
3 if(i<rev[i])
4 std::swap(a[i],a[rev[i]]);
5 for(int i=1;i<len;i<<=1){
6 Complex wn=Complex(cos(PI/i),opt*sin(PI/i));
7 int step=i<<1;
8 for(int j=0;j<len;j+=step){
9 Complex w=Complex(1,0);
10 for(int k=0;k<i;k++,w=w*wn){
11 Complex x=a[j+k];
12 Complex y=w*a[j+k+i];
13 a[j+k]=x+y;
14 a[j+k+i]=x-y;
15 }
16 }
17 }
18 }
预处理部分大概这样:
1 for(int i=0;i<bln;i++){
2 rev[i]=(rev[i>>1]>>1)|((i&1)<<(bct-1));
3 }
然后我们就成功地给 $FFT$ 加了个速OwO
总结
$FFT$ 作为快速计算卷积的算法, 在现在OI中较少考察高精度计算的情况下也很有用, 常常用于加速转移方程为卷积形式的动态规划问题.
而 $FFT$ 还有一个数论版本 $NTT$ , 利用的是模意义下的原根而不是单位复根. 我可能会在后续的博文中介绍关于 $NTT$ 的部分.
参考资料:
Introduction To Algorithms , 3rd Edition , Episode 30 , Polynomials and the FFT
如果看完觉得有点意思好像明白了什么的话请移步右下角点击"推荐"来滋磁博主的劳动OwO
如果你觉得文章中有错误也欢迎在评论区指正OwO