上次更新是在十月份？我到底摸了多久

卷积

对于两个序列 \(\{x\}, \{y\}​\)，定义它们的卷积 \(\{x \ast y\}​\) 为： \[ (x \ast y)_k = \sum_{j = 0}^k x_jy_{k- j} \]

单位根

定义复平面单位圆上的 \(n​\) 次单位根 \[ \omega_{n} = e^{\frac{2\pi i}{n}} = \cos\left(\frac{2\pi}{n}\right) + i\sin\left(\frac{2\pi}{n}\right) \] 易证其有如下的性质：

  1. \(\omega_{2n}^{2k} = \omega_n^k\)
  2. \(\omega_n^k = \omega_n^{k +n}\)
  3. \(\omega_n^k = - \omega_{n}^{k + n / 2}\)

DFT

对于一长度为 \(n\) 的序列 \(\{x\}\)，定义其离散傅里叶变换（DFT）\(\mathcal{F}\{x\} = \{X\}\) 为： \[ \mathcal{F}\{x\}_k = X_k = \sum_{j = 0}^{n - 1} x_i\omega_n^{jk} \] 则有其离散傅里叶逆变换（IDFT）为： \[ x_k = \frac{1}{n}\sum_{j = 0}^{n - 1} X_i\omega_n^{-jk} \] （标准定义里面 \(\omega_n^{jk}​\) 和 \(\omega_n^{-jk}​\) 应该互换，但是不影响）

暴力证明：考虑 \(\mathcal{F}\{x_n\}\) 的计算式： \[ \begin{pmatrix} 1 &1 &1 &\cdots &1 \\ 1 &\omega_n &\omega_n^2 &\cdots &\omega_n^{n - 1} \\ 1 &\omega_n^2 &\omega_n^4 &\cdots &\omega_n^{2(n-1)} \\ \vdots &\vdots &\vdots &\ddots &\vdots \\ 1 &\omega_n^{n - 1} &\omega_{n}^{2(n-1)} &\cdots &\omega_n^{(n-1)^2} \end{pmatrix} \begin{pmatrix} x_1 \\x_2 \\x_3 \\ \vdots \\x_{n - 1} \end{pmatrix} = \begin{pmatrix} X_1 \\X_2 \\X_3 \\ \vdots \\X_{n - 1} \end{pmatrix} \] 设左边矩阵为 \(V​\)，接下来从充分性上证明 \(V^{-1}_{jk} = \frac{\omega_n^{-jk}}{n}​\)： \[ \left(VV^{-1}\right)_{jk} = \sum_{l = 0}^{n - 1}\omega_n^{jl}\cdot\frac{\omega_n^{-lk}}{n} = \frac{1}{n}\sum_{l = 0}^{n - 1}\omega_n^{l(j - k)} \] 显然当 \(j = k​\) 时有 \(\left(VV^{-1}\right)_{jk} = 1​\)，而 \(j \neq k​\) 时，由等比数列求和公式： \[ \begin{aligned} \frac{1}{n}\sum_{l = 0}^{n - 1}\omega_n^{l(j - k)} &= \frac{1}{n} \cdot \frac{1 - \omega_n^{n(j - k)}}{1 - \omega_n^{j - k}} = 0 \end{aligned} \] 因此 \(VV^{-1} = I​\)，得证。

上面的 DFT 有如下意义：

  1. 将 \(\{x_n\}​\) 作为一组正交基 \(\{\omega_n^0, \omega_n^1, \cdots\}, \{\omega_n^0, \omega_n^{2}, \cdots\}​\) 的系数，计算和向量。（在此意义下逆变换就是将向量进行分解，因此逆变换的意义显然，无需上述证明）
  2. 将 \(\{x_n\}\) 作为以个多项式的系数，计算该多项式在 \(\omega_n^0, \omega_n^1,\cdots\) 点上的取值。

DFT 是线性的，即和的 DFT 等于 DFT 的和，且常数倍的 DFT 是 DFT 的常数倍。

同时 DFT 的另一个重要性质被称为卷积性，即 \[ \mathcal{F}\{x\}_k\mathcal{F}\{y\}_k = \mathcal{F}\{x \ast y\}_k \] 证明：把式子写出来： \[ \left(\sum_{j = 0}^{n - 1} x_j \omega_n^{jk}\right)\left(\sum_{j = 0}^{n - 1} y_j \omega_n^{jk}\right) \overset{?}{=} \sum_{j = 1}^{n - 1}\left(\sum_{l = 0}^jx_ly_{j - l}\right)\omega_n^{jk} \] 显然左边式子展开对应的任意一项 \(x_ay_b\omega_n^{(a + b)k}\) 都唯一地出现在右式 \(j = a+b\) 内层求和的展开式中，因此两个式子总是等价的。

DFT 将卷积 \(O(n^2)\) 的操作转化为 \(O(n)\) 的操作，因此可以用于加速卷积。

FFT

因为 DFT 对应点相乘这一操作本身是 \(O(n)​\) 的，因此我们只要找到快速的 DFT/IDFT 算法，就可以避免 \(O(n^2)​\) 地计算卷积。FFT 就是这样一种基于分治的快速 DFT 算法，时间复杂度为 \(O(n\log n)​\)。

首先假设我们需要计算一个序列 \(\{x\}\) 的 DFT（假设 \(n\) 是 \(2\) 的整数次幂），我们按照下标的奇偶性将这个序列分为 \(\{a\}, \{b\}\) 两部分，其中 \(a_j = x_{2j}\) 且 \(b_j = x_{2j + 1}\)。那我们有： \[ \begin{aligned} \mathcal{F}\{x\}_k &= \sum_{j = 0}^{n - 1} x_j\omega_n^{jk} \\ &= \sum_{j = 0}^{\frac{n}{2} - 1} \left[(x_{2j}\omega_n^{2jk} + x_{2j + 1}\omega_n^{(2j + 1)k}\right] \\ &= \sum_{j = 0}^{\frac{n}{2} - 1} x_{2j}\omega_n^{2jk} + \omega_n^k \sum_{j = 0}^{\frac{n}{2} - 1} x_{2j + 1}\omega_n^{2jk} \\ &= \sum_{j = 0}^{\frac{n}{2} - 1} a_j\omega_{n / 2}^{jk} + \omega_n^k \sum_{j = 0}^{\frac{n}{2} - 1} b_j\omega_{n / 2}^{jk} \\ &= \mathcal{F}\{a\}_k + \omega_{n}^k\mathcal{F}\{b\}_k \end{aligned} \] 看起来好像没有问题？我们最后一步推导其实是有一点小问题的，因为当 \(\mathcal{F}\{a\}\) 依据我们的定义只有 \(\frac{n}{2}\) 项，因此上面的推导只能让我们计算 \(0 \le k < \frac{n}{2}\) 的值。对于 \(k \ge \frac{n}{2}\) 的情况要如何处理呢？

首先，注意到： \[ \omega_n^{2jk} = \omega_n^{2jk - jn} = \omega{_n^{2j(k - n / 2)}} \] 因此我们稍微变一下： \[ \begin{aligned} \mathcal{F}\{x\}_k &= \sum_{j = 0}^{n - 1} x_j\omega_n^{jk} \\ &= \sum_{j = 0}^{\frac{n}{2} - 1} \left[x_{2j}\omega_n^{2j(k - n /2)} + x_{2j + 1}\omega_n^{2j(k - n /2) + k}\right] \\ &= \sum_{j = 0}^{\frac{n}{2} - 1} x_{2j}\omega_n^{2j(k - n /2)} + \omega_n^k \sum_{j = 0}^{\frac{n}{2} - 1} x_{2j + 1}\omega_n^{2j(k - n /2)} \\ &= \sum_{j = 0}^{\frac{n}{2} - 1} a_j\omega_{n / 2}^{j(k - n /2)} + \omega_n^k \sum_{j = 0}^{\frac{n}{2} - 1} b_j\omega_{n / 2}^{j(k - n /2)} \\ &= \mathcal{F}\{a\}_{k - n / 2} - \omega_{n}^{k - n / 2}\mathcal{F}\{b\}_{k - n / 2} \end{aligned} \] （其实本质上就是 DFT 的序列其实是可以周期拓展的）

这使得我们得以计算 \(\frac{n}{2} \le k < n\) 时 DFT 的值。因此，最为朴素的 FFT 实现就是基于递归不断将 \(\{x_n\}\) 进行分治，并通过上式对于分治结果进行合并。时间复杂度自然就是 \(O(n\log n)\)。

然而事实一定没有那么简单。递归写法的 FFT 常数很大，很大程度上是因为合并本身不是 in-place 的，需要额外申请空间。更快的 in-place FFT 实现，Cooley–Tukey 算法，是基于这样一个思路的：

假设每次分出来的 \(\{a\}\) 向左走，\(\{b\}\) 向右走，那递归版的 FFT 就会画出一棵很漂亮的二叉树，其中叶子节点对应的下标和序列本身的下标顺序是有很大差别的。我们可以预先快速计算出叶子节点的顺序，然后直接从叶子开始就地一层一层地合并上去。

那叶子节点的顺序遵循怎样的规律呢？其实我们依据奇偶分组本质上是在依据二进制的末位分组，在这样的思想下，我们脑补一下发现这就是叶子节点的下标把原来序列的下标的二进制翻转一下（低位变高位，高位变低位）（即 \(\{00, 01,10, 11\}\) 到叶子节点变成 \(\{00, 10, 01, 11\}\)）（可以结合画图理解）。这样的话我们维护一个从高位到低位的二进制加法器就行了。

代码如下：

#include <cstdio>
#include <complex>
#include <cmath>

using namespace std;
typedef complex<double> cplx;
const double PI = 3.1415926535897932384626;

void fft(int len, cplx *x, cplx *X, int dir = 1) {
    for (int i = 0, k = 0; i < len; i++) {
        X[i] = x[k];
        int y = len >> 1;
        while (y & k) k ^= y, y >>= 1; // 维护倒序加法器
        k |= y; 
    }
    for (int s = 2; s <= len; s <<= 1) {
        cplx w_(cos(2 * PI / s), dir * sin(2 * PI / s));
        for (int l = 0; l < len; l += s) {
            int mid = l + (s >> 1);
            cplx w = 1;
            for (int i = 0; i < s >> 1; i++) { // 逐层向上合并（一次蝴蝶操作）
                cplx t = w * X[mid + i];
                X[mid + i] = X[l + i] - t;
                X[l + i] = X[l + i] + t;
                w *= w_; 
            }
        }
    }
}