计算机系学生应当如何正确地卷

Cliggord Geertz 在 1963 年出版的工作 Agricultural Involution 让内卷（Involution）一词获得了它的现代意义。语义漂浮之中，“内卷”至今已当之无愧地成了现代汉语的一个热词。

又根据身边统计学调查，如何正确地对待 Involution 几乎已经成为了计算机系学生最关心的问题。而本文将讨论一个更强的问题：计算机系的学生应该如何正确对待 Convolution？

先解释笑话：Con- 作为前缀往往表达 thoroughly 之意，如 con-clude。那么，在这个意义上 con-volution 就可以理解为彻底的卷，因此对待 convolution 是比对待 involution 更强的问题。尽管实际上 convolution 作卷积讲时 con- 词缀应当取 with, be together 之意，类似于 con-nect。

---

一、开卷有益：卷积及其计算

卷积（Convolution）是一种函数运算，它将两个函数高度耦合并返回一个同时具有两个函数特征的新函数。作为计算机系的学生，我们尤其关注离散的情形，也就是可形式化地表达为

$f_n = \sum\limits_{i \in D} g_i \cdot h_{n - i}$

的所谓离散卷积。

这一形式最直观地呈现在多项式乘法中。考虑多项式 $A(x) = \sum\limits_{j = 0}^{n - 1} a_j x^j$ 与 $B(x) = \sum\limits_{j = 0}^{n - 1} b_j x^j$，其乘积应当满足 $C(x) = \sum\limits_{j = 0}^{2n - 2} c_j x^j = A(x) \cdot B(x)$。考虑 $A, B$ 中的每一对 $a_k x^k$ 与 $b_l x^l$ 相乘都会在 $C$ 中贡献 $a_kb_l x^{k + l}$ 一项。那么，如果我们想要求出 $c_j$ 的系数，也就只需枚举 $j$ 的成分有多少来自 $A$。于是：

$c_j = \sum\limits_{k = 0}^{j} a_k \cdot b_{j - k}$

容易发现多项式 $C$ 的系数形呈 $A$，$B$ 系数的卷积。

我们很自然地想要知道这个式子的计算方法。最朴素的办法就是对于每一个 $c_j$ 展开求和式，复杂度是不言自明的 $O(n^2)$。不过早在 1965 年就出现了一个复杂度为 $O(n \log n)$ 的算法：快速傅里叶变换（Fast Fourier transform, FFT）。

快速傅里叶变换来源于一个最简单的 insight：除了系数表示法以外，多项式还可以用点值表示。譬如对于多项式 $A(x) = 5x^3 - 2x + 1$，系数序列 $(5, -2, 1)$ 自然构成了一种表示；但同时我们知道，如果知道了 $n$ 多项式经过的 $n$ 个点，也可以唯一确定这一多项式——从而点列 $\{ (-1, -2), (0, 1), (1, 4) \}$ 也可以构成多项式 $A$ 的一个表示。这意味着什么呢？回忆我们对于多项式乘积的规定：

$C(x) = \sum\limits_{j = 0}^{2n - 2} c_j x^j = A(x) \cdot B(x)$

实际上可以拆分成三个等式，其意义分别为：

1. $C(x) = \sum\limits_{j = 0}^{2n - 2} c_j x^j$：乘积 $C$ 存在一种系数表示 $(c_0, \dots, c_{2n - 2})$。
2. $C(x) = A(x) \cdot B(x)$：对于 $A$ 的一个点值表示 $(x, A(x))$ 与 $B$ 的一个点值表示 $(x, B(x))$，它们的简单积 $(x, A(x)B(x))$ 一定构成 $C$ 的一个点值表示。
3. $\sum\limits_{j = 0}^{2n - 2} c_j x^j = A(x) \cdot B(x)$：系数表示与点值表示等价。

Eureka！这意味着：如果我们在多项式 $A$ 和 $B$ 上各采样出 $2n$ 个点，就可以在线性时间内得到 $C$ 的完整点值表示。不过想要真正求出 $c_j$ 的值，我们仍然面临一个主要的困难：如何实现多项式系数表示与点值表示之间的快速转化？这就引出了下述两个方法。他们分别被称作离散傅里叶变换（DFT）与离散傅里叶逆变换（IDFT），同时使用这两种方法以求解卷积的方式就成为快速傅里叶变换（FFT）。代码见附录。

傅里叶变换：从系数表示到点值表示

考虑一个 $n-1$ 次多项式。我们想要在多项式采样出 $n$ 个点，最自然的想法是考虑 $n$ 个不同的数作为自变量取值，这需要 $O(n)$ 时间。之后依次代入多项式中求解，每次多项式求值又需要 $O(n)$ 时间。此时总计复杂度是 $O(n^2)$ 的，实在不能接受。算法的瓶颈在于多项式求值的过程：而如果能适当的选取自变量，在求值时能复用一些结果就能实现更高效的算法。

这件事在实数域上似乎是不可能的——或许真的是这样。不过如果我们把论域扩充到复数域上，情况似乎就发生了改变。考虑取出方程 $\omega ^n = 1$ 的 $n$ 个原根（Plural Principal Root，即单位复根）作为自变量。我们知道 $e^{\pi i} = 1$ 那么记 $\omega_n^1 = e^{\pi i / n}$ ，$\omega_n^0, \omega_n^1, \dots, \omega_n^{n-1}$ 则能取到这 $n$ 个不同的原根。 不妨假设 $n = 2^k$（下文同。对于其余情况我们总可以通过在多项式高次补 0 的方式凑出 2 的整次幂），那么单位根就有两个重要的性质：

• $\omega_n^{2k} = \omega_{n/2}^k$
• $\omega_n^{k + n/2} = -\omega_n^k$

这里性质 1 的证明只需简单的指数变换；性质 2 的证明则可以从左右两式的平方入手。性质 1 提示我们可以从 $\omega_{n/2}$ 系列的 $\frac{n}{2}$ 个单位负根直接对应 $\frac{n}{2}$ 个 $\omega_n$ 系列的单位负根。与此同时性质 2 似乎又可以从已知的这一半单位根得到另一半单位根。如果在多项式求值时能利用这一性质，似乎就能指数减少 $n$ 的规模，得到一个高效的算法了。因此我们考虑将给定的多项式 $A(x) = \sum\limits_{j=0}^{n - 1} a_j x^j$ 变形为：

$A(x) = (a_0 + a_2 x^2 + \dots + a_{2k - 2} x^{2k - 2}) + (a_1 x + a_3 x^3 + \dots + a_{2k-1} x^{2k-1})$

我们记：

$\begin{aligned} F(x) &:= \sum\limits_{i = 0}^{k - 1} a_{2i} x^i = a_0 + a_2 x + a_4 x^2 + \dots + a_{2k - 2} x^{k - 1}\\ G(x) &:= \sum\limits_{i = 0}^{k - 1} a_{2i + 1} x^i = a_1 + a_3 x + a_5 x^2 + \dots + a_{2k - 1} x^{k - 1} \end{aligned}$

则有 $A(x) = F(x^2) + xG(x^2)$。将单位根代入，可以得到

$\begin{aligned} A(\omega_n^k) &= F(\omega_n^{2k}) + \omega_n^k \cdot G(\omega_n^{2k}) \\ &= F(\omega_{n/2}^k) + \omega_n^k \cdot G(\omega_{n/2}^k) \\[0.3cm] A(\omega_n^{k + n/2}) &= F(\omega_n^{2k}) - \omega_n^k \cdot G(\omega_n^{2k}) \\ &= F(\omega_{n/2}^k) - \omega_n^k \cdot G(\omega_{n/2}^k) \end{aligned}$

这意味着 $A(\omega_n)$ 与 $A(\omega_{n/2})$ 之间也存在简单直接的关系。于是在计算时只需递归地减少 $n$ 的规模就能以 $O(n \log n)$ 的复杂度计算出所有 $\omega_n$ 系列单位根对应的函数值了。这一过程可以由蝴蝶变换在相同渐进复杂度下获得更好的性能表现，这段内容我们留由附录叙述。

傅里叶逆变换：从点值表示到系数表示

要从点值表示求得系数表示其实很容易。考虑给定的点值表示 $\{ (\omega_n^0, y_0), (\omega_n^1, y_1), \dots, (\omega_n^{n - 1}, y_{n-1}) \}$ 对应了系数表示 $A(x) = \sum\limits_{j = 0}^{n - 1} a_j x^j$，那么就有

$\begin{bmatrix} y_{n - 1} \\ y_{n - 2} \\ \vdots \\ y_0 \end{bmatrix} = \begin{bmatrix} 1 & 1 & 1 & \dots & 1 \\ 1 & \omega_n^1 & \omega_n^2 & \dots & \omega_n^{n - 1} \\ 1 & \omega_n^2 & \omega_n^4 & \dots & \omega_n^{2(n - 1)} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & \omega_n^{n - 1} & \omega_n^{2(n - 1)} & \dots & \omega_n^{(n-1)^2} \\ \end{bmatrix} \begin{bmatrix} a_0 \\ a_1 \\ a_2 \\ \vdots \\ a_{n - 1} \end{bmatrix}$

记该式为 $\vec{y} = X \vec{a}$，这里 $X_{i, j} = \omega_n^{(i - 1)(j - 1)}$。

显然 $X$ 呈范德蒙德形 $V(1, \omega_n, \dots, \omega_n^{n-1})$，于是其行列式 $\det(X) = \prod\limits_{0 \le j < k < n-1} (\omega_n^j - \omega_n^k)$ 在 $\omega_n$ 互异时非零，于是 $X$ 可逆（这一结论与线性方程组 $\vec{y} = X \vec{a}$ 有唯一解等价，后者又内蕴于系数表示与点值表示的唯一对应性）。至此，通过 $\vec{a} = \left( X^{-1} \right) \vec{y}$ 就可以求出系数向量了。接下来考虑如何求解 $X^{-1}$。

这里我们考虑用拉格朗日插值法还原多项式并作系数对照。由插值法有：

$A(x) = \sum\limits_{i = 0}^{n-1} y_i \prod\limits_{i \ne j} \dfrac{x - x_j}{x_i - x_j}$

这里 $x_i$ 是第 $i$ 个复单位根。记 $f_i(x) = \prod\limits_{i \ne j} x - x_j$，即 $A(x) = \sum\limits_{i = 0}^{n-1} y_i \dfrac{f_i(x)}{f_i(x_i)}$。考虑到这里的 $x_j$ 即 $\omega_{n}^j$ 均为方程 $\omega ^n = 1$ 的解，于是 $x^n - 1 = \prod\limits_{i = 0}^{n -1} (x - x_i)$，从而有

$\begin{aligned} f_i(x) &= \prod\limits_{i \ne j} x - x_j = \dfrac{1}{x - x_i} \prod\limits_{i = 0}^{n -1} (x - x_i) \\ &= \dfrac{x^n - 1}{x - x_i} = x^{n-1} + x_i x^{n-2} + x_i^2 x^{n-3} + \dots + x_i^{n - 1} \end{aligned}$

那么可以得到：

$\dfrac{f_i(x)}{f_i(x_i)} = \dfrac{1}{n \cdot x_i^{n-1}} \left( x^{n-1} + x_i x^{n-2} + x_i^2 x^{n-3} + \dots + x_i^{n - 1} \right)$

对照该式与 $\left( X^{-1} \right) \vec{y}$ ，容易发现 $\left(X^{-1}\right)_{i, j}$ 显然就是上式的 $j-1$ 次项系数，于是有

$\left(X^{-1}\right)_{i, j} = \dfrac{x_i^{n - j}}{n \cdot x_i^{n-1}} = \dfrac{1}{n} x_i^{1-j} = \dfrac{1}{n} \omega_n^{-(i-1)(j-1)} = \dfrac{1}{n} X_{i, j}^{-1}$

回忆 $X_{i, j} = \omega_n^{(i-1)(j-1)}$，那么 $X^{-1}$ 也就是 $X$ 施简单线性变换后除以 $n$ 得到的。从而使用 $\vec{a} = \left( X^{-1} \right) \vec{y}$ 计算系数向量时，$\vec{a}$ 也可以视为 $Xy$ 的结果反置后除以 $n$ 得到的。这时 $Xy$ 是一个符合上一节形式的式子，可以用快速傅里叶变换求得。于是傅里叶逆变换最后转化为傅里叶变换，总的时间复杂度达到了 $O(n \log n)$。

二、手不释卷：更多形式的卷积

层层加码：递推的卷积

通常形式的卷积往往要计算形如

$c_n = \sum\limits_{i = 0}^{n} a_i \cdot b_{n - i}$

这样的式子。这里 $a, b, c$ 完全独立，可以使用朴素的 FFT 直接完成快速计算。但在实际应用中，我们可能会关心这样的一类式子：

$f_n = \sum\limits_{i = 0}^{n - 1} f_i \cdot g_{n - i}$

也就是数列 $f$ 的递推式中，第 $n$ 项是通过前序的各项 $f_i$ 卷积得来的。假设我们想要计算这一数列的第 $n$ 项，似乎就需要做 $n$ 次 FFT，复杂度为 $O(n^2 \log n)$。但实际上，我们有一种被称作分治 FFT 的做法可以做到 $O(n \log^2 n)$ 的复杂度。

考虑对于需要求的 $f_l, \dots, f_r$（记作区间 $[l, r]$） 分治。划分子问题为求解区间 $[l, m]$ 与 $(m, r]$，这里 $m = (l + r) / 2$。在求解 $(m, r]$ 时会依赖 $[l, m]$ 段求得的值。不过我们容易发现区间 $[l, m]$ 中所有值对于区间 $(m, r]$ 中任意一个元素 $f_x$ 的贡献都是相同的，为

$w_m = \sum\limits_{i = l}^{m} f_i \cdot g_{x - i}$

因此考虑先递归求解 $[l, m]$ 中的所有元素，之后再考虑求解 $(m, r]$ 时首先给各元素加上 $w_m$ 再递归求解即可。又求 $w_m$ 时使用 FFT 复杂度为 $O(n \log n)$，总的复杂度应满足

$T(n) = 2T(\dfrac{n}{2}) + O(n \log n)$

由主定理显然有 $T(n) \in O(n \log n)$。代码实现见附录。

美好的有限域：快速数论变换

快速数论变换(Number-theoretic transform, NTT) 是快速傅里叶基于数论在任意有限域上的实现。简单说来，快速傅里叶变换涉及大量浮点复数运算，计算机在处理时需要额外花费，而且可能会损失精度。如果我们规定需要计算的各个数值都在某个有限域中，则可以通过 NTT 得到更好的效率与精度。

这里我们考虑在模 $P$ 乘法群 $\langle \mathbb{Z}_P,~\times \rangle$ 中的情形，这里 $P$ 是一个形式特定的素数，我们将稍后讨论对 $P$ 的限制。形式地说，就是我们希望快速求解这样的式子：

$c_n \equiv \sum\limits_{i = 0}^{n} a_i \cdot b_{n - i} \mod P$

首先我们要找出复数域中的原根 $\omega_n = e^{\pi i / n}$ 在群 $\mathbb{Z}_P$ 中的对应。抽象代数给我们提供了一个名字与形式都相同的对应——有限域中的原根（生成元）。

原根 $g$ 是满足 $g^{|\mathbb{Z}_P|} \equiv 1 \mod P$ 且 $g^1, g^2, \dots, g^{|\mathbb{Z}_P|}$ 互异的元素。由于乘法群的阶满足 $|\mathbb{Z}_P| = \varphi(P)$，原根也可以等价定义为满足 $|g| = \varphi(P)$ 的元素。这里 $\varphi$ 是欧拉函数，$|g|$ 表示元素 $g$ 的阶，$|\mathbb{Z}_P|$ 表示群 $\mathbb{Z}_P$ 的阶。

假设我们已经求出了 $\mathbb{Z}_P$ 的某个原根 $g$，为了将其利用于求卷积的过程中，我们需要找出 $n$ 个互不相同的元素并能在形式上与 FFT 用到的几个单位复根性质一致。一种可能的构造是令 $g_n^i \equiv \left( g ^{\frac{P - 1}{n}} \right)^i \mod P$。此时：

• $g_n^{2k} \equiv \left( g ^{\frac{P - 1}{n}} \right)^{2k} \equiv \left( g ^{\frac{P - 1}{n/2}} \right)^k \equiv g_{n/2}^k \mod P$

• $g_n^{k + n/2} \equiv \left( g ^{\frac{P - 1}{n}} \right)^{k + n/2} \equiv g_n^k \cdot g ^{\frac{P - 1}{2}} \mod P$

由于 $P$ 是素数，考虑费马小定理有 $\left( g^{\frac{P - 1}{2}} \right)^2 \equiv g^{P - 1} \equiv 1 \mod P$。而我们知道 $1 = g^{P-1} \ne g^{\frac{P-1}{2}}$，于是就有 $g^{\frac{P-1}{2}} \equiv -1 \mod P$，从而 $g_n^{k + n/2} = - g_n^k$

从而我们可以认为原根 $g$ 生成的 $g_n^i$ 符合所需性质，可以直接在计算卷积时使用。我们就此得到了 NTT 算法。

但这里仍然有几个显著的问题。首先我们必须保证群 $\mathbb{Z}_P$ 存在原根，还须保证 $\dfrac{P - 1}{N}$ 为整数。

模 $P$ 乘法群存在原根的充要条件是 $P$ 形呈 $2, 4, p^\alpha, 2p^\alpha$（$p$ 是奇素数，$\alpha$ 是正整数）；由于在计算卷积时我们往往通过补零的方式将 $N$ 扩大为 2 的整次幂，使得 $(P-1) / N$ 为整数就需要 $P$ 形呈 $a2^b + 1$。这两个条件使得能使用本节介绍的 NTT 方法的乘法群非常稀少。常见的 $P$ 可能取 $998244353 = 119 \times 2^{23} + 1$ 或 $1004535809 = 479 \times 2^{21} + 1$。二者的原根都可取 $g = 3$。

不过 NTT 是可以扩展到任意模数有限域上的。一般说来会选取几个符合条件 $P$ 做多次 NTT 后用中国剩余定理（Chinese remainder theorem, CRT）合并。本节对此不作展开。

关于原根相关代数性质的证明见附录。

三、席卷八荒：无处不在的卷积

本章介绍一些卷积的应用场景。

整数乘法

一般认为整数乘法的复杂度是 $O(n^2)$ 的，$n$ 是数位个数。但整数可以视作特殊形式的多项式，因此可以用 $O(n \log n)$ 的时间通过 FFT 计算得到。

随机游走

题目描述：考虑一款 RPG 游戏，你可以给装备升级。若当前装备在第 $i$ 级，那么可以花费 $c_i$ 的代价进行一次升级。升级成功的概率为 $p_i$，装备等级加一；否则有 $(1 - p_i) w_j / \sum_{k=1}^i w_k$ 的概率降为 $i - j$ 级。现在你有一个 0 级的装备，试求出将其升级到 $n$ 级所需的期望花费。

这是一个经典的随机游走（Random walk, or the Drunkard’s walk）模型，它更具体地符合马尔科夫形式。本题的一个简单情况是有偏向的随机游走——即每次升级时增加等级的期望 $e_i$ 总是正数的情况。这时由期望的线性可加性，答案显然为

$E = \sum\limits_{i = 0}^{n - 1} \dfrac{c_i}{e_i} = \sum\limits_{i = 0}^{n - 1} \dfrac{c_i}{p_i - \dfrac{(1-p_i)}{\sum_{k=1}^i w_k}\sum\nolimits_{j = 1}^{i} j\cdot w_j }$

但实际情况往往没有这么简单。比如最经典的期望为 0 时的随机游走模型：考虑 $p_0 = 1, p_1 = 0.5$ 与 $c_0 = 1， c_1 = 5$ 且 $w_1 = 1$ 的情况。此时等级 0 向等级 1 的升级需要 1 的代价而必然成功，等级 1 向等级 2 升级时需要 5 的代价，有 $0.5$ 的概率成功；有 $0.5$ 的概率掉回 0 级。由定义，此时的答案应该是

$\begin{aligned} E &= (1 + 5) \times \dfrac{1}{2} + (1 + 5 + 1 + 5) \times \left( \dfrac{1}{2} \right)^2 + \dots \\ &= \sum_{i = 1}^{\infty} \dfrac{6i}{2^i} = 12 \end{aligned}$

计算中甚至需要使用无穷级数！这对于算法设计来说是相当困难的。

那么我们应该如何计算本题的期望呢？考虑前文提到的马尔科夫模型，用向量 $\vec{v_k} = \begin{bmatrix} p_n & p_{n-1} & \dots & p_0 \end{bmatrix}^T$ 表示经过 $k$ 次升级后装备等级的概率分布情况，$p_i$ 指示装备此时是 $i$ 级的概率。于是矩阵

$X = \begin{bmatrix} (1 - p_{n}) \cdot \frac{w_1}{\sum w_k} & p_{n-1} & 0 & 0 & \dots & 0 \\ (1 - p_{n}) \cdot \frac{w_2}{\sum w_k} & (1 - p_{n - 1}) \cdot \frac{w_1}{\sum w_k} & p_{n - 2} & 0 &\dots & 0 \\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots \\ (1 - p_{n}) \cdot \frac{w_n}{\sum w_k} &(1 - p_{n-1}) \cdot \frac{w_{n-1}}{\sum w_k} & (1 - p_{n-2}) \cdot \frac{w_{n - 2}}{\sum w_k} & \dots & \dots & 0 \end{bmatrix}$

就构成了 $\vec{v_{k + 1}} = X \vec{v_k}$ 的状态转移矩阵。于是就可以考虑通过不断迭代至某个精度上限的方式计算本题所需的期望。但我们当然有更好的办法。

考虑到矩阵做乘法时只可能涉及到期望的和运算（积运算是关于概率的），按这种状态转移的思路，考虑设出 $f_i$ 表示目前等级为 $i$ 时要升级到 $n$ 所需的期望花费，那么有转移：

$f_i = c_i + p_i \cdot f_{i + 1} + \dfrac{1 - p_i}{\sum_{k=1}^i w_k}\sum_{j = 0}^{i - 1} w_{i - j} \cdot f_{j}$

该式内涵于上述的 $\vec{v_{k + 1}} = X \vec{v_k}$ 中。进一步整理，有

$f_{i + 1} = \dfrac{1}{p_i} \left( f_i - c_i - \dfrac{1 - p_i}{\sum_{k=1}^i w_k}\sum_{j = 0}^{i - 1} w_{i - j} \cdot f_{j} \right)$

构成了一个含有卷积运算的递推式。且依次排列前 $n$ 个递推式可以构成一个下三角形的方程组。而观察方程组形式，我们可以断言有： $f_i = k_i f_0 + b_i$ 成立。从而

$\begin{aligned} & k_0 = 1 \\ & b_0 = 0 \\ & k_{i + 1} = \dfrac{1}{p_i} \left( k_i - \dfrac{1 - p_i}{\sum_{k=1}^i w_k}\sum_{j = 0}^{i - 1} w_{i - j} \cdot f_{j} \right ) \\ & b_{i + 1} = \dfrac{1}{p_i} \left( b_i - c_i - \dfrac{1 - p_i}{\sum_{k=1}^i w_k}\sum_{j = 0}^{i - 1} w_{i - j} \cdot f_{j} \right ) \\ \end{aligned}$

而这组式子这是可以用前文中提到的分治 FFT快速求解的。

汉明距离

题目描述：给定 01 比特串 $S$, $T$。问至少要修改 $T$ 中的多少个比特位才能让 $T$ 成为 $S$ 的子串。

问题即求 $S$ 的任一长度为 $|T|$ 的子串与 $T$ 之间的汉明距离。于是答案可以形式化为

$\min\limits_{0 \le i \le |S| - |T|} \sum_{j = 1}^{|T|} S_{i + j} \oplus T_j$

这里 $\oplus$ 表示异或运算。观察到表达式很像卷积形式，我们考虑令 $T^{-1}$ 表示串 $T$ 翻转后的结果，令

$c_i = \sum_{j = 1}^{|T|} S_{i + j} \oplus {\left(T^{-1}\right)}_{|T| - 1 -j}$

答案即求 $\{c_i\}$ 的最小值。假如这里的异或能换成乘法，就可以用 FFT 快速求解了。所以我们如此将它变为乘法：

$c_i = \left (\sum_{i = j + k} S_j(1 - T_k) \right ) + \left ( \sum_{i = j + k} (1 - S_j)T_k \right)$

从而就可以用 FFT 在 $O(n \log n)$ 时间内求出答案。

烷烃计数

摸了。有空补上。

图像处理与卷积神经网络

不是很会。会了补上。

附录

FFT的蝴蝶变换

蝴蝶变换又称位逆序置换（bit-reversal permutation），本质上是将递归方法用循环改写。下举出一例说明该过程：

譬如考虑递归求解 $\{\omega_n^0, \omega_n^1, \omega_n^2, \omega_n^3, \omega_n^4, \omega_n^5, \omega_n^6, \omega_n^7\}$ 的过程，

第一次递归得子问题： $\{\omega_n^0, \omega_n^2, \omega_n^4, \omega_n^6\},~\{\omega_n^1, \omega_n^3, \omega_n^5, \omega_n^7\}$

第二次递归得子问题： $\{\omega_n^0, \omega_n^4\},~\{\omega_n^2, \omega_n^6\},~\{\omega_n^1, \omega_n^5\},~ \{\omega_n^3, \omega_n^7\}$

第三次递归得子问题：$\{\omega_n^0\},~\{\omega_n^4\},~\{\omega_n^2\},\{\omega_n^6\},~\{\omega_n^1\},~\{\omega_n^5\},~\{\omega_n^3\},\{\omega_n^7\}$

我们发现根 $\omega_n^k$ 在最后一次递归时出现的位置序号就是 $k$ 的二进制表示反置后的位置。如 $\omega_n^1$ 的 001 反置后为 100，他也就出现在最终序列中下标为 4 的位置上。根据这个性质就可以直接得到最终递归的序列，并通过非递归实现求解原问题。

可以察觉到蝴蝶变换的本质，实际来源于在 $i$ 次划分子问题时，算法会依据二进制表达的低 $i$ 位进行。在这个理解上证明直接使用归纳法即可。

FFT的代码实现

下给出一个基于 C++ STL 的 FFT 实现，这里我们使用了蝴蝶变换。

#include <bits/stdc++.h>
using namespace std;

typedef complex<double> Complex;
const double Pi = acos(-1.0);
const int N = 1e7 + 10;

Complex F[N], G[N];
int rev[N];

void FFT(Complex *a, int n, int op) {
    for (int i = 0; i != n; ++i) {
        rev[i] = rev[i >> 1] >> 1;
        if (i & 1) rev[i] |= (n >> 1);
    }

    for (int i = 0; i != n; ++i) {
        if (i < rev[i]) {
            swap(a[rev[i]], a[i]);
        }
    }

    for (int mid = 2; mid <= n; mid <<= 1) {
        Complex wn(cos(2 * Pi / mid), op * sin(2 * Pi / mid));
        for (int i = 0; i < n; i += mid) {
            Complex w(1, 0);
            for (int j = i; j < i + mid / 2; j++, w *= wn) {
                Complex x = a[j], y = w * a[j + mid / 2];
                a[j] = x + y;
                a[j + mid / 2] = x - y;
            }
        }
    }
}

int main() {
    double x;
    int n, m, limit = 1;

    scanf("%d%d", &n, &m);
    for (int i = 0; i <= n; i++) {
        scanf("%lf", &x);
        F[i].real(x);
    }
    for (int i = 0; i <= m; i++) {
        scanf("%lf", &x);
        G[i].real(x);
    }
    while (limit < n + m + 1) {
        limit <<= 1;
    }
    
    FFT(F, limit, 1);
    FFT(G, limit, 1);
    for (int i = 0; i != limit; ++i) {
        F[i] = F[i] * G[i];
    }
    FFT(F, limit, -1);
    
    for (int i = 0; i <= n + m; i++) {
        printf("%d ", (int)std::round(F[i].real()));
    }
    puts("");
    return 0;
}

分治FFT实现

这里仅给出关键的递归求解函数。

void solve(int l, int r) {
    if (l == r) return;

    int mid = (l + r) >> 1;
    solve(l, mid);

    for (int i = l; i <= mid; i++) {
        A[i - l] = f[i];
    }
    for (int i = 1; i <= r - l; i++) {
        B[i - 1] = g[i];
    }

    int lim = 1, cnt = 0;
    while (lim <= r - l + mid - l + 1) {
        lim <<= 1;
        cnt++;
    }
    for (int i = 0; i < lim; i++) {
        R[i] = (R[i >> 1] >> 1) | ((i & 1) << cnt - 1);
    }

    FFT(A, 1, lim);
    FFT(B, 1, lim);
    for (int i = 0; i < lim; i++) {
        A[i] = A[i] * B[i];
    }
    FFT(A, -1, lim);

    for (int i = mid - l; i <= r - l - 1; i++) {
        f[i + l + 1] = f[i + l + 1] + A[i];
    }

    for (int i = 0; i < lim; i++) {
        A[i] = B[i] = 0;
    }
    
    solve(mid + 1, r);
}

原根

摸了。有空补上。

NTT的代码实现

inline int ADD(int x, int y) { return x + y >= P ? x + y - P : x + y; }

inline int SUB(int x, int y) { return x - y < 0 ? x - y + P : x - y; }

void get_root() {
    int lim = 1, cnt = 0;
    while (lim <= 2 * n) lim <<= 1, cnt++;
    for (int mid = 1; mid < lim; mid <<= 1) {
        int W = qpow(3, (P - 1) / (mid << 1));
        wn[mid] = 1;
        for (int j = 1; j < mid; j++) {
            wn[mid + j] = 1ll * wn[mid + j - 1] * W % P;
        }
    }
}

int qpow(int a, int x) {
    int ans = 1;
    while (x) {
        if (x & 1) ans = 1ll * ans * a % P;
        a = 1ll * a * a % P;
        x >>= 1;
    }
    return ans;
}

void NTT(int *a, int type, int lim) {
    for (int i = 0; i < lim; i++) {
        if (i < R[i]) {
            swap(a[i], a[R[i]]);
        }
    }

    for (int mid = 1; mid < lim; mid <<= 1)
        for (int i = 0; i < lim; i += (mid << 1))
            for (int j = 0; j < mid; j++) {
                int x = a[i + j], y = 1ll * a[i + mid + j] * wn[mid + j] % P;
                a[i + j] = ADD(x, y);
                a[i + mid + j] = SUB(x, y);
            }
    if (type == -1) {
        int INV = qpow(lim, P - 2);
        for (int i = 0; i < lim; i++) {
            a[i] = 1ll * a[i] * INV % P;
        }
        reverse(a + 1, a + lim);
    }
}