当前位置：首页 > news >正文

NTT快速数论变换：从原理到C++实现与优化

news 2026/7/16 11:47:52

1. 项目概述：为什么我们需要NTT？

在算法竞赛和高性能计算领域，多项式乘法是一个绕不开的核心问题。无论是大整数乘法、字符串匹配，还是生成函数相关的计数问题，最终都归结为两个多项式系数的卷积运算。最直观的算法是 O(n²) 的双层循环，这在数据规模稍大时（比如 n=10⁵）就完全不可行了。

这时，快速傅里叶变换（FFT）登场了。它利用复数单位根的性质，将卷积运算的复杂度降到了 O(n log n)，堪称魔法。然而，FFT 有一个天生的“阿喀琉斯之踵”：浮点数精度误差。当你处理整数系数、特别是模意义下的运算时，FFT 计算过程中的舍入误差可能会在取整时导致错误结果。虽然可以通过调整 eps 或使用 long double 来缓解，但在需要绝对精确的场合（比如模 998244353 下的卷积），这始终是个隐患。

于是，快速数论变换（NTT）应运而生。它可以说是 FFT 在整数域上的“孪生兄弟”，将复数域上的单位根替换为有限域（通常是模一个大质数 p 的剩余系）下的“原根”，从而在保持 O(n log n) 时间复杂度的同时，实现了完全精确、无误差的整数卷积。对于 ACM-ICPC、Codeforces 等竞赛选手，或是任何需要在模意义下进行快速多项式运算的开发者来说，掌握 NTT 是一项必备的硬核技能。

这篇文章，我将从一个 C++ 实现者的角度，带你彻底吃透 NTT。我们不只讲“怎么做”，更要讲清楚“为什么这么做”，并分享我在实战中积累的模板代码、调试技巧和性能优化心得。无论你是正在备赛的 OIer，还是对高性能数学计算感兴趣的工程师，这篇文章都能让你从原理到实现，一站式搞定 NTT。

2. 核心原理：从FFT到NTT的数学桥梁

要理解 NTT，必须先理解 FFT。FFT 的核心思想是“分治”，而分治得以进行的关键，在于单位根ω_n = e^(2πi/n) 所具备的一系列完美性质：

消去引理：ω_{dn}^{dk} = ω_n^k
折半引理：ω_n^{k + n/2} = -ω_n^k （当 n 为偶数时）
求和引理：Σ_{j=0}^{n-1} (ω_n^k)^j 在特定条件下为 n 或 0。

这些性质使得我们可以将规模为 n 的 DFT 分解为两个规模为 n/2 的 DFT，从而实现 O(n log n) 的快速变换。

NTT 的精妙之处在于，它在模素数 p 的剩余系中，找到了具有完全类似性质的数——原根。

2.1 原根与单位根的对应关系

设 p 是一个素数。根据数论中的原根定理，存在一个整数 g，使得 g, g², …, g^{p-1} 在模 p 下遍历 1 到 p-1 的所有整数。这个 g 就是模 p 的一个原根。

现在，假设我们需要进行长度为 n 的 NTT，并且 n 是 2 的幂次（这是 Cooley-Tukey 蝶形算法的要求）。NTT 能够成立的一个关键前提是：n 必须能整除 p-1。因为只有这样，我们才能在模 p 的剩余系中，找到一个阶恰好为 n 的元素。

具体操作如下：

令gn = g^((p-1)/n) mod p。
那么，序列{gn^0, gn^1, gn^2, ..., gn^(n-1)}在模 p 下，就构成了一个“n 次单位根群”。
这个gn完美模拟了 FFT 中 ω_n 的角色：
- gn^n ≡ 1 (mod p)（对应 ω_n^n = 1）
- gn^(n/2) ≡ -1 (mod p)（对应 ω_n^{n/2} = -1，这是折半引理的关键）
- gn^(k + n/2) ≡ -gn^k (mod p)（折半引理）
- 消去引理同样成立。

这样一来，FFT 中所有基于复数单位根的推导和蝶形运算结构，都可以原封不动地移植到模 p 的整数世界，只需要把复数乘法换成模 p 乘法，把复数加法换成模 p 加法。这就是 NTT 的理论基础。

2.2 模数的选择：为什么是998244353？

你会在几乎所有竞赛代码中看到P = 998244353。这不是偶然，而是精心选择的结果。一个好的 NTT 模数 p 需要满足以下几个条件：

p 是素数：这是原根存在的前提。
p-1 包含大量因子 2：因为我们的变换长度 n 必须是 2 的幂。p-1中因子 2 越多，我们能进行的最大变换长度n_max就越大。998244353 = 119 * 2^23 + 1，这意味着它最多支持长度为2^23的 NTT，足以应对绝大多数题目。
p 的大小适中：既要足够大，使得两个模数内的数相乘不会溢出常用的整数类型（如 64 位），又要便于计算。998244353约等于 1e9，其平方约 1e18，在 64 位有符号整数范围内，乘法结果可以用1ll * a * b % P安全计算。
存在一个较小的原根：g = 3就是 998244353 的一个原根。计算gn = pow_mod(g, (p-1)/n, p)非常方便。

其他常见的 NTT 友好模数还有：

1004535809 = 479 * 2^21 + 1,g = 3
469762049 = 7 * 2^26 + 1,g = 3
167772161 = 5 * 2^25 + 1,g = 3

当结果可能超过单个模数的范围时，可以使用任意模数 NTT（MTT），通常通过使用两个或多个上述模数分别进行 NTT，然后用中国剩余定理（CRT）合并结果。这是后话，但了解模数家族很有必要。

实操心得：在写代码时，我习惯将模数P和原根G定义为全局常量。对于 998244353，它的一个常用原根是 3，而(P-1)的最大 2 的幂因子是1 << 23。记住这个“23”，它在确定最大变换长度时很有用。

3. 算法实现：迭代Cooley-Tukey蝶形算法详解

理解了原理，我们来看如何用 C++ 实现一个高效、实用的 NTT。我们将采用迭代版的 Cooley-Tukey 算法，它比递归版更快，且避免了递归开销。

3.1 核心操作：蝶形变换

蝶形变换是 FFT/NTT 分治过程的核心操作。对于 NTT，其蝶形变换公式如下：

// 已知当前单位根 gn，对位置 i 和 i+k 进行操作 (k 是当前半长度) u = a[i]; v = 1ll * a[i + k] * gn % P; // 注意这里用 1ll 防止乘法溢出 a[i] = (u + v) % P; a[i + k] = (u - v + P) % P; // 保证结果非负

这个操作将两个点的值线性组合，是构成整个快速变换的基本单元。

3.2 位逆序置换

迭代算法需要先将输入序列按照“位逆序”重新排列。什么是位逆序？对于一个长度为 n（2 的幂）的序列，下标 i 的二进制表示反转后，就得到了它应该放置的新位置r[i]。例如，n=8 时：i=1 (001)的逆序是4 (100)。i=3 (011)的逆序是6 (110)。

这个预处理是为了让迭代过程能够正确地“合并”子问题的解。有一个高效递推公式可以计算r[i]：

for (int i = 0; i < lim; ++i) { r[i] = (r[i >> 1] >> 1) | ((i & 1) ? (lim >> 1) : 0); }

其中lim是变换长度 n。这个公式的意思是：i的逆序等于i/2的逆序右移一位，然后根据i的最低位决定是否补上最高位的 1。

3.3 完整迭代NTT过程

结合蝶形变换和位逆序置换，我们可以写出正向 NTT（求值）的过程：

const int P = 998244353, G = 3; // 模数和原根 // 快速幂，用于计算逆元 int qpow(int x, int y) { int res = 1; while (y) { if (y & 1) res = 1ll * res * x % P; x = 1ll * x * x % P; y >>= 1; } return res; } void ntt(int *a, int lim, int opt) { // opt=1 为正变换，opt=-1 为逆变换 // 1. 位逆序置换 for (int i = 0; i < lim; ++i) { if (i < r[i]) swap(a[i], a[r[i]]); // 每个对只交换一次 } // 2. 迭代进行蝶形变换 for (int m = 2; m <= lim; m <<= 1) { // m 是当前合并的子问题规模 int k = m >> 1; // 半长度 // 计算当前层的主单位根：gn = g^((P-1)/m) // 如果 opt==-1（逆变换），则使用 gn 的逆元 int gn = qpow(G, (P - 1) / m); if (opt == -1) gn = qpow(gn, P - 2); // 费马小定理求逆元 for (int i = 0; i < lim; i += m) { // 遍历每一块 int g = 1; // 当前单位根的幂 for (int j = 0; j < k; ++j) { // 对块内的每一对进行操作 int u = a[i + j]; int v = 1ll * a[i + j + k] * g % P; a[i + j] = (u + v) % P; a[i + j + k] = (u - v + P) % P; g = 1ll * g * gn % P; // 更新到下一个单位根 } } } // 3. 如果是逆变换，需要除以 n if (opt == -1) { int inv_lim = qpow(lim, P - 2); // lim 的逆元 for (int i = 0; i < lim; ++i) { a[i] = 1ll * a[i] * inv_lim % P; } } }

逐段解析：

位逆序置换：通过交换，将数据排列成迭代计算需要的顺序。判断i < r[i]是为了避免重复交换。
迭代合并：最外层循环m从 2 开始，每次翻倍，模拟自底向上的合并过程。对于每一层，我们计算该层对应的主单位根gn。内层循环遍历每一个长度为m的块，并对块内前后两半进行蝶形运算。变量g从gn^0 = 1开始，每次乘以gn，依次得到gn^1, gn^2, ...。
逆变换处理：逆变换（INTT）有两个不同点：一是使用的单位根是原单位根的逆元（gn^{-1}），这通过qpow(gn, P-2)实现；二是最终结果需要乘以n^{-1}，即lim的逆元。

注意事项：这里有一个非常关键的细节，也是新手容易出错的地方——逆变换时单位根的处理。正向变换用gn = g^((P-1)/m)，逆向变换必须用其逆元gn_inv = g^(-(P-1)/m)。由于在模运算中，求逆元就是求幂的P-2次方，所以代码中直接用qpow(gn, P-2)得到。这个点如果搞反，结果会是错误的。

4. 实战应用：多项式乘法与卷积

NTT 最直接的应用就是计算两个多项式 A(x) 和 B(x) 的乘积 C(x) = A(x) * B(x)。设 A 的系数数组为a[]（次数为 n-1），B 的系数数组为b[]（次数为 m-1），那么 C 的次数为 n+m-2，共有 n+m-1 个系数。

卷积定理告诉我们，时域（系数）的卷积等于频域（点值）的乘积。因此，NTT 求解多项式乘法的步骤如下：

确定长度：找到大于等于n+m-1的最小的 2 的幂，作为变换长度lim。这是为了满足 Cooley-Tukey 算法对长度的要求。
初始化：将数组a和b的长度扩充到lim，多余的位置补 0。同时计算好位逆序数组r[]。
正向变换：对a和b分别进行ntt(..., lim, 1)，得到它们在“频域”的点值表示A'和B'。
点值相乘：在频域，卷积就是对应点值相乘：C'[i] = 1ll * A'[i] * B'[i] % P。
逆向变换：对C'进行ntt(..., lim, -1)，将其变换回“时域”，即得到卷积结果系数数组c。

下面是一个完整的卷积函数示例：

// 假设已经定义了 qpow, ntt 函数以及全局变量 r[] void convolve(int *a, int n, int *b, int m, int *c) { // 1. 确定变换长度 lim int lim = 1; while (lim < n + m - 1) lim <<= 1; // 2. 初始化位逆序数组 for (int i = 0; i < lim; ++i) { r[i] = (r[i >> 1] >> 1) | ((i & 1) ? (lim >> 1) : 0); } // 3. 准备扩展数组（这里为了清晰，使用了动态分配，实际比赛可能用全局数组） int *A = new int[lim](); // 初始化为0 int *B = new int[lim](); memcpy(A, a, n * sizeof(int)); memcpy(B, b, m * sizeof(int)); // 4. 正向NTT ntt(A, lim, 1); ntt(B, lim, 1); // 5. 点值相乘 for (int i = 0; i < lim; ++i) { A[i] = 1ll * A[i] * B[i] % P; } // 6. 逆向NTT ntt(A, lim, -1); // 7. 输出结果 memcpy(c, A, (n + m - 1) * sizeof(int)); delete[] A; delete[] B; }

性能与优化提示：

避免动态内存分配：在算法竞赛中，频繁的new/delete或vector构造开销很大。通常的做法是预先分配好足够大的全局数组（例如int A[MAXL], B[MAXL]），然后每次卷积时复用它们。
长度计算：while (lim < n + m - 1) lim <<= 1;这个循环可以用__builtin_clz（计算前导零）来优化，但可读性会下降。在不是极端卡常的情况下，这个循环足够了。
清零操作：扩展的部分（A[n..lim-1]）必须清零。使用()初始化或memset确保无误。

5. 关键细节与边界处理

实现一个健壮的 NTT，需要注意许多细节。这些细节往往决定了代码在边界情况下是否能正确运行。

5.1 变换长度的对齐

NTT 要求长度是 2 的幂。我们通过while (lim < need) lim <<= 1;来找到合适的lim。这里的need通常是n + m - 1。但有一个特殊情况：如果n或m为 0（空多项式），那么need可能为负数（当 n=0, m=0时，n+m-1 = -1）。因此，更安全的写法是：

int need = n + m - 1; if (need <= 0) { /* 处理空多项式情况 */ } int lim = 1; while (lim < need) lim <<= 1;

5.2 原地运算与临时变量

观察蝶形运算的代码：

int u = a[i + j]; int v = 1ll * a[i + j + k] * g % P; a[i + j] = (u + v) % P; a[i + j + k] = (u - v + P) % P;

这里必须使用临时变量u和v。如果直接写成：

a[i + j] = (a[i + j] + 1ll * a[i + j + k] * g) % P; // 错误！ a[i + j + k] = (a[i + j] - 1ll * a[i + j + k] * g + P) % P; // 此时 a[i+j] 已改变！

第二行计算时，a[i+j]已经被第一行修改了，导致结果错误。这是一个经典的“原地运算”陷阱。

5.3 逆变换的缩放因子

在逆变换的最后，需要乘以lim的逆元。这是因为正向变换可以看作乘以一个范德蒙德矩阵V，其元素是单位根的幂。而V^{-1} = (1/n) * V^H（在复数域是共轭转置，在模域是逆序和逆元）。所以逆变换后，每个元素都放大了n（即lim）倍，需要缩回去。

公式推导：设正向变换为Y = V * X，则X = V^{-1} * Y = (1/n) * V^H * Y。在模 p 下，1/n就是n^{-1} mod p。

务必记住：如果你忘记在逆变换时除以lim，得到的结果将是正确结果的lim倍（模 p 下）。这通常是调试时第一个要检查的点。

5.4 负数的处理

在模运算中，(u - v) % P可能得到负数。为了保证结果在[0, P-1]范围内，我们总是写成(u - v + P) % P。虽然在某些编译器或架构下，负数取模可能得到我们想要的结果，但为了可移植性和代码清晰，显式地+P是更好的习惯。

6. 模板代码与性能优化

经过上面的讨论，我们可以整合出一个竞赛中常用的、经过优化的 NTT 模板。这个模板包含了预处理逆元、预处理单位根等常见优化。

#include <bits/stdc++.h> using namespace std; const int MAXL = 1 << 21; // 支持的最大长度，2^21 const int P = 998244353, G = 3; // 模数，原根 int qpow(int x, int y) { int res = 1; for (; y; y >>= 1, x = 1ll * x * x % P) if (y & 1) res = 1ll * res * x % P; return res; } int r[MAXL], inv[MAXL + 5]; // 预处理逆元，可选，用于需要频繁逆变换的场景 void init_inv() { inv[1] = 1; for (int i = 2; i <= MAXL; ++i) inv[i] = 1ll * (P - P / i) * inv[P % i] % P; } // 预处理单位根，避免重复计算 qpow，显著优化性能 int omega[MAXL], omegaInv[MAXL]; void init_roots(int n) { // 计算主单位根及其逆元 int gn = qpow(G, (P - 1) / n); int gn_inv = qpow(gn, P - 2); omega[0] = omegaInv[0] = 1; for (int i = 1; i < n; ++i) { omega[i] = 1ll * omega[i - 1] * gn % P; omegaInv[i] = 1ll * omegaInv[i - 1] * gn_inv % P; } } void ntt(int *a, int lim, int *w) { // w 是预处理的单位根数组 for (int i = 0; i < lim; ++i) if (i < r[i]) swap(a[i], a[r[i]]); for (int m = 2; m <= lim; m <<= 1) { int k = m >> 1; int step = lim / m; // 单位根数组的步长 for (int i = 0; i < lim; i += m) { for (int j = 0; j < k; ++j) { int u = a[i + j]; int v = 1ll * a[i + j + k] * w[j * step] % P; // 查表 a[i + j] = (u + v) % P; a[i + j + k] = (u - v + P) % P; } } } } // 封装的卷积函数，使用全局数组避免动态分配 int A[MAXL], B[MAXL]; void convolve(int *a, int n, int *b, int m, int *c) { if (n == 0 || m == 0) { fill(c, c + n + m - 1, 0); return; } int need = n + m - 1; int lim = 1, l = 0; while (lim < need) lim <<= 1, l++; // 初始化位逆序 for (int i = 0; i < lim; ++i) r[i] = (r[i >> 1] >> 1) | ((i & 1) << (l - 1)); // 拷贝数据到全局数组并清零多余部分 copy(a, a + n, A); fill(A + n, A + lim, 0); copy(b, b + m, B); fill(B + m, B + lim, 0); // 预处理单位根（可以提前为常见长度预计算，这里每次计算） init_roots(lim); // 正向变换 ntt(A, lim, omega); ntt(B, lim, omega); // 点乘 for (int i = 0; i < lim; ++i) A[i] = 1ll * A[i] * B[i] % P; // 逆变换，使用逆单位根表 ntt(A, lim, omegaInv); // 缩放 int inv_lim = qpow(lim, P - 2); for (int i = 0; i < need; ++i) c[i] = 1ll * A[i] * inv_lim % P; }

优化点解析：

预处理单位根：在init_roots中，我们一次性计算出当前长度lim对应的所有单位根幂次，存储在数组omega和omegaInv中。在蝶形运算时直接查表w[j * step]，避免了在最内层循环中调用昂贵的qpow或进行乘法g = 1ll * g * gn % P。这是 NTT 最重要的优化之一，性能提升显著。
位逆序计算优化：r[i] = (r[i>>1]>>1) | ((i&1) << (l-1))是更标准的写法，其中l = log2(lim)。它与之前介绍的公式等价，但更简洁。
全局数组：使用全局数组A[MAXL], B[MAXL]避免了每次卷积时的内存分配，也便于清零操作。
边界处理：增加了对n或m为 0 情况的处理。

踩坑记录：我曾经在预处理单位根时犯过一个错误：step的计算应该是lim / m，而不是n / m或其他。step的意义是，对于长度为m的块，其使用的单位根在预处理数组中的间隔。因为预处理数组omega的长度是lim，存储了ω_lim^0, ω_lim^1, ..., ω_lim^(lim-1)。而当前层需要的单位根是ω_m^j = ω_lim^(j * (lim/m))，所以步长是lim/m。这个细节错了，整个变换结果就会乱掉。

7. 常见问题与调试技巧

即使理解了原理，实现 NTT 时也难免遇到各种 bug。下面是我在多年实践中总结的一些常见问题和调试方法。

7.1 结果不正确

这是最令人头疼的问题。可以按照以下步骤排查：

检查逆变换的缩放：这是最高频的错误。确保在ntt(..., -1)后，每个系数都乘以了lim的逆元。你可以用一个简单的例子测试：A(x)=1, B(x)=1，卷积结果应该是C(x)=1。如果得到的是lim，那就是忘记缩放了。
检查单位根：确保正向变换使用gn = g^((P-1)/m)，逆向变换使用其逆元。可以用小数据测试：对数组[1, 2, 3, 4]做一次 NTT 正变换，紧接着做一次逆变换，看是否能还原。
检查模运算和溢出：所有乘法都必须用1ll * a * b % P的形式，防止中间结果溢出int。加法减法后要立即取模，保证值在[0, P-1]范围内。
检查长度和清零：确保lim是 2 的幂且大于等于n+m-1。检查数组扩展部分是否被正确清零（特别是复用全局数组时）。
检查位逆序：可以打印出r[]数组检查。对于lim=8，正确的r[]应该是[0, 4, 2, 6, 1, 5, 3, 7]。

7.2 性能瓶颈

如果你的 NTT 运行太慢，可以考虑：

使用预处理单位根：如前所述，这是最有效的优化。
减少取模次数：在蝶形运算中，(u + v) % P和(u - v + P) % P可以优化。因为u和v都在[0, P-1]内，所以u+v在[0, 2P-2]内。可以这样写：
```
int t = u + v; a[i + j] = t >= P ? t - P : t; // 条件判断比取模快 t = u - v; a[i + j + k] = t < 0 ? t + P : t;
```
但要注意，这种优化可能会受编译器影响，有时现代编译器对%的优化已经很好。在竞赛极端优化时可以使用。
使用更大的模数：如果题目允许，使用constexpr模数，编译器可能进行更好的优化。
循环展开：手动或通过编译器指令展开最内层循环。例如，在k较小（如 4, 8）时，可以展开蝶形运算。但这会牺牲代码可读性。

7.3 小数据测试用例

这里提供几个用于验证 NTT 正确性的简单测试：

测试1：单位多项式

int a[] = {1}; int b[] = {1}; // 期望结果 c = {1}

测试2：常数多项式

int a[] = {5}; int b[] = {7}; // 期望结果 c = {35}

测试3：线性多项式

int a[] = {1, 2}; // 1 + 2x int b[] = {3, 4}; // 3 + 4x // 卷积： (1*3) + (1*4 + 2*3)x + (2*4)x^2 = {3, 10, 8}

测试4：随机数据对拍写一个朴素的 O(n²) 卷积函数，用随机生成的小规模数据（n, m < 100）与你的 NTT 结果对比。这是最可靠的验证方法。

7.4 内存与缓存优化

对于非常大的变换（例如lim > 2^20），缓存不命中会成为性能瓶颈。虽然迭代版 NTT 的访问模式已经比递归版更友好（主要是顺序访问），但仍可以进一步优化：

分块：将蝶形运算分块，使得每一块的数据能更好地留在缓存中。但这会大大增加代码复杂度，通常只在极其追求性能的库（如 FFTW）中实现。
使用std::vector<int>并预留空间：如果你使用vector，确保使用reserve()预留足够容量，避免多次重分配。
对齐内存：使用alignas或编译器扩展确保数组按缓存行对齐，但这属于高级优化，收益需要实测。

对于绝大多数竞赛和日常应用，前面给出的优化模板已经足够快。我的经验是，在 Codeforces 上，一个优化良好的 NTT 模板处理n=m=10^5的数据可以在 200ms 内完成，完全满足时限要求。

8. 扩展与变种

掌握了基础的 NTT 后，你可能会遇到一些变种或扩展需求。

8.1 任意模数NTT（MTT）

当模数不是 NTT 友好模数（如1e9+7）时，我们不能直接使用 NTT，因为1e9+7-1 = 10^9+6的因子 2 太少。这时需要使用任意模数 NTT。常见的方法有：

三模数NTT：选取三个 NTT 友好模数（如上述的 998244353, 1004535809, 469762049），分别做三次 NTT，然后用中国剩余定理（CRT）合并结果。由于三个模数的乘积大于可能的结果最大值，可以唯一确定结果。
拆系数FFT：将每个系数拆成a = x * M + y的形式（M 约等于 sqrt(P)），将多项式乘法转化为 4 次 FFT（复数域），最后合并。这种方法利用 FFT 的高精度浮点数计算，再取模还原。实现复杂，且要注意精度问题。

三模数 NTT 更常见，因为它完全在整数域运算，没有精度风险。但需要实现三个不同模数下的 NTT，以及 CRT 合并，代码量较大。