[学习笔记] 多项式与快速傅里叶变换(FFT)基础

时间 2019-11-13

标签学习笔记多项式快速傅里叶变换 fft 基础繁體版

原文原文链接

==== €€£ WARNING ====html

这篇博文内容相对偏少, 已经在后续博文中扩充.
算法

你们能够看个人最新博文 [学习笔记&教程] 信号, 集合, 多项式, 以及各类卷积性变换 (FFT,NTT,FWT,FMT)数组

==== ====函数

引入

可能有很多OIer都知道FFT这个神奇的算法, 经过一系列玄学的变化就能够在 $O(nlog(n))$ 的总时间复杂度内计算出两个向量的卷积(或者多项式乘法/高精度乘法), 而代码量却很是小. 博主一年半前曾经因COGS的一道叫作"神秘的常数 $\pi$"的题目而去学习过FFT, 可是基本就是照着板子打打完并不知道本身在写些什么鬼畜的东西OwO学习

不过...博主这几天忽然照着算法导论本身看了一遍发现本身彷佛忽然意识到了什么OwO而后就打了一道板子题还1A了OwO再加上午考试差点AK以及日更频率即将不保因而就有了这篇博文233优化

博主在写这篇文字的过程当中也发现了很多本身以前忽略的FFT细节, 感受对FFT彷佛有了更深入的理解, 但愿想学习FFT的读者能认真看完这篇文章并仔细体味FFT的优雅性质.spa

本文可能并不会介绍关于使用FFT进行信号处理的相关信息, 主要介绍在OI中的应用即快速卷积.code

Rush了很久终于写出来了QAQorm

那么如今, 就让咱们变玄学为科学了解 $FFT$ 背后的工做原理与它的优雅性质htm

前置知识: 多项式的表示

多项式的系数表达

相信你们都知道对于一个次数界为 $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\]

由于次数界折半了因此代入的是 $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]} \\ &= 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