[多项式学习笔记]快速傅里叶变换（FFT）与快速数论变换（NTT）

快速傅里叶变换（FFT）与快速数论变换（NTT）

卷积是多项式的基本运算，而快速傅里叶变换（FFT）可以加速卷积的过程。由于 FFT 引入了三角函数，因此在计算的过程中会产生浮点误差。为了规避误差，我们提出了以数论为基础的快速数论变换（NTT）。它以 FFT 作为基础，但是用原根代替了复平面上的单位根。在只涉及整数的多项式卷积中，我们常常使用 NTT 代替 FFT。

本文首先介绍了卷积的概念，然后讲解 FFT 的原理及实现，再进一步讲解 NTT。

卷积

卷积是一种运算。对于两个数列 \(a\) 和 \(b\)，定义它们的卷积 \(a \ast b\) 为一个新序列 \(c\)，满足 \(c_{k} = \sum_{i = 0}^{k} a_{i}b_{k - i}\)。也就是说，从 \(a\) 和 \(b\) 中各选择一个元素，如果它们的下标之和为 \(k\)，则它们的积就会被加到 \(c_{k}\) 中。

我们之所以定义这种运算，是因为它有明显的组合意义。考虑这样的例子：我们有两个骰子，掷出骰子之后可以得到 \(1 \sim 6\) 中的一个数。由于骰子的质量不均匀，所以得到每个数的概率不相同。具体而言我们分别用 \(a_i\) 和 \(b_i\) 来表示两个骰子掷出 \(i\) 的概率，用 \(c_{k}\) 表示同时掷出两个骰子，得到的数之和为 \(k\) 的概率。那么显然有 \(c_{k} = \sum_{i = 1}^{k - 1} a_{i}b_{k - i}\)——这个式子和卷积的形式相同！也就是说，\(c\) 实际上就是 \(a\) 和 \(b\) 的卷积。（实际上两者求和的上下限有点差别，因为 \(a\) 和 \(b\) 的下标都从 \(1\) 开始。但如果把 \(a_{0}\) 和 \(b_{0}\) 都视作 \(0\)，并修改求和上下限为 \(0\) 和 \(k\)，则两者形式完全相同。）

类似“掷两个骰子”这样的问题产生的组合意义非常常见：本质上来说，如果可以用两种方式完成某个事情，用第一种方式产生 \(i\) 的贡献的方案数为 \(a_i\)，第二种方式为 \(b_i\)，那么两个序列的卷积中第 \(k\) 项 \(c_k\) 就是接连使用两种方式，产生的总贡献为 \(k\) 的方案数。这也就解释了卷积运算为什么重要。

卷积也可以定义在多项式函数上：两个多项式 \(f(x)\) 和 \(g(x)\) 的卷积也是多项式，记为 \((f \cdot g)(x)\)——这实际上就是多项式的乘积。具体来说，给定 \(n\) 次多项式 \(f(x) = \sum_{i = 0}^{n} a_{i}x^{i}\) 和 \(m\) 次多项式 \(g(x) = \sum_{i = 0}^{m} b_{i}x_{i}\)，则两者的卷积为：

\[(f \cdot g)(x) = \sum_{i = 0}^{n} \sum_{j = 0}^{m} a_{i}b_{j}x^{i + j} = \sum_{i = 0}^{n + m} c_{i}x^{i} \]

从序列过渡到多项式是自然的。因为如果我们把序列 \(a\) 中的第 \(i\) 项 \(a_{i}\) 作为多项式的 \(x^{i}\) 项系数，那么一个数列就对应唯一一个多项式。在解题过程中，很多时候多项式不会直接出现，但我们可以把序列转化成多项式。

以上是卷积最常见的形式，但我们还有其它的卷积。把卷积的形式一般化：对两个数列 \(a\) 和 \(b\)，定义其卷积中第 \(k\) 项为

\[c_{k} = \sum_{i \oplus j = k} a_{i}b_{j} \]

其中 \(\oplus\) 表示整数上的二元运算。如果 \(\oplus\) 是加法，这就是上文中卷积的形式。为了区分，也可以称作加法卷积。也可以把 \(\oplus\) 改为其它运算，常见的有减法、\(\max\)、位运算（与/或/异或）等。但无论怎么改，其本质都是：从 \(a\) 和 \(b\) 中各选择一个元素，如果它们下标的 \(\oplus\) 运算结果为 \(k\)，则把它们的卷积加入到 \(c_{k}\) 中。

显然，对于上文提到的这些运算，都可以暴力地在 \(O(n^{2})\) 时间内计算卷积。（这里认为两个序列的长度数量级相同，都视作 \(n\)。）对于不同的卷积，有不同的优化方法。本文介绍的 FFT 用于加速加法卷积，可以在 \(O(n \log n)\) 的时间内计算两个序列进行加法卷积的结果（或者说计算两个多项式的乘积，这两者本质相同。）而减法卷积可以在 \(O(n)\) 时间内转化成加法卷积的形式，因此也可以使用 FFT 优化。

多项式的点值表示

一个 \(n\) 次多项式可以由它的 \((n + 1)\) 项系数唯一确定，这种表示方法称为多项式的系数表示。除此之外多项式还有一种表示方法：点值表示。我们有如下结论：

已知一个 \(n\) 次多项式在 \((n + 1)\) 个不同的点的点值，则这个多项式可以唯一确定。

使用拉格朗日插值可以构造出这样的一个多项式，证明了其存在性。 而唯一性可以由 代数基本定理 导出：假设有两个 \(n\) 次多项式 \(p(x)\) 和 \(q(x)\) 都经过这 \(n + 1\) 个点，则它们的差 \(p(x) - q(x)\) 是不超过 \(n\) 次的多项式，且有至少 \(n + 1\) 个零点。根据代数基本定理，非零的 \(n\) 次多项式只有不超过 \(n\) 个零点，所以 \(p(x) - q(x)\) 必须是零多项式，也就是说 \(p(x) = q(x)\)。

据此结论，对于一个 \(n\) 次多项式 \(f\)，可以任选其图象上 \(n + 1\) 个不同的点来表示它。记 \(\tau(f)\) 为 \(f\) 的点值表示，则 \(f(x)\) 可以被写作

\[\tau(f) = \{(x_{0}, f(x_{0})), (x_{1}, f(x_{1})), \cdots, (x_{n}, f(x_{n})) \} \]

其中 \(x_{0..n}\) 两两不同。

点值表示的优越性在于：如果已知 \(f(x)\) 和 \(g(x)\) 的点值表示，则可以在 \(O(n)\) 的时间复杂度内计算二者的卷积（而非系数表示法的 \(O(n^{2})\)）。假设二者的点值表示分别为

\[\tau(f) = \{(x_{0}, f(x_{0})), (x_{1}, f(x_{1})), \cdots, (x_{n + m}, f(x_{n + m})) \} \\ \tau(g) = \{(x_{0}, g(x_{0})), (x_{1}, g(x_{1})), \cdots, (x_{n + m}, g(x_{n + m})) \} \]

则

\[\tau(f \cdot g) = \{(x_{0}, f(x_{0})g(x_{0})), (x_{1}, f(x_{1})g(x_{1})), \cdots, (x_{n + m}, f(x_{n + m})g(x_{n + m})) \} \]

正确性显然。

需要注意的是，由于 \((f \cdot g)(x)\) 是 \((n + m)\) 次多项式，所以需要 \((n + m + 1)\) 个点的点值来表示它。为此，\(f\) 和 \(g\) 都需要 \((n + m + 1)\) 个点值来表示，尽管它们原本的次数可能小于 \(n + m\)。但这并没有什么问题，因为把一个多项式的次数提高，可以看作高次项的系数都是 \(0\)，不改变多项式的值。

现在我们看到了优化卷积的时间复杂度的一种方案：由于在点值表示法下，计算两个多项式的卷积非常快，如果能快速地把一个多项式的系数表示法转换成点值表示法，再将点值表示法转换成系数表示法，就可以优化卷积的时间复杂度了。为了快速地把一个多项式在系数表示法和点值表示法之间转换，我们还需要引入一个概念：

单位根

在复平面上有一个单位圆，把它 \(n\) 等分，每个 \(n\) 等分点可以看作一个复数，统称为 \(n\) 次单位根，因为它们都是方程 \(\omega^{n} = 1\) 的解。

记辐角为正且最小的 \(n\) 次单位根为主 \(n\) 次单位根 \(\omega_{n}\)，可以计算出 \(\omega_{n} = \cos \dfrac{2\pi}{n} + i\sin \dfrac{2\pi}{n}\)，但除了在代码实现中需要用到，我们并不关心它具体的值。其它单位根都可以用 \(\omega_{n}\) 表示：逆时针方向，第 \(k\) 个 \(n\) 次单位根为 \(\omega_{n}^{k}\)。因为 \(\omega_{n}\) 的辐角为 \(\theta = 2k\pi/n\)，则一个复数乘上 \(\omega\) 就相当于把这个复数旋转了 \(\theta\)，旋转几次就能得到第几个单位根。

我们需要用到单位根的如下几个性质：

  1. （消去引理）\(\omega_{rn}^{rk} = \omega_{n}^{k}\)。相当于可以同时“约掉”单位根的上下指标，实际上它和分数的约分也确实有关系，因为

    \[\omega_{rn}^{rk} = \cos \dfrac{2rk\pi}{rn} + i\sin \dfrac{2rk\pi}{rn} = \omega_{n}^{k} \]

    • 推论：\(\omega_{n}^{k + n / 2} = -\omega_{n}^{k}\)。这是因为 \(\omega_{n}^{n / 2} = -1\)。

  2. （求和引理）

    \[\sum_{i = 0}^{n - 1} (\omega_{n}^{k})^{i} = \begin{cases} 0, n \not\mid k \\ n, n \mid k \end{cases} \]

    当 \(n \mid k\) 时，\(\omega_{n}^{k} = 1\)，因此和式的值为 \(n\)。否则直接代入等比数列求和公式即可证明：

    \[\sum_{i = 0}^{n - 1} (\omega_{n}^{k})^{i} = \dfrac{1 - \omega_{n}^{nk}}{1 - \omega_{n}^{k}} = 0 \]

除此之外，定义本原单位根为满足 \(\gcd(k, n) = 1\) 的单位根 \(\omega_{n}^{k}\)。例如 \(\omega_{n}\) 本身就是本原单位根。之所以定义本原单位根，是因为它们满足一个重要的性质：对任意 \(0 \le i < n\)，\((\omega_{n}^{k})^{i}\) 的取值两两不同。从数论的角度来说，就是要证明对 \(0 \le i < j < n\)，\(ik \not\equiv jk \pmod{n}\)。由于 \(\gcd(k, n) = 1\)，可以直接消去 \(k\)，得到 \(i \not\equiv j \pmod{n}\)，这显然成立。在 FFT 的过程中会看到这个性质的重要性。

系数表示法转点值表示法（DFT）

有了这些前置知识后就可以进入正题了。快速卷积的算法可以分为三步：

  1. 把 \(f(x)\) 和 \(g(x)\) 转换成点值表示法 \(\tau(f)\) 和 \(\tau(g)\)。（DFT）
  2. 计算 \(\tau(f \cdot g)\)。
  3. 把 \(\tau(f \cdot g)\) 转成系数表示法，从而得到答案 \((f \cdot g)(x)\)。（IDFT）

第一步称为离散傅里叶变换（DFT），FFT 是一种快速实现 DFT 的算法。

现在我们想要计算 \(n\) 项多项式 \(f(x)\) 的点值表示法 \(\tau(f)\)。（注意 \(n\) 是项数而非次数，和上文不同）需要保证 \(n\) 是 \(2\) 的非负整数次幂，如果不是，则补上若干个系数为 \(0\) 的高次项。

我们需要快速计算 \(f(x)\) 在 \(n\) 个点的点值。下面将会看到，如果选择计算 \(f(x)\) 在 \(\omega_{n}^{0}, \omega_{n}^{1}, \omega_{n}^{2}, \cdots, \omega_{n}^{n - 1}\) 这 \(n\) 个点处的点值，由于单位根的良好性质，可以在 \(O(n \log n)\) 的时间内得到结果。

把 \(f(x)\) 的奇数次项和偶数次项的系数提取出来，构建两个新的多项式：

\[f_{0}(x) = a_{0} + a_{2}x + a_{4}x^{2} + \cdots + a_{n - 2}x^{n / 2 - 1} \\ f_{1}(x) = a_{1} + a_{3}x + a_{5}x^{2} + \cdots + a_{n - 1}x^{n / 2 - 1} \]

则 \(f(x) = f_{0}(x^{2}) + xf_{1}(x^{2})\)。

代入 \(x = \omega_{n}^{k}\)（\(k < n / 2\)）可得

\[\begin{aligned} f(\omega_{n}^{k}) &= f_{0}(\omega_{n}^{2k}) + \omega_{n}^{k} f_{1}(\omega_{n}^{2k}) \\ &= f_{0}(\omega_{n / 2}^{k}) + \omega_{n}^{k} f_{1}(\omega_{n / 2}^{k}) \end{aligned} \]

代入 \(x = \omega_{n}^{k + n / 2}\)（\(k < n / 2\)）可得

\[\begin{aligned} f(\omega_{n}^{k + n / 2}) &= f_{0}(\omega_{n}^{2k + n}) + \omega_{n}^{k + n / 2} f_{1}(\omega_{n}^{2k + n}) \\ &= f_{0}(\omega_{n}^{2k} \cdot \omega_{n}^{n}) - \omega_{n}^{k} f_{1}(\omega_{n}^{2k} \cdot \omega_{n}^{n}) \\ &= f_{0}(\omega_{n / 2}^{k}) - \omega_{n}^{k} f_{1}(\omega_{n / 2}^{k}) \end{aligned} \]

（这里体现了为什么 FFT 中多项式的项数必须是 \(2\) 的非负整数次幂：只有当 \(n\) 是偶数时，才能恰好划分成两个大小相等的子问题，这样带来了优美的结构。为了每次递归时 \(n\) 都是偶数，必须一开始就让多项式的项数必须是 \(2\) 的非负整数次幂。）

可以发现，如果知道 \(f_{0}(\omega_{n / 2}^{k})\) 和 \(f_{1}(\omega_{n / 2}^{k})\)，就可以在 \(O(1)\) 的时间内计算 \(f(x)\) 在 \(x = \omega_{n}^{k}\) 和 \(x = \omega_{n}^{k + n / 2}\) 的点值。要对 \(0 \le k < n\) 求出所有 \(f(x)\) 在 \(x = \omega_{n}^{k}\) 的点值，就需要对所有 \(0 \le k < n / 2\) 求出所有 \(f_{0}(x)\) 和 \(f_{1}(x)\) 在 \(x = \omega_{n / 2}^{k}\) 的点值。

这就把原问题转化成了两个规模减小了一半的子问题：我们原本要求出 \(n\) 项多项式 \(f(x)\) 在 \(n\) 个点的点值，现在要求出两个 \(n / 2\) 项多项式 \(f_{0}(x)\) 和 \(f_{1}(x)\) 在 \(n / 2\) 个点的点值。\(n = 1\) 是边界条件：此时有 \(f(\omega_{1}^{0}) = f(1) = a_{0}\)。

按照上述过程递归计算就可以得到待求多项式的点值表示法。可以用递归式 \(T(n) = 2T(n / 2) + O(n)\) 表示时间复杂度，解得 \(T(n) = O(n \log n)\)。

点值表示法转系数表示法（IDFT）

进行 DFT 之后，已知多项式 \(f(x)\) 的点值表示法（其中 \(n\) 是 \(2\) 的非负整数次幂）

\[\tau(f) = \{(\omega_{n}^{0}, y_{0}), (\omega_{n}^{1}, y_{1}), \cdots, (\omega_{n}^{n - 1}, y_{n - 1})\} \]

现在我们要求出它的系数表示法。

从点值反推多项式的系数，这个过程就是插值。我们已经知道拉格朗日插值可以在 \(O(n^{2})\) 的时间内解决这个问题，但还不够快。为了加速这个过程，我们仍然要用到单位根的性质。这个过程称为逆离散傅里叶变换（IDFT）。

这里直接给出结论：把点值数组代入 DFT ，并把 \(\omega_{n}\) 替换成 \(\omega_{n}^{-1}\)，然后把结果数组中的每一项都除以 \(n\)，就能得到 \(f(x)\) 的系数表示法。这就是 IDFT 的过程。

证明如下：

构造多项式

\[F(x) = \sum_{i = 0}^{n - 1} y_{i}x^{i} \]

我们要表示出它在 \(\omega_{n}^{0}, \omega_{n}^{-1}, \omega_{n}^{-2}, \cdots, \omega_{n}^{-(n - 1)}\) 处的点值。

代入得

\[\begin{aligned} F(\omega_{n}^{-k}) &= \sum_{i = 0}^{n - 1} y_{i} \omega_{n}^{-ik} \\ &= \sum_{i = 0}^{n - 1} f(\omega_{n}^{i}) \omega_{n}^{-ik} \\ &= \sum_{i = 0}^{n - 1} \left( \omega_{n}^{-ik} \cdot \sum_{j = 0}^{n - 1} a_{j}\omega_{n}^{ij} \right) \\ &= \sum_{i = 0}^{n - 1} \sum_{j = 0}^{n - 1} a_{j}\omega_{n}^{i(j - k)} \\ &= \sum_{j = 0}^{n - 1} a_{j} \sum_{i = 0}^{n - 1} (\omega_{n}^{j - k})^{i} \end{aligned} \]

根据单位根的求和引理，只有在 \((j - k) \mid n\) 时，\(\sum_{i = 0}^{n - 1} (\omega_{n}^{j - k})^{i}\) 才不为 \(0\)。唯一使 \((j - k) \mid n\) 成立的情况是 \(j = k\)，因此

\[F(\omega_{n}^{-k}) = n \cdot a_{k} \]

也就是说，\(F(x)\) 在 \(\omega_{n}^{0}, \omega_{n}^{-1}, \omega_{n}^{-2}, \cdots, \omega_{n}^{-(n - 1)}\) 处的点值表示法为

\[\tau(F) = \{(\omega_{n}^{0}, n \cdot a_{0}), (\omega_{n}^{-1}, n \cdot a_{1}), \cdots, (\omega_{n}^{-(n - 1)}, n \cdot a_{n - 1})\} \]

这个形式印证了一开始说的：用 \(\omega_{n}^{-1}\) 替代 \(\omega_{n}\) 进行 DFT，把结果数组中的每一项都除以 \(n\)，就能得到 \(f(x)\) 的系数表示法。下面进一步说明其正确性。

注意到 DFT 的正确性在于单位根所具有的性质（例如满足消去引理等），而 \(\omega_{n}^{-1}\) 也是单位根，也满足这些性质。除此之外，\(\omega_{n}\) 还是本原单位根，这保证了它的 \(0\) 到 \(n - 1\) 次幂两两不同，所以不会产生相同的点值。而 \(\omega_{n}^{-1} = \omega_{n}^{n - 1}\)，\(n\) 和 \(n - 1\) 显然互质，所以 \(\omega_{n}^{-1}\) 也是本原单位根。既然 \(\omega_{n}^{-1}\) 满足 \(\omega_{n}\) 所具有的所有重要性质，用 \(\omega_{n}^{-1}\) 进行 IDFT 也就是正确的，

FFT 的递归实现

根据现有的推导，我们可以写出 DFT 和 IDFT 的递归实现版本。我们用 C++ 自带的复数类（complex）来存储多项式的系数、点值和单位根。

我们用欧拉公式计算单位根：\(\omega_{n} = e^{\frac{2\pi}{n}} = \cos \dfrac{2\pi}{n} + i\sin \dfrac{2\pi}{n}\)。计算 \(\omega_{n}^{-1}\) 时只需把 \(\sin\) 前的正号改成负号。在函数中传入一个值代表 \(\sin\) 前的符号，就可以用一个函数同时实现 DFT 和 IDFT。

下面的代码实现中，传入的 s 数组代表系数（在 IDFT 中则是点值），求出点值（在 IDFT 中是系数）后原址存储在 s 数组中。

void solve(vector<complex<double>> &s, int n, int o) {
    if(n == 1) return;
    vector<complex<double>> s0(n / 2), s1(n / 2);
    for(int i = 0; i < n; i += 2) {
        s0[i / 2] = s[i], s1[i / 2] = s[i + 1]; // 提取奇数次项和偶数次项
    }
    solve(s0, n / 2, o), solve(s1, n / 2, o);
    complex<double> wi(1, 0), wn(cos(2 * pi / n), o * sin(2 * pi / n));
    for(int i = 0; i < n / 2; i++) { // 合并结果
        s[i] = s0[i] + wi * s1[i];
        s[i + n / 2] = s0[i] - wi * s1[i];
        wi *= wn;
    }
}

void FFT(vector<complex<double>> &s, int n) { solve(s, n, 1); }
void IFFT(vector<complex<double>> &s, int n) { solve(s, n, -1); }

更高效的 FFT 实现

迭代实现（蝶形变换）

递归调用会产生额外的时间开销，我们希望优化这一点。下面介绍 FFT 的迭代实现。

以 \(n = 8\) 为例，分析 FFT 递归树的结构：

（此处省略一张图）

如果能把系数按其出现在叶子节点中的顺序重新排布，就可以自底向上地合并结果。

为了得到每个系数在叶子节点中第几个出现，我们应用位逆序置换。我们有如下结论：把 \(i\) 看作 \(\log_2(n)\) 位二进制数，系数 \(a_{i}\) 在叶子节点中的下标为其二进制表示的翻转。例如 \(n = 8\) 时，\(5 = (011)_{2}\) 最终的位置是 \((110)_{2} = 6\)。

证明：设 \(p_i\) 是 \(i\) 最终的位置。一开始令 \(p_i = 0\)，然后逐步计算 \(p_i\)。如果 \(i\) 二进制下的最低位为 \(0\)，则 \(i\) 被分到左半部分，对应 \(p_i\) 最高位为 \(0\)；否则被分到右半部分，对应 \(p_i\) 最高位为 \(1\)。依次类推，每次递归都相当于把 \(i\) 的最低位移到到 \(p_i\) 的下一位中，最终递归结束时 \(p_i\) 就是 \(i\) 二进制表示的翻转。

记 \(R(x)\) 表示把 \(n\) 位二进制数 \(x\) 翻转得到的结果，显然可以在 \(O(n \log n)\) 的时间内求出，但有简单的 \(O(n)\) 递推方法：

\[R(x) = \begin{cases} 0, x = 0 \\ \lfloor R(\lfloor x / 2\rfloor) / 2 \rfloor + (x \bmod 2) \cdot (n / 2) \end{cases} \]

我们考虑把 \(x\) 分成 \(\lfloor x / 2 \rfloor\) 和 \(x \bmod 2\) 两部分，分别翻转再合并结果。前者翻转的结果是 \(R(\lfloor x / 2 \rfloor)\)，但由于 \(\lfloor x / 2 \rfloor\) 的最高位为 \(0\)，所以翻转后在末尾多了一个 \(0\)，应该去掉，对应 \(\lfloor R(\lfloor x / 2\rfloor) / 2 \rfloor\)。再加上 \((x \bmod 2) \cdot (n / 2)\) 即可得到递推式。

应用位逆序置换后，自底向上合并结果，即可实现 FFT。

三次变两次

设 \(f(x) = \sum_{k} a_{k} x^{k}\)，\(g(x) = \sum_{k} b_{k} x^{k}\)。要求两者的卷积，可以用以下方法把原先调用三次 FFT 优化到只调用两次 FFT：

构造多项式 \(h(x) = f(x) + ig(x)\)，则

\[h^{2}(x) = (f(x) + ig(x))^{2} = (f^{2}(x) - g^{2}(x)) + 2(f(x) + g(x))i \]

把 \(h^{2}(x)\) 的 \(k\) 次项系数的虚部除以二，就可以得到 \((f \cdot g)(x)\) 的 \(k\) 次项系数。因此，我们直接把 \(f(x)\) 和 \(g(x)\) 分别放到 \(h(x)\) 的实部和虚部中，再平方，就可以得到 \(f(x)\) 和 \(g(x)\) 的卷积。这样只需做一次 DFT，因此把 FFT 的总次数降到了两次。

需要注意的是，虽然这种方法减小了常数，但是增大了误差，特别是在计算数本身就比较大的时候，需要谨慎使用。（例如 「ZJOI2014」力 这题，使用两次 FFT 连样例都过不去。）

速度比较

实现版本 时间/ms
递归版本，三次 FFT 1638
迭代版本，三次 FFT 872
迭代版本，两次 FFT 537
迭代版本，三次 NTT 992

无精度损失的多项式乘法：快速数论变换（NTT）

FFT 最大的缺点在于计算单位根时引入了浮点运算。即使多项式的系数都是整数，用 FFT 计算时还是会产生误差。为了规避这一点，我们引入了基于数论的快速数论变换（NTT），它的思想和 FFT 本质相同。只是把复数单位根替换成了原根。

本质上来说，FFT 要找一个数 \(z\)，满足 \(z^{n} = 1\) 且当 \(0 \le k < n\) 时 \(z^{i}\) 两两不同。在复数域中，我们已经知道 \(n\) 次本原单位根就满足这样的性质。在整数中我们也想找到这样的值。如果不取模，这样的值显然是不存在的，所以下面我们只考虑模 \(p\) 的数域。

为了找到这个值，下面简要介绍一些数论概念：

阶

若正整数 \(a\)，\(p\) 互质且 \(p > 1\)，定义 \(a\) 模 \(p\) 的阶为：最小的整数 \(n\)，满足 \(a^{n} \equiv 1 \pmod{p}\)，记作 \(\delta_{p}(a)\)。

阶的性质：对于 \(0 \le i < \delta_{p}(a)\)，\(a^{i} \bmod p\) 两两不同。（注意，这就是我们要找的性质！）

证明：假设存在 \(0 \le j < k < \delta_{p}(a)\) 满足 \(a^{j} \equiv a^{k} \pmod{p}\)，则 \(a^{k - j} \equiv 1 \pmod{p}\)，与阶的定义矛盾。

原根

若 \(\delta_{p}(g) = \varphi(p)\)，则称 \(g\) 是模 \(p\) 的原根。注意，原根可能有多个，也可能不存在。

定理：对于素数 \(p\)，一定存在模 \(p\) 的原根。证明略。

NTT 的实现

在 NTT 的过程中，模数通常是形如 \(p = a \cdot 2^{k} + 1\) 的质数，且 \(2^{k}\) 需要比 FFT 中多项式的长度要大。这么做是为了保证存在我们要找的那个整数：由于 \(p\) 是素数，所以一定存在原根 \(g\)，满足 \(g^{\varphi(p)} 1 \equiv \pmod{p}\)。又由于 \(\varphi(p) = p - 1\)，所以 \(g^{p - 1} \equiv 1 \pmod{p}\)。由于 \(p = a \cdot 2^{k} + 1\)，所以 \(p - 1\) 是 \(2\) 的非负整数次幂的倍数。设 \(n\) 是 FFT 中多项式的长度，定义 \(g_{n} = g^{(p - 1) / n}\)，那么 \(g\) 就是我们寻求的用来替代复数单位根的数：它满足 \(g_{n}^{n} = g^{p - 1} \equiv 1 \pmod{p}\)，再根据阶的定义，对 \(0 \le k < n\)，\(g_{n}^{k} \bmod p\) 两两不同。

因此，在 DFT 的过程中把 \(\omega_{n}\) 替换成 \(g_{n}\) 就可以得到 NTT。而 \(\omega_{n}^{-1}\) 之于 \(\omega_{n}\) 就相当于 \(g_n^{-1}\) 之于 \(g_{n}\)，所以在 IDFT 的过程中把 \(\omega_{n}^{-1}\) 替换成 \(g\) 的逆元 \(g_{n}^{-1}\) 即可。

以下是几个常见的 NTT 质数及其原根：

  • \(p = 998244353 = 7 \times 17 \times 2^{23} + 1\)，\(g = 3\)。
  • \(p = 1004535809 = 479 \times 2^{21} + 1\)，\(g = 3\)。
  • \(p = 4179340454199820289 = 29 \times 2^{57} + 1\)，\(g = 5\)。

例题

P1919 【模板】高精度乘法 | A*B Problem 升级版

FFT 的直接应用。对于整数的十进制表示 \(\overline{a_{n = 1}\cdots a_{1}a_{0}}\)，把它看作 \(n\) 项多项式 \(f(x) = \sum_{i = 0}^{n - 1} a_{i}x^{i}\) 代入 \(x = 10\) 得到的点值，那么两个整数相乘就可以转化成两个多项式相乘，要求的整数之积就是多项式的积代入 \(x = 10\) 的点值，可以用 FFT 加速这个过程。

求出多项式的积以后，并不需要真正代入 \(x = 10\) 求点值，只需要把它的各项系数进位以后输出即可。（需要注意的是系数不一定小于 \(10\)。）

ABC392G Fine Triplets

这道题看似和多项式没有任何关系，但是可以转化成多项式问题以后用 FFT 加速求解。（实际上应该说是转化成了生成函数问题。）

改写题目中的式子为 \(2B = A + C\)，我们要做的就是对所有整数 \(k\)，求有多少个数对 \((a, b)\) 满足 \(a + b = 2k\)。构造多项式 \(f(x) = \sum_{a \in S} x^{a}\)，设 \(c_{k} = [x^{k}]f^{2}(x)\)，则 \(c_{2k}\) 就是从 \(S\) 中选择两个数，使得和为 \(2k\) 的方案数。（我们可以把多项式的卷积的组合意义看作两个集合的笛卡尔积）。但是这包含了 \(k + k = 2k\) 这种不合法的情况，而且把每种方案统计了两次，在考虑这些影响之后得到答案为 \({\large\sum\limits}_{k} [k \in S]\dfrac{c_{2k} - 1}{2}\)。用 FFT/NTT 加速卷积即可。AC 记录

「ZJOI2014」力

首先

\[E_{j} = \sum_{i = 1}^{j - 1} \dfrac{q_i}{(i - j)^{2}} - \sum_{i = j+ 1}^{n} \dfrac{q_i}{(i - j)^{2}} \]

我们的目的是把式子推成两个多项式的卷积的形式，然后就可以用 FFT 加速。所谓“卷积” 的形式就是：和式中有两项相乘，两项的下标之和为某定值 \(k\)，且和式的下界为 \(0\)，上界为 \(k\)。

分母的 \((i - j)^{2}\) 不好处理，但我们不用处理，直接设 \(g_{i} = 1 / i^{2}\)，则

\[E_{j} = \sum_{i = 1}^{j - 1} q_{i}g_{j - i} - \sum_{i = j+ 1}^{n} q_{i}g_{i - j} \]

左边和式非常接近卷积的形式，只是上下界有点差别。令 \(q_i = g_i = 0\)，那么我们可以直接把上下界写成 \(0\) 和 \(j\) 而不改变结果。

对于右边的和式，为了方便，先把下界改成 \(j\)，这也不影响结果。然后要把求和下界改成 \(0\)。先改写成 \(\sum_{j \le i \le n} q_{i}g_{i - j}\)，为了把求和下界改成 \(0\)，用 \(i + j\) 替换 \(i\)，则和式变成 \(\sum_{j \le i + j \le n} q_{i + j}g_{i}\)，整理得 \(\sum_{i = 0}^{n - j} q_{i + j}g_{i}\)。现在已经是差卷积的形式。设 \(p\) 是 \(q\) 翻转以后的序列，则和式变为 \(\sum_{i = 0}^{n - j} p_{n - i - j}g_{i}\)，这就是 \(p\) 和 \(g\) 卷积的第 \(n - j\) 项。到此就做完了，时间复杂度 \(O(n \log n)\)。AC 记录

需要指出，本题的数比较大，两次 FFT 会产生较大的误差，应该用三次 FFT。

[AH2017/HNOI2017] 礼物

AC 记录

[ICPC 2021 Macao R] Pass the Ball!

题目的样例解释有点骗人。相比于考虑第 \(i\) 个孩子在 \(k\) 轮后拿着第几个球，我们更应该考虑第 \(i\) 个球 \(k\) 轮后在哪个孩子手上，因为研究球传递的轨迹更直观。

记题目给出的排列 \(p\) 为 \(\sigma\)。设第 \(i\) 个球 \(k\) 轮后在第 \(c_{i}\) 个孩子手上，则 \(c_{i} = \sigma^{k}(i)\)， 因此 \(\operatorname{ans}_{k} = {\large\sum}_{i = 1}^{n} i \cdot \sigma^{k}(i)\)。套路性地把排列拆分成置换环。先假设整个排列只有一个置换环。在环上考虑传递小球的过程：想象每个点上有一个球，一开始球的编号和节点的编号相同。每传递一次，所有球就在环上前进一次。任选环上一个点做为起点，记 \(a_{i}\) 表示环上第 \(i\) 个点的编号。那么 \(k\) 轮后编号为 \(a_{i}\) 的球在编号为 \(a_{i + k}\) 的点上。（这里 \(i + k\) 是环意义上的加法：例如 \(n + 1\) 表示从第 \(n\) 个点向前走一步，停在第 \(1\) 个点。此外假设 \(0 \le k < n\)，否则可以直接对 \(n\) 取模）考虑环上的加法比较麻烦，我们断环成链，把 \(a\) 复制一份，就有：

\[\operatorname{ans}_{k} = \sum_{i = 1}^{n} a_{i}a_{i + k} \]

这就是差卷积的形式。记 \(b\) 表示 \(a_{1..n}\) 翻转以后的序列，就可以转成加法卷积：

\[\operatorname{ans}_{k} = \sum_{i = 1}^{n} b_{n - i + 1}a_{i + k} \]

也就是说 \(\operatorname{ans}_{k}\) 就是 \(a\) 和 \(b\) 卷积的第 \(n + 1 + k\) 项。用 NTT 加速卷积即可。需要注意的是本题的答案是 \(O(n^{3})\) 量级，最大约为 \(10^{15}\)，用 FFT 会产生较大的精度误差。这里我们使用一个经典技巧：虽然题目没有让我们取模，但是我们知道答案的上界，如果我们选择一个大于上界的数作为模数，即使计算的过程中超出了模数，最终也能得到正确答案，这是显然的。所以我们选择一个大于 \(10^{15}\) 的质数进行 NTT。（我使用的是 \(p = 27 \times 2^{56} + 1\)，原根 \(g = 5\)。）计算过程中需要用到 __int128。

接下来考虑有多个置换环的情况。这也是简单的：置换环之间互不影响，所以分别计算每个置换环的答案再加起来即可。由于置换环的大小之和为 \(n\)，所以 NTT 的总时间复杂度也为 \(O(n \log n)\)。但是对于一个询问，要枚举所有的置换环来统计，置换环较多的时候时间复杂度为 \(O(nq)\)。想想我们是怎么统计答案的：对于长度为 \(m\) 的置换环，询问第 \(k\) 轮时，我们在答案中加上卷积的第 \(m + 1 + (k \bmod m)\) 项系数。对于相同的 \(m\)，加上的是卷积中的相同项的系数。这启示我们把相同大小的置换环的贡献加在一起处理。又由于置换环的大小之和为 \(n\)，根据经典结论，大小不同的置换环个数只有 \(O(\sqrt{n})\) 个，所以可以做到单次询问 \(O(\sqrt{n})\)。

综上所述，我们在 \(O(n \log n + q \sqrt{n})\) 的时间内解决了此问题。AC 记录

关于大小不同置换环的个数只有 \(O(\sqrt{n})\) 的证明：

\[\sum_{k = 1}^{2\lceil \sqrt{n} \rceil} k > \sum_{k = \lceil \sqrt{n} \rceil + 1}^{2\lceil \sqrt{n} \rceil} > \lceil \sqrt{n} \rceil^{2} \ge n \]

参考资料

posted @ 2025-03-31 18:12  DengStar  阅读(342)  评论(0)    收藏  举报