【笔记/模板】多项式Ⅰ：快速傅里叶变换和快速数论变换

前言

这一节用于讲解拉格朗日插值法（Lagrange Polynomial）和快速傅里叶变换（Fast Fourier Transform），但是含有前置知识，因此有大量学过的知识可以直接跳过，存在大量证明会放出相应的链接。

多项式

基本概念

我们将形如 $\textstyle{\sum a_nx^n}$ 的有限项相加式子成为多项式，记作 $f(x) = \textstyle{\sum_{i = 0}^n a_i x^i}$。

如果是形如 $\textstyle{\sum_{i = 0}^\infty a_i x^i}$ 的无限项相加式，则称为形式幂级数，后面应该会讲到。

• $a_i$ 被称为第 $i$ 次项的系数，记作 $[x ^ i]f(x) := a_i$，对于 $\forall a_i(i \gt n)$，均为 $0$。

• 多项式的系数非零的最高次项的次数为该多项式的度，也为该多项式的次数，记作 $\deg f$。

• 使得 $f(x) = 0$ 的所有 $x$​ 成为多项式的根。

• 如果 $\forall a_i(0 \lt i \le n), a_i \in \mathbb{R}$，我们称 $f$ 为实系数多项式，否则称为复系数多项式。

以上都是初中知识范畴里的，可以直接跳过。

代数基本定理

任何非零一元 $n$ 次复系数多项式恰好有 $n$ 个可重合复数根。

证明有兴趣详见代数基本定理的证明，证明需要的知识不在我们学习之内。但是这个定理对于后面的多项式插值以及像 FFT 等变换很重要的推动作用。

多项式的两种表示方法

第一种也是最常见的一种方法：系数表示法。这种方式直接表示出了多项式 $f(x)$ 的每一项系数，因此任意一点的点值我们都可以在线性时间内求出。同样，这种方式也方便进行计算机的存储，让次数和系数一一对应。

第二种同样十分重要：点值表示法。我们选取多个不同的 $x_i$，将 $x_i$ 带入进 $f(x)$ 得出对应的点值 $(x_0, y_0), (x_1, y_1), (x_2, y_2), \dots$。问题在于：应该选取多少个不同的 $x_i$ 才可以唯一的确定一个多项式 $f$ 呢？

至于如何确定，这就相当于解一个 $n$ 元线性方程组，一共多少组才能确定唯一解。

先给出结论：$n$ 个点值唯一确定的多项式的度数为 $n - 1$。

一个感性的理解是：$n - 1$ 次多项式通过系数表示法需要 $n$ 个系数表述，因此 $n$ 个点值表述出的多项式也应为 $n - 1$ 次的。

详细的证明采用范德蒙德矩阵，具体地：对于任意 $0 \le j \lt n, \textstyle{\sum_{i = 0}^{n - 1} a_ix_j^i = y_j}$，则该系数矩阵为：

$$\begin{bmatrix} 1 & x_0 & x_0^2 & \cdots & x_0^{n-1} \\ 1 & x_1 & x_1^2 & \cdots & x_1^{n-1} \\ 1 & x_2 & x_2^2 & \cdots & x_2^{n-1} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & x_{n-1} & x_{n-1}^2 & \cdots & x_{n-1}^{n-1} \end{bmatrix}$$

其中 $x_i$ 互不相同，根据范德蒙德矩阵的性质，矩阵行列式不为 $0$， 该方程组存在唯一解。同样地，这也证明了 $n$ 个点值不可能唯一确定度数至少为 $n$​ 的多项式。

这两种表示法显然可以互相转换，我们把从系数表示法转换为点值表示法的过程称为求值（Evaluation）。

相反地，我们把逆过程成为插值（Interpolation）。

基本运算

系数表示法下，设两个多项式 $f(x) = \textstyle{\sum_{i = 0} ^ n a_i x^i}$ 和 $g(x) = \textstyle{\sum_{i = 0} ^ m b_i x^i}$，它们的度数分别为 $n$ 和 $m$。

• 多项式之间可以进行加法，定义 $h := f + g = \textstyle{\sum_{i = 0} ^ n a_i x^i} + \textstyle{\sum_{i = 0} ^ m b_i x^i}$，因此

  $$h = \sum_{i = 0}^{\max(n, m)} (a_i + b_i) x^ i$$

  可以得到 $\deg h = \deg (f + g) = \max(\deg f, \deg g)$。

• 多项式之间可以进行乘法，定义 $h := f \cdot g = (\textstyle{\sum_{i = 0} ^ n a_i x^i}) \times (\textstyle{\sum_{i = 0} ^ m b_i x^i})$，因此

  $$h = \sum_{i = 0}^{n + m} \left ( \sum_{j = 0}^i a_j b_{i - j}\right ) x^i$$

  可以得到 $\deg h = deg(f \cdot g) = \deg f + \deg g$。

因此在系数表示法下，我们可以在 $\mathcal{O}(\max(n, m))$ 的情况下进行加法，在 $\mathcal{O}(mn)$​ 的情况下进行乘法。

点值表示法下，设 $f(x):(x_0, y_0), (x_1, y_1),\dots,(x_n, y_n)$，和 $g(x):(x_0, y_0'), (x_1, y_1'),\dots,(x_m, y_m')$，它们的度数分别为 $n,m(n \le m)$。

• 它们之间进行加法，$h(x) := f(x) + g(x) = (f + g)(x)$，因此点值对应相加，表示为：

  $$h(x):(x_0, y_0 + y_0'), (x_1, y_1 + y_1'),\dots,(x_n, y_n + y_n'),\dots,(x_m, y_m + y_m')$$

• 它们之间进行乘法，$h(x) := f(x)g(x) = (fg)(x)$，因此点值对应相乘，表示为：

  $$h(x):(x_0, y_0 y_0'), (x_1, y_1 y_1'),\dots,(x_n, y_n y_n'),\dots,(x_m, y_m y_m')$$

  值得注意的是，因为 $\deg h = \deg f + \deg g$，因此进行乘法时要选出一共 $n + m + 1$ 个点值相乘才可以唯一确定 $h(x)$。

因此在点值表示法下，我们可以在线性的时间复杂度下完成加法和乘法运算。

复数与单位根

FFT 算法的最巧妙之一在于它很好的结合了复数的大量性质，本小节直接忽略复数的基础知识。

复平面

任意一个复数 $z \in \mathbb{Z}, z = a + bi$，我们可以将其放在平面直角坐标系上观察，此时 $z$ 的坐标为 $(a, b)$。这样很多特征就很明显了：

• 定义复数 $z$ 的模 $|z| := \sqrt{a ^ 2 + b ^ 2}$，在复平面上的几何意义表示该点到原点的距离，以下记为 $r$。
• 定义复数 $z$ 的幅角 $\arg z := \theta$，满足 $\tan \theta = \frac{b}{a}$，即原点到 $z$ 的连线与 $x$ 轴的夹角。

根据高中所学，$a = r \cos \theta, b = r\sin \theta$，因此可以表示复数的三角形式：$z = r(\cos \theta + i \sin \theta)$。

直接引入棣莫弗定理：复数相乘，模长相乘，幅角相加。

$$z1 = r1(\cos \theta_1 + i \sin \theta_1), z2 = r2(\cos \theta_2 + i \sin \theta_2) \\ z1 \times z2 = r1 \times r2 \left ( \cos(\theta_1 + \theta_2) + i \sin(\theta_1 + \theta_2) \right )$$

这个定理较为常见，但对于单位根的性质十分重要。

单位根

不妨钦定 $|z| = 1$，此时所有满足条件的 $z$ 都在单位圆上，考虑对其进行幂运算，根据棣莫弗定理不难发现，$z ^ k = \cos ^ k \theta + i \sin ^ k \theta$，仍然在单位圆上。并且由于幅角相加，$z ^ k$ 可以看作在复平面上从 $(1, 0)$ 开始以 $z$ 的幅角 $\theta$ 旋转了 $k$ 次。

如果将这个单位圆 $n$ 等分，那么这 $n$ 个等分点的 $n$ 次方均为 $(1, 0)$，看似神奇，其实证明是十分显然的。我们假设这 $n$ 个等分点为 $P_0, P_1, \dots,P_{n - 1}$，其中 $P_1 ^ n$ 刚好绕单位元一圈，而 $P_i(1 \lt i \lt n)$ 均饶了单位元整数圈。

写出它们的复平面坐标：$P_k = (\cos \frac{2k \pi}{n}, i \sin \frac{2k \pi}{n})$，直接将 $\omega_n := P_1$，发现 $P_k = \omega_n^k = P_1 ^ k$。我们将这些所有这些 $P_k$ 称之为关于 $n$ 的单位根。

这些点的性质如下：

• 存在循环，且循环周期为 $n$：$\omega_n^k = \omega_n^{kt}(t \gt 0)$，$\omega_n^k = \omega_n^{k + nt}(0 \le k \lt n, t \gt 0)$。
• $w_n^k = -w_n^{k + \frac{n}{2}}(0 \le k \lt \frac{n}{2}, n \text { is even})$​。
• $\forall w_n^k(0 \le k \lt n)$，这 $n$ 个复数共同组成了 $n$ 次方程 $x^n = 1$ +的 $n$ 个根，因为 $(w_n^k)^n = 1(0 \le k \lt n)$。

着重注意第二、三性质，对于后续问题很重要。

本原单位根

本原单位根指的是集合 $\{ \omega_n ^k \ | \ 0 \le k \lt n,n \perp k \ \}$ 中的元素。特殊之处在于通过其中任意一个都可以得出其余所有的单位根，因此无论如何，$w_n^1$ 必然是本原单位根。

拉格朗日插值法

前面说到过，系数表示法和点值表示法可以进行转换，其中取值可以在 $\mathcal{O}(\deg f)$ 内完成，然而插值就十分复杂了，因为我们要通过高斯消元解出一个 $n$ 元线性方程组，时间复杂度 $\mathcal{O}(n^3)$。

而拉格朗日插值法提供了一种特殊的方法，能让我们在 $\mathcal{O}(n ^ 2)$ 甚至线性的复杂度下通过点值表示法求出多项式的系数表示法。

算法过程

由于点值表示法的点值可加性，我们可以把一个有 $n + 1$ 组点值的多项式 $f$ 看作 $n + 1$ 个多项式相加， 可以表示为：

$$f_i(x) = \begin{cases} y_i, & \text{if } x = x_i \\ 0, & \text{otherwise} \end{cases}, \quad 0 \leq i \leq n$$

将这些多项式相加即为得到的原来多项式，容易证明这是对的，并且生成的多项式是唯一的。

那么如何构造出这些 $f_i(x)$ 呢，我们借用 CRT 的思想，将所有 $\textstyle{\prod_{j \neq i} x - x_j}$ 插入多项式，同时保证 $x = x_i$ 时的值为 $y_i$，我们插入 $\textstyle{\prod_{j \neq i} \frac{1}{x_i - x_j}}$ 加入多项式当中。因此，我们成功构造出来了，即

$$f_i(x) = y_i \prod_{j \neq i} \frac{x - x_j}{x_i - x_j}$$

因此合成后可以得到目标函数：

$$f = \sum_{i = 0}^n f_i = \sum_{i = 0}^n \left ( y_i \prod_{j \neq i} \frac{x - x_j}{x_i - x_j} \right )$$

至于如何求出 $f$ 的各个系数，我们预处理出辅助函数 $F(x) := \textstyle{\prod_{i = 0}^n (x - x_i)}$ 存好每项的系数，对于每个 $f_i$ 花费 $\mathcal{O}(n)$ 代价暴力除掉 $x - x_i$ 得到 $f_i$ 的各个项系数，最后相加即可，因此时间复杂度为 $\mathcal{O}(n ^ 2)$。如果只用来求特定点上的值，将 $x$ 带入即可，时间复杂度不变。

特殊情况

如果给出的点对是连续的，例如 $\forall i(0 \le i \lt n), x_i = i$，满足很多性质，我们可以直接做到 $\mathcal{O}(n)$。

$$f = \sum_{i = 0} ^ n f_i = \sum_{i = 0} ^ n y_i \prod_{j \neq i} \frac{x - j}{i - j}$$

分子是经典的 $\textstyle{\prod_{j = 0}^n x - j}$ 除去一个 $x - i$，使用前缀后缀积维护即可。

分母因为连续同样可以预处理，分正负阶乘相乘即可，最终公式为：

$$\begin{aligned} f = \sum_{i = 0} ^ n f_i &= \sum_{i = 0} ^ n y_i \frac{pre_{i - 1} \times suf_{i + 1}}{\prod_{j = 0} ^ {i - 1}(i - j) \times \prod_{j = i + 1} ^ {n} (i - j)} \\ &= \sum_{i = 0} ^ n \frac{pre_{i - 1} \times suf_{i + 1}}{i! \times (-1)^{n - i + 1} \times (i - n + 1)!} \end{aligned}$$

由于 $\mathcal{O}(1)$ 求出每一个 $f_i$，总共的复杂度为 $\mathcal{O}(n)$。

代码

P4781 【模板】拉格朗日插值 - 洛谷

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20

signed main()
{
    read(n, K);
    for (int i = 1; i <= n; i ++) read(x[i], y[i]);

    int ans = 0;
    for (int i = 1; i <= n; i ++)
    {
        int multi = 1, div = 1;
        for (int j = 1; j <= n; j ++)
        {
            if (i == j) continue;
            multi = multi * (K - x[j] + mod) % mod;
            div = div * (x[i] - x[j] + mod) % mod;
        }
        ans = (ans + multi * y[i] % mod * ksm(div, mod - 2) % mod) % mod;
    }
    cout << ans << '\n';
    return 0;
}

快速傅里叶变换（FFT）

如果想从纯数学知识或者通信方面了解 FFT，请尝试先了解离散型傅里叶变换（DFT，Discrete Fourier Transform）的概念，FFT 是 DFT 的改进和转化（等我会了再说）。

而如果像我一样从计算机算法方面，只需要知道 FFT 可以用来优化多项式乘法就行了。

给定多项式 $f$ 和 $g$，求出多项式 $h = f \times g$ 的各个项的系数。朴素做法很简单，时间复杂度为 $\mathcal{O}(n ^ 2)$ 的。而快速傅里叶变换算法让这个时间复杂度降为了 $\mathcal{O}(n \log n)$。

算法方向

前文提到过，系数表示法下的多项式乘法只能做到平方级别，然而点值表示法可以线性完成。这启发我们，如果我们对于多项式 $f,g$，花费一定代价分别求出 $\deg f + \deg g$ 个点值对，将这些点对对应相乘，此时乘法的时间复杂度就降到了线性了。但是这种考虑不够充分，因为我们要预先对两个多项式分别求 $\deg f + \deg g + 1$ 次求值，然而每一次求值都花费 $\mathcal{O}(n)$ 的时间，因此时间复杂度仍然为 $\mathcal{O}(n ^ 2)$。

不过很有启发性，我们能否找到一种方法，使得预先的求值可以做到 $\mathcal{O}(n\operatorname{polylog}(n))$ 呢？

第一次折半

有的。

考虑一个极其普通的二次多项式 $f(x) = x ^ 2$，我们自然需要三个点来唯一表示它，例如 $x = \{ 0, 1, 2\}, f(x) = \{ 0, 2, 4\}$，这样我们需要计算一共 $O(n ^ 2)$ 次。但如果我们取值为 $x = \{-1, 0, 1\}$ 呢？不难发现 $f(-1) = f(-(1)) = f(1) = 1$，因为该多项式为偶函数并且 $-1,1$ 互为相反数，它们的点值也相同。就这样，我们通过取特殊的 $x$ 从而省略掉了一次的点值计算。

不仅是偶函数，奇函数也如此，例如 $f(x) = x ^ 3 + x$ 为奇函数，其中 $f(-1) = f(-(1)) = -f(1)$，因此也可以简化。

既然发现了这个性质，我们推广到普遍的多项式 $f(x) = \textstyle{\sum_{i = 0}^n a_ix^i}$，为了让其保持偶函数或者奇函数的特征，我们拆分这个多项式，将偶数项放在一起，奇数项发在一起，即：

$$f(x) = \sum_{i = 0}^n a_ix^i \\ f_e(x) := a_0x^0 + a_2x^2 + a_4x^4 + \dots,f_o(x) := a_1x^1 + a_3x^3 + a_5x^5 + \dots \\ \therefore f(x) = f_e(x) + f_o(x)$$

根据高中所学，显然 $f_e(x)$ 和 $f_o(x)$ 分别满足偶函数和奇函数的性质。我们先随意选取共 $\lceil \frac{n}{2}\rceil$ 对互为相反数的 $x$ 作为所求的点值，然后分别带入 $f_e$ 和 $f_o$ 当中，求出相应值后带回原多项式中：

$$\begin{cases} f(x) = f_e(x) + f_o(x) \\ f(-x) = f_e(x) - f_o(x) \end{cases}$$

讨论此时的时间复杂度，$f_e$ 和 $f_o$ 项数减半，同时选取的点值减半，因此此时的复杂度为 $\mathcal{O}(2 \times \lceil \frac{n}{2}\rceil \lceil \frac{n}{2}\rceil) = \mathcal{O}(\frac{n ^ 2}{2}) = \mathcal{O}(n ^ 2)$。总共的时间复杂度没有变化，但是实际的常数减半了。我们猜想如果这个过程如果可以递归下去，那么只需要递归共 $\log n$ 层，时间复杂度可以优化至 $T(n) = 2 \times T(\frac{n}{2}) + \mathcal{O}(n) = \mathcal{O}(n \log n)$。

继续递归

在第一层递归时，我们选取了互为相反数作为点值，从而将求解多项式 $f(x)$ 点值一分为二，划分为了求 $f_e(x)$ 和 $f_o(x)$ 的两个子问题，我们接着考虑。

对于偶函数 $f_e(x) = a_0 + a_2x^2 + a_4x^4 + \dots$，因为都为偶数次项，我们定义 $f_e'(x) := f_e(x ^ 2) = a_0 + a_2x + a_4x^2 + \dots$，这个形式和最初的 $f(x)$ 很像，但是此时我们加入一对相反数 $(p, -q)$，它们在 $f_e'(x) = f_e(x ^ 2)$ 由于平方的影响导致 $p ^ 2 = q ^ 2$，这对相反数不成立。

换句话说，目前我们要找到一对不同数 $(p, -p)$ 满足 $p ^ 2 = -(-p) ^ 2$，显然实数域不存在，因此考虑复数内的 $i$，显然 $i^2 = - (-i)^2$。如图所示 ，设 $f(x) = x ^ 3 + x ^ 2 - x - 1$， 我们目前所在第二层，选取两对相反数 $(x1, -x1), (x2, -x2)$ 同时满足 $(x_1^2, x_2^2)$ 互为相反数，因此选择了 $(1, i, -1, -i)$。在第一层时，我们选择 $(1, i, -1, -i)$，在第二层的 $f_e, f_o$ 计算中我们选择 $(1^2, i^2)$，在最后一层中因为只有一项，并且这一项的 $x = 1 ^ 4 = i ^ 4$​。

继续拓展至 $n$ 次多项式，注意到递归是一个满二叉树的形式，因此我们补满空余的取值，令 $n \gets 2 ^ {\lceil \log n \rceil}$​，复杂度不变。

如何选择 $(x_1, -x_1, x_2, -x_2, x_3, -x_3, x_4, -x_4)$ 满足互不相同且有如上性质呢？答案是单位根，因为它们均要满足 $x_k^n = 1$。

不仅如此，$w_n^k = -w_n^{k + \frac{n}{2}}(0 \le k \lt \frac{n}{2}, n \text { is even})$ 这一性质很好的让我们确定了 $(x1, x2, x3, x4)$ 和它们对应的相反数，它们这两个点在单位圆上关于原点对称。

当 $w_n^k$ 递归至下一层时，根据棣莫弗定理，$w_n^k \to w_n^{k^2} = (w_n^k)^2$，而 $(w_n^k) ^ 2 = w_{\frac{n}{2}}^k$，也就是说明单位根的幅角扩大一倍，而决定的顺序不会发生改变，并且仍然存在它的相反数，如图所示。

像这样，最开始我们给多项式设定的 $n$ 个值就应当是单位圆上关于 $n$ 的 $n$ 个单位根 $\omega_0, \omega_1, \dots, \omega_{n - 1}$，递归到最深一层后都会变成 $1$，这就是为什么最深一层取值只需要 $\mathcal{O}(n)$ 次，而一共 $\log n$ 层每一层花费 $\mathcal{O}(n)$ 进行归并处理。

大致流程

根据单位根的特殊性质证明了可行性，我们再梳理一遍 FFT 的大致流程。

1. 为方便递归，同时保持对称性质，令 $L \gets 2 ^ {\lceil \log n + m + 1 \rceil}$，多项式系数中默认补 $0$。

2. 设当前递归的多项式为 $border$ 项，如果 $border = 1$ 则结束递归，返回 $a_0$ 即可。

3. 否则处理出所有的单位根 $w_{border}^k (0 \le k \lt \frac{border}{2})$，对于当前的多项式 $f‘(x) = f_e'(x) + f_o'(x)$，我们有：

   $$f’(w_{border}^k) = f_e'(w_{border}^{2k}) + w_{border}^k \times f_o'(w_{border}^{2k})$$

   计算出所有 $0 \lt k \lt \frac{border}{2}$ 的 $f'(\omega_{border}^{k})$，由于 $w_n^k = -w_n^{k + \frac{n}{2}}$，我们可以通过这个得出所有的 $f'(\omega_{border}^{k + \frac{border}{2}})$，具体地：

   $$\begin{aligned} f'(w_{border}^{k + \frac{border}{2}}) &= f'_e(w_{border}^{2k + border}) + w_{border}^{k + \frac{border}{2}} \times f'_o(w_{border}^{2k + border}) \\ &= f_e'(w_{border}^{2k} \times w_{border}^{border}) - w_{border}^k \times f_o'(w_{border}^{2k} \times w_{border}^{border}) \\ &= f_e'(w_{border}^{2k}) - w_{border}^k \times f_o'(w_{border}^{2k}) \end{aligned}$$

   通过这两个等式：

   $$\begin{cases} f'(w_{border}^k) = f_e'(w_{border}^{2k}) + w_{border}^k \times f_o'(w_{border}^{2k}) \\ f'(w_{border}^{k + \frac{border}{2}}) = f_e'(w_{border}^{2k}) - w_{border}^k \times f_o'(w_{border}^{2k}) \end{cases} (0 \le k \lt \frac{border}{2})$$

因此猜想成立，时间复杂度：$T(n) = 2 \times T(\frac{n}{2}) + \mathcal{O}(n) = \mathcal{O}(n \log n)$，实现代码采取递归方式：

1
2
3
4
5
6
7
8
9
10
11
12
13

void FFT(int border, vector<Complex> &A, int type)    // A 中为原多项式，存出点值后递归返回
{
    if (border == 1) return;
    vector<Complex> F(border / 2), G(border / 2);    // 新建两个多项式存储 f_e 和 f_o
    for (int i = 0; i < border; i += 2) F[i >> 1] = A[i], G[i >> 1] = A[i + 1];    // 分奇偶存储
    
    FFT(border / 2, F, type), FFT(border / 2, G, type);    // 递归调用
    Complex base = Complex(cos(2.0 * PI / border), sin(2.0 * PI / border) * type), cur = Complex(1, 0);
        // base = w_1, cur 表示当前是哪一个单位根

    for (int i = 0; i < border / 2; i ++, cur = cur * base)    // cur = w_k
        A[i] = F[i] + G[i] * cur, A[i + (border / 2)] = F[i] - G[i] * cur;    // 最终公式可知
}

迭代优化

虽然我们的递归算法很好地实现了 FFT，将时间复杂度控制在了 $\mathcal{O}(n \log n)$ 以内，但是每一层新建的多项式以及递归的大常数很大程度地影响了代码的运行效率。注意到代码的结构和线段树类似，那么我们能否找到一种方法，使得算法可以通过循环从下往上迭代而非递归？

有的，我们称之为蝴蝶变换。

由于从下往上递归，原先的递归版本中递归前的这么一句话：

1

for (int i = 0; i < border; i += 2) F[i >> 1] = A[i], G[i >> 1] = A[i + 1];    // 分奇偶存储

要求我们在迭代前预处理出来，换句话说，对于原多项式系数序列（这里钦定 $n = 16$）：

$$(a_0, a_1, a_2, a_3, a_4, a_5, a_6, a_7, a_8, a_9, a_{10}, a_{11}, a_{12}, a_{13}, a_{14}, a_{15}) \Rightarrow \\ (a_0, a_2, a_4, a_6, a_8, a_{10}, a_{12}, a_{14}), (a_1, a_3, a_5, a_7, a_9, a_{11}, a_{13}, a_{15}) \Rightarrow \\ \dots \Rightarrow \\ (a_0, a_8), (a_2, a_{10}), (a_4, a_{12}), (a_6, a_{14}), (a_1, a_9), (a_3, a_{11}), (a_5, a_{13}), (a_7, a_{15})$$

重排成以上序列。

$$\begin{array}{|c|c|c|c|} \hline i \text{ (十进制)} & i \text{ (二进制)} & r[i] \text{ (十进制)} & r[i] \text{ (二进制)} \\ \hline 1 & 0001 & 8 & 1000 \\ 2 & 0010 & 4 & 0100 \\ 3 & 0011 & 12 & 1100 \\ 4 & 0100 & 2 & 0010 \\ 5 & 0101 & 10 & 1010 \\ 6 & 0110 & 6 & 0110 \\ 7 & 0111 & 14 & 1110 \\ 8 & 1000 & 1 & 0001 \\ 9 & 1001 & 9 & 1001 \\ 10 & 1010 & 5 & 0101 \\ 11 & 1011 & 13 & 1101 \\ 12 & 1100 & 3 & 0011 \\ 13 & 1101 & 11 & 1011 \\ 14 & 1110 & 7 & 0111 \\ 15 & 1111 & 15 & 1111 \\ \hline \end{array}$$

列出表格后凭借惊人的注意力发现 $i$ 和 $r_i$ 的二进制表示恰好反转，因此可以 $\mathcal{O}(n)$ 预处理出 $r_i$ 和序列顺序。

1
2

for (int i = 0; i < T; i ++) R[i] = (R[i >> 1] >> 1) | (i & 1 ? T >> 1 : 0);
for (int i = 0; i < T; i ++) if (i < R[i]) swap(A[i], A[R[i]]);

这就是蝴蝶变换，同样迭代实现的 FFT 算法也可以完成了：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21

void FFT(vector<Complex> &A, int type)
{
    vector<int> R(T);
    for (int i = 0; i < T; i ++) R[i] = (R[i >> 1] >> 1) | (i & 1 ? T >> 1 : 0);
    for (int i = 0; i < T; i ++) if (i < R[i]) swap(A[i], A[R[i]]);    // 蝴蝶变换

    for (int cur = 1; cur < T; cur <<= 1)    // cur = border / 2
    {
        Complex base(cos(PI / cur), sin(PI / cur) * type);    // 当前的单位根

        static vector<Complex> comp(T + 1);
        comp[0] = Complex(1, 0);
        for (int i = 1; i < cur; i ++) comp[i] = comp[i - 1] * base;    // 预处理，防止重复计算
        for (int i = 0; i < T; i += cur << 1)    // 当前处理区间的左下标
            for (int j = 0; j < cur; j ++)    // 当前区间的第 k 个数
            {
                Complex tempx = A[i | j], tempy = A[i | j | cur] * comp[j];
                A[i | j] = tempx + tempy, A[i | j | cur] = tempx - tempy;
            }    // 这里的或运算因为没有冲突，所以也是加法，位运算常数更小
    }
}

以上过程就称为 对长度为 $n$ 的多项式 $f$ 做快速傅里叶变换。

快速傅里叶逆变换（IFFT）

结束了么，并没有。我们在 $\mathcal{O}(n \log n)$ 内求出了多项式的点值，花费 $\mathcal{O}(n)$ 进行点值乘法，那么我们还要将结果的多项式插值回去。但是先前介绍的拉格朗日插值法需要 $\mathcal{O}(n)$ 的时间复杂度，而且单位根导致我们的点值不是连续的。因此接下来引入快速傅里叶逆变换（Inverse Fast Fourier Transform），也就是 FFT 的逆操作。

讲到 IFFT 了就不得不引入范德蒙德矩阵了，作为和多项式紧密联系的矩阵，我们通过矩阵乘法可以求出多项式的点值：

$$P(x) = p_0 + p_1x^1 +p_2x^2+ \dots + p_{n - 1}x^{n - 1}\\ \begin{bmatrix} P(x_0) \\ P(x_1) \\ P(x_2) \\ \vdots \\ P(x_{n-1}) \end{bmatrix} = \begin{bmatrix} 1 & x_0 & x_0^2 & \cdots & x_0^{n-1} \\ 1 & x_1 & x_1^2 & \cdots & x_1^{n-1} \\ 1 & x_2 & x_2^2 & \cdots & x_2^{n-1} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & x_{n-1} & x_{n-1}^2 & \cdots & x_{n-1}^{n-1} \end{bmatrix} \begin{bmatrix} p_0 \\ p_1 \\ p_2 \\ \vdots \\ p_{n-1} \end{bmatrix}$$

下面令：

$$\mathcal{F} = \begin{bmatrix} 1 & x_0 & x_0^2 & \cdots & x_0^{n-1} \\ 1 & x_1 & x_1^2 & \cdots & x_1^{n-1} \\ 1 & x_2 & x_2^2 & \cdots & x_2^{n-1} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & x_{n-1} & x_{n-1}^2 & \cdots & x_{n-1}^{n-1} \end{bmatrix}$$

既然范德蒙德矩阵可以用来求多项式点值，那么它的逆 $\mathcal{F}^{-1}$ 应该可以求出多项式的插值了。

考虑求出 $\mathcal{F}^{-1}$ 的每个元素，先使用拉格朗日插值法暴力求出每个元素，可以得到：

$$f = \sum_{i = 0}^n f_i = \sum_{i = 0}^n \left ( y_i \prod_{j \neq i} \frac{x - x_j}{x_i - x_j} \right ) \\ L_i := \prod_{j \neq i} \frac{x - x_j}{x_i - x_j}$$

其中 $L_i$ 为基函数，满足：

$$L_i(x) = \begin{cases} 1 & \text{if } x = x_j \\ 0 & \text{otherwise} \end{cases}$$

不难发现，这 $n$ 个基函数的值一次排列形成的矩阵刚好是单位矩阵 $\mathcal{I}$，因此可以得到：

$$\mathcal{F} \times \mathcal{F}^{-1} = \mathcal{I}$$

解出来得到：

$$(\mathcal{F}^{-1})_{ij} = \frac{ c_j }{ (x_i - x_j)c_i } \quad \text{for } i \neq j, \quad (\mathcal{F}^{-1})_{ii} = \frac{1}{c_i} \\ c_i = \prod_{\substack{0 \leq k \leq n-1 \\ k \neq i}} (x_i - x_k) = (x_i - x_0)(x_i - x_1) \cdots (x_i - x_{i-1})(x_i - x_{i+1}) \cdots (x_i - x_{n-1})$$

$c_i$ 其实就是基函数 $L_i(x)$ 的分母。

看似没有什么特性，但是我们把单位根代入进去，这就是 FFT 的矩阵：

$$\omega := \omega_{N}, N := border\\ \mathcal{F} = \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}$$

算出它的逆 $\mathcal{F}^{-1}$，惊奇的发现是：

$$\mathcal{F}^{-1} = \frac{1}{N} \begin{bmatrix} 1 & 1 & 1 & \cdots & 1 \\ 1 & \omega^{-1} & \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}$$

也就是说，我们把 FFT 算法中的所有 $\omega$ 替换为 $\omega^{-1}$，由于这里单位根的模长为 $1$，因此 $\omega^{-1} = \bar{\omega}$，最后结果再乘上 $\frac{1}{N}$ 就可以得到 IFFT 的结果了。

如果你看不懂上述矩阵，我们可以逆推出来：

对于多项式 $f(x) = \textstyle{\sum_{i = 0}^{n - 1}} a_ix^i$，我们对其做一次 FFT 得到多项式 $g(x) = \textstyle{\sum_{i = 0}^{n - 1} b_ix^i}$，其中 $b_i = \textstyle{\sum_{k = 0}^{n - 1}a_k (\omega_n^{k})^i}$，紧接着，我们将所有 FFT 中的 $\omega$ 换成 $\omega^{-1}$ 对多项式 $g$ 再进行一次 FFT 得到多项式 $h(x)$：

$$h(x) = \sum_{i = 0}^{n - 1} c_ix^i \\ \begin{aligned} c_i &= \sum_{k = 0}^{n - 1} b_k (\omega_n^{-k})^i \\ &= \sum_{k = 0}^{n - 1} \left ( \sum_{j = 0}^{n - 1} a_j(\omega_n^j)^k \right )(\omega_n^{-k})^i \\ &= \sum_{j = 0}^{n - 1} \left ( a_j \sum_{k = 0}^{n - 1} (\omega_{n}^{j - i})^k \right ) \end{aligned}$$

同理求出可以，通过拉格朗日插值法的基函数也可以证明：当 $i = j$ 时这个式子才会对 $c_i$ 产生贡献，因此

$$c_i = n a_i \Rightarrow a_i = \frac{c_i}{n}$$

我们在求完之后对整个 FFT 得到的序列除以 $n$ 就可以得到 IFFT 转换的结果了。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21

void FFT(vector<Complex> &A, int type)    // type 取 1 或者 -1 实现了这种功能
{
    vector<int> R(T);
    for (int i = 0; i < T; i ++) R[i] = (R[i >> 1] >> 1) | (i & 1 ? T >> 1 : 0);
    for (int i = 0; i < T; i ++) if (i < R[i]) swap(A[i], A[R[i]]);

    for (int cur = 1; cur < T; cur <<= 1)
    {
        Complex base(cos(PI / cur), sin(PI / cur) * type);

        static vector<Complex> comp(T + 1);
        comp[0] = Complex(1, 0);
        for (int i = 1; i < cur; i ++) comp[i] = comp[i - 1] * base;
        for (int i = 0; i < T; i += cur << 1)
            for (int j = 0; j < cur; j ++)
            {
                Complex tempx = A[i | j], tempy = A[i | j | cur] * comp[j];
                A[i | j] = tempx + tempy, A[i | j | cur] = tempx - tempy;
            }
    }
}

快速多项式乘法

这自然是 FFT 和 IFFT 在信息竞赛最常见的应用，大致流程和代码已经给出。

P3803 【模板】多项式乘法（FFT） - 洛谷 代码：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41

int T = 1, n, m;

void FFT(vector<Complex> &A, int type)
{
    vector<int> R(T);
    for (int i = 0; i < T; i ++) R[i] = (R[i >> 1] >> 1) | (i & 1 ? T >> 1 : 0);
    for (int i = 0; i < T; i ++) if (i < R[i]) swap(A[i], A[R[i]]);

    for (int cur = 1; cur < T; cur <<= 1)
    {
        Complex base(cos(PI / cur), sin(PI / cur) * type);

        static vector<Complex> comp(T + 1);
        comp[0] = Complex(1, 0);
        for (int i = 1; i < cur; i ++) comp[i] = comp[i - 1] * base;
        for (int i = 0; i < T; i += cur << 1)
            for (int j = 0; j < cur; j ++)
            {
                Complex tempx = A[i | j], tempy = A[i | j | cur] * comp[j];
                A[i | j] = tempx + tempy, A[i | j | cur] = tempx - tempy;
            }
    }
}

signed main()
{
    read(n, m);
    while (T <= n + m) T <<= 1;
    vector<Complex> F(T), G(T);
    
    for (int i = 0, x; i <= n; i ++) read(x), F[i].a = x;
    for (int i = 0, x; i <= m; i ++) read(x), G[i].a = x;
    FFT(F, 1), FFT(G, 1);
    
    vector<Complex> H(T);
    for (int i = 0; i < T; i++) H[i] = F[i] * G[i];
    FFT(H, -1);

    for (int i = 0; i <= n + m; i++) write((int)round(H[i].a / T), ' ');
    return 0;
}

快速数论变换（FNTT）

阅读前请自行了解 【笔记】阶和原根 | Thy’s Blog。

FNTT 的优势

在了解了 FFT 之后，不难发现虽然算法的时间复杂度达到了 $\mathcal{O}(n \log n)$，但是在具体实现的时候用到了复数和 double 类型。这显然是非常不好的，在 $n$​ 非常大或者值域很大时，会产生精度错误，并且大量的复数乘法给算法产生了大常数。

那么有什么其它的东西可以代替复数和单位根，精度高，运算快，满足竞赛需求，并且允许我们像类似复数运算那样完成 FFT 算法呢？

有的，兄弟有的，它叫做原根，被广泛用于快速数论变换（Fast Number Theoretic Transform）。

原根

这里再提一下原根是什么：在模数域 $\mathbb{Z}_m$ 下，如果有一个数 $g \in \mathbb{Z}_m$ 满足 $g \perp m$ 且 $\delta_m(g) = \varphi(m)$，我们将 $g$ 称之为 $m$​​ 的原根。

我们定义 $\omega_n := g^{\varphi(m) / n}$ 定义模数域下的单位根。

那原根为什么可以用来替代复数域中的单位根呢，它有如下几个性质：

• $\omega_n$ 的最小正周期为 $n$，它可以生成一共 $n$ 个模意义下互不相同的数。

• $\omega_n$ 必然存在逆元，因为 $\omega_n \perp m$，$\omega_n^k$ 的逆元为 $\omega_n^{-k}$。

• $\omega_n$ 具有周期性质，即 $\omega_n^k = \omega_n^{k + n}$，因此 $\omega_n^{-k} = \omega_n^{n - k}$。

• 具有迭代传递性质，即 $\omega_n^k = \omega_{2n}^{2k}$。

• 中心对称性质：

  $$\omega_n^n \equiv \omega_n^{n / 2}\omega_n^{n / 2} \equiv 1 \pmod m \Rightarrow \omega_n^{n / 2} \equiv -1 \pmod m \Rightarrow \omega_n^{k + n / 2} \equiv -\omega_n^k \pmod m \Rightarrow \omega_n^k = -\omega_n^{k + n / 2}$$