FFT 学习笔记
前置知识
- 复数运算
- 单位根
简介
快速傅里叶变换（Fast Fourier Transform，FFT）是一种用于在 \(O(n \log n)\) 时间内计算离散傅里叶变换（Discrete Fourier Transform，DFT）的算法。在算法竞赛中，常用于快速计算多项式卷积。
朴素算法
例如有两个多项式，\(F(x),G(x)\)。
设 \(F(x)\) 次数为 \(n\)，\(G(x)\) 次数为 \(m\)。
显然，要求他们的卷积 \(H(x) = F(x) \ast G(x)\)，可以用 \(O(nm)\)（即 \(n^2\) 级别）的复杂度解决：
\[ H(k) = \sum_{i=0}^{k} F(i)G(k-i) \]
但，这实在太慢了，而目前看来都没有什么很好的做法，我们不妨改变一下多项式的表示方法。
点值表示法
平时我们用的 \(F(x) = a_0 + a_1x + a_2x^2 + \cdots + a_{n-1}x^{n-1}\) 被称作 系数表示法，它可以用一个向量来描述：
\[ \mathbf{a} = \begin{bmatrix} a_0 \\ a_1 \\ \vdots \\ a_{n-1} \end{bmatrix} \]
现在，我们可以取 \(n\) 个互异的 \(x\) 值（记为 \(x_0, \dots, x_{n-1}\)），并计算出对应的 \(F(x)\) 值，那么，就构成了 \(n\) 组 \((x_i, F(x_i))\) 的点对。这就是多项式的 点值表示法，它有许多特点：
- 唯一性：对于次数 \(< n\) 的多项式，任取 \(n\) 组互异点对，多项式唯一确定。
- 方便计算：
- 加法：\((F+G)(x_i) = F(x_i) + G(x_i)\)，复杂度 \(O(n)\)。
- 乘法：\((F\cdot G)(x_i) = F(x_i)G(x_i)\)，复杂度 \(O(n)\)。
瓶颈分析
虽然点值表示法下的乘法很快，但是：
- 系数转点值（求值）：如果我们随意取 \(n\) 个互异的 \(x\) 值，代入计算一个点的值需要 \(O(n)\)（秦九韶算法），计算 \(n\) 个点则总共需要 \(O(n^2)\) 的复杂度。
- 点值转系数（插值）：从点值推回系数，一般需要 \(O(n^2)\) 的复杂度（拉格朗日插值）。
如果只用朴素的方法，总复杂度依然是 \(O(n^2)\)，并没有变快。 我们需要一种特殊的取点方式，使得“求值”和“插值”都能在 \(O(n \log n)\) 内完成，这就是 FFT 的作用。
离散傅里叶变换 (DFT)
专业地说：
离散傅里叶变换（Discrete Fourier transform，DFT）是傅里叶变换在时域和频域上都呈离散的形式。 —— OI Wiki
在算法竞赛中，我们可以简单理解为：它是一种选取特殊的点（单位根），将 系数表示法 转到 点值表示法 的变换。
它是一种 线性 变换，可以用矩阵乘法来表示。
设系数向量为 \(\mathbf{a}\)，变换后的点值向量为 \(\mathbf{y}\)，则 \(\mathbf{y} = V \mathbf{a}\)：
\[ \begin{bmatrix}y_0 \\y_1 \\y_2\\\vdots \\y_{n-1}\end{bmatrix}=\begin{bmatrix}1 & 1 & 1 & \cdots & 1 \\1 & \omega & \omega^2 & \cdots & \omega^{n-1} \\1 & \omega^2 & \omega^4 & \cdots & \omega^{2(n-1)} \\\vdots & \vdots & \vdots & \ddots & \vdots \\1 & \omega^{n-1} & \omega^{2(n-1)} & \cdots & \omega^{(n-1)(n-1)}\end{bmatrix}\begin{bmatrix}a_0 \\a_1 \\a_2 \\\vdots \\a_{n-1}\end{bmatrix} \]
其中，\(\omega = \omega_n = e^{i\frac{2\pi}{n}}\)（\(n\) 次单位根）。
快速傅里叶变换 (FFT)
快速傅里叶变换（Fast Fourier Transform，FFT），并不是一种新的变换，而是利用单位根的性质（周期性、对称性），高效计算上述 DFT 矩阵乘法的算法。本质思想是分治。
分治法
FFT 的核心思想是分治。我们就 \(n\) 为 \(2\) 的整数幂的情况进行讨论（若不足则补 0）。
我们的目标是计算 \(x\) 取所有单位根 \(\omega_n^k\) (\(k=0, \dots, n-1\)) 时 \(F(x)\) 的值。
对于多项式： \[ F(x) = \sum_{i=0}^{n-1} a_i x^i \]
我们将系数按照下标的 奇偶性 分类：
- 偶数项系数：\(a_0, a_2, a_4, \dots\)
- 奇数项系数：\(a_1, a_3, a_5, \dots\)
将 \(F(x)\) 重写为： \[ F(x) = (a_0 + a_2 x^2 + \dots + a_{n-2}x^{n-2}) + (a_1x + a_3x^3 + \dots + a_{n-1}x^{n-1}) \]
在奇数项部分提公因式 \(x\)： \[ F(x) = (a_0 + a_2 x^2 + \dots) + x(a_1 + a_3 x^2 + \dots) \]
如果我们定义两个新的多项式 \(G(x)\) 和 \(H(x)\)，大小仅为原来的一半（项数为 \(n/2\)）： \[ \begin{align*} G(x) &= a_0 + a_2 x + a_4 x^2 + \dots = \sum_{i=0}^{n/2 - 1} a_{2i} x^i \\ H(x) &= a_1 + a_3 x + a_5 x^2 + \dots = \sum_{i=0}^{n/2 - 1} a_{2i+1} x^i \end{align*} \]
那么原多项式就可以表示为： \[ F(x) = G(x^2) + x \cdot H(x^2) \]
这是 FFT 最关键的一步。 我们可以利用单位根的性质，将求 \(F\) 的 \(n\) 个点值的问题，转化为求 \(G\) 和 \(H\) 的 \(n/2\) 个点值的问题。
对于 \(k < \frac{n}{2}\)，我们要计算 \(F(\omega_n^k)\) 和 \(F(\omega_n^{k + n/2})\)：
代入 \(x = \omega_n^k\)： 利用性质 \(\omega_n^{2k} = \omega_{n/2}^k\)（折半引理）： \[ \begin{align*} F(\omega_n^k) &= G((\omega_n^k)^2) + \omega_n^k \cdot H((\omega_n^k)^2) \\ &= G(\omega_n^{2k}) + \omega_n^k \cdot H(\omega_n^{2k}) \\ &= G(\omega_{n/2}^k) + \omega_n^k \cdot H(\omega_{n/2}^k) \end{align*} \]
代入 \(x = \omega_n^{k + n/2}\)： 利用性质 \(\omega_n^{2(k+n/2)} = \omega_n^{2k+n} = \omega_n^{2k}\omega_n^n = \omega_{n/2}^k\) 和 \(\omega_n^{n/2} = -1\)： \[ \begin{align*} F(\omega_n^{k + n/2}) &= G(\omega_n^{2k + n}) + \omega_n^{k + n/2} \cdot H(\omega_n^{2k+n}) \\ &= G(\omega_{n/2}^k) - \omega_n^k \cdot H(\omega_{n/2}^k) \end{align*} \]
也就是说，只要知道 \(G(\omega_{n/2}^k)\) 和 \(H(\omega_{n/2}^k)\)，就能推导出 \(F(\omega_n^k)\) 和 \(F(\omega_n^{k + n/2})\)，即：
\[ \begin{cases} F(\omega_n^k) &= G(\omega_{n/2}^k) + \omega_n^k H(\omega_{n/2}^k) \\ F(\omega_n^{k + n/2}) &= G(\omega_{n/2}^k) - \omega_n^k H(\omega_{n/2}^k) \end{cases} \]
注意到分治左右两边长度应当相等，所以长度只能是 \(2\) 的幂，不足要补齐。
另外单位根是复数，由于 STL 提供的 Complex 常数巨大，需要手写。
倍增法实现
假设你轻松的写出了分治法的 FFT，你会发现你 TLE 了。究其原因是递归的常数太大了，我们考虑不用递归。
我们考虑递归到最底层时，系数数组的顺序是怎样的，以 \(8\) 次多项式为例：
- 第一层：\((x_0,x_1,x_2,x_3,x_4,x_5,x_6,x_7)\)
- 第二层：\((x_0,x_2,x_4,x_6),(x_1,x_3,x_5,x_7)\)
- 第三层：\((x_0,x_4),(x_2,x_6),(x_1,x_5),(x_3,x_7)\)
- 第四层：\((x_0),(x_4),(x_2),(x_6),(x_1),(x_5),(x_3),(x_7)\)
观察其二进制表示：
0 (000)->0 (000)1 (001)->4 (100)2 (010)->2 (010)3 (011)->6 (110)可以发现，最终序列的下标恰好是原序列下标的 二进制位翻转 形式，我们称这个变换为 位逆序置换（bit-reversal permutation）。
简单证明
每一层奇偶拆分，等价于根据当前下标的 最低有效位 (LSB) 进行分组。LSB 为 0 的分到左边，为 1 的分到右边。经过 \(\log n\) 次拆分，原下标的最高位最终决定了它在新序列中的最低位，反之亦然。
补充
这里补充如何快速知道其二进制位的翻转形式，设 \(r_x\) 为 \(x\) 二进制翻转之后的数。
显然，\(r_0 = 0\)，我们从小到大求 \(r\)。
我们只要把 \(x\) 右移一位，翻转，再右移一位，就能得到 \(x\) 二进制下除了个位的翻转结果。
- 个位如果是 \(0\)，反转后最高位为 \(0\)。
- 个位如果是 \(1\)，反转后最高位为 \(1\)，因此加上 \(2^{len -1}\)。
即：
\[ R(x) =\lfloor \frac{R(\lfloor\frac{x}{2}\rfloor)}{2} \rfloor + (x \bmod 2) \cdot 2^{len -1} \]
容易发现我们可以 \(O(n)\) 获得最后的系数顺序，我们就可以从下往上依次合并，我们还需解决合并问题，回到这个式子：
\[ \begin{cases} F(\omega_n^k) &= G(\omega_{n/2}^k) + \omega_n^k H(\omega_{n/2}^k) \\ F(\omega_n^{k + n/2}) &= G(\omega_{n/2}^k) - \omega_n^k H(\omega_{n/2}^k) \end{cases} \]
观察使用位逆序置换后：
- \(G(\omega_{n/2}^k)\) 存在数组下标为 \(k\) 的位置，\(H(\omega_{n/2}^k)\) 存在数组下标为 \(k + n/2\) 的位置。
那么我们可以：
- 将 \(F(\omega_n^k)\) 存在数组下标为 \(k\) 的位置，\(F(\omega_n^{k + n/2})\) 存在数组下标为 \(k + n/2\) 的位置。
这样我们就可以直接覆盖，而不用额外数组，我们称这种方法为 蝶形运算。
快速傅里叶逆变换 (IFFT)
快速傅里叶逆变换（Inverse Fast Fourier Transform，IFFT）是 FFT 的逆过程，其目标是将 点值表示法 还原为 系数表示法。
从线性代数角度看，DFT 是一个线性变换，可以表示为 \(\mathbf{y} = V \mathbf{a}\)，其中 \(\mathbf{a}\) 是系数向量，\(\mathbf{y}\) 是点值向量，\(V\) 是 DFT 的变换矩阵：
\[ V_{jk} = \omega_n^{jk} \quad \text{其中 } \omega_n = e^{i\frac{2\pi}{n}} \]
那么，IFFT 就是找到一个逆矩阵 \(V^{-1}\)，使得 \(\mathbf{a} = V^{-1} \mathbf{y}\)。
幸运的是，这个 DFT 矩阵的逆矩阵具有一个非常优美的形式。设 \(V^{-1}\) 是一个矩阵，其元素为：
\[ (V^{-1})_{jk} = \frac{1}{n} \omega_n^{-jk} \quad \text{其中 } \omega_n^{-1} = e^{-i\frac{2\pi}{n}} \]
也就是说，逆矩阵 \(V^{-1}\) 只需要将原矩阵 \(V\) 中的单位根 \(\omega_n\) 替换为其共轭复数 \(\omega_n^{-1}\)，然后将整个矩阵除以 \(n\) 即可。
证明
基于单位根的一个关键性质，有时被称为 求和引理：
\[ \sum_{k=0}^{n-1} (\omega_n^j)^k = \begin{cases} n & \text{if } j \equiv 0 \pmod n \\ 0 & \text{otherwise} \end{cases} \]
发现相乘为单位矩阵。
实现
由于 FFT 和 IFFT 的转移矩阵只有正负不一样，所以我们在编写 FFT 函数时，可以多加形参 type=1/-1 来区分 FFT 和 IFFT，给单位根取不同的符号。
记得 IFFT 最后要除以 \(n\)。
1 |
|