1. 项目概述为什么我们需要NTT在算法竞赛和高性能计算领域多项式乘法是一个绕不开的核心问题。无论是大整数乘法、字符串匹配还是生成函数相关的计数问题最终都归结为两个多项式系数的卷积运算。最直观的算法是 O(n²) 的双层循环这在数据规模稍大时比如 n10⁵就完全不可行了。这时快速傅里叶变换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 为偶数时求和引理Σ_{j0}^{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 1gn^(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 3469762049 7 * 2^26 1,g 3167772161 5 * 2^25 1,g 3当结果可能超过单个模数的范围时可以使用任意模数 NTTMTT通常通过使用两个或多个上述模数分别进行 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 和 ik 进行操作 (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 位逆序置换迭代算法需要先将输入序列按照“位逆序”重新排列。什么是位逆序对于一个长度为 n2 的幂的序列下标 i 的二进制表示反转后就得到了它应该放置的新位置r[i]。 例如n8 时i1 (001)的逆序是4 (100)。i3 (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) { // opt1 为正变换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-1B 的系数数组为b[]次数为 m-1那么 C 的次数为 nm-2共有 nm-1 个系数。卷积定理告诉我们时域系数的卷积等于频域点值的乘积。因此NTT 求解多项式乘法的步骤如下确定长度找到大于等于nm-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可能为负数当 n0, m0时nm-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[ij] 已改变第二行计算时a[ij]已经被第一行修改了导致结果错误。这是一个经典的“原地运算”陷阱。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[i1]1) | ((i1) (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 的幂且大于等于nm-1。检查数组扩展部分是否被正确清零特别是复用全局数组时。检查位逆序可以打印出r[]数组检查。对于lim8正确的r[]应该是[0, 4, 2, 6, 1, 5, 3, 7]。7.2 性能瓶颈如果你的 NTT 运行太慢可以考虑使用预处理单位根如前所述这是最有效的优化。减少取模次数在蝶形运算中(u v) % P和(u - v P) % P可以优化。因为u和v都在[0, P-1]内所以uv在[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::vectorint并预留空间如果你使用vector确保使用reserve()预留足够容量避免多次重分配。对齐内存使用alignas或编译器扩展确保数组按缓存行对齐但这属于高级优化收益需要实测。对于绝大多数竞赛和日常应用前面给出的优化模板已经足够快。我的经验是在 Codeforces 上一个优化良好的 NTT 模板处理nm10^5的数据可以在 200ms 内完成完全满足时限要求。8. 扩展与变种掌握了基础的 NTT 后你可能会遇到一些变种或扩展需求。8.1 任意模数NTTMTT当模数不是 NTT 友好模数如1e97时我们不能直接使用 NTT因为1e97-1 10^96的因子 2 太少。这时需要使用任意模数 NTT。常见的方法有三模数NTT选取三个 NTT 友好模数如上述的 998244353, 1004535809, 469762049分别做三次 NTT然后用中国剩余定理CRT合并结果。由于三个模数的乘积大于可能的结果最大值可以唯一确定结果。拆系数FFT将每个系数拆成a x * M y的形式M 约等于 sqrt(P)将多项式乘法转化为 4 次 FFT复数域最后合并。这种方法利用 FFT 的高精度浮点数计算再取模还原。实现复杂且要注意精度问题。三模数 NTT 更常见因为它完全在整数域运算没有精度风险。但需要实现三个不同模数下的 NTT以及 CRT 合并代码量较大。8.2 快速沃尔什变换FWTNTT 解决的是卷积c_k Σ_{ijk} a_i * b_j。还有一类问题是位运算卷积例如与卷积c_k Σ_{ijk} a_i * b_j或卷积c_k Σ_{i|jk} a_i * b_j异或卷积c_k Σ_{i^jk} a_i * b_j解决这类问题需要快速沃尔什变换FWT。有趣的是FWT 的算法结构与 FFT/NTT 非常相似也是基于分治和线性变换只是变换矩阵不同。如果你深刻理解了 NTT 的蝶形运算结构学习 FWT 会容易得多。8.3 多项式求逆、exp、ln等NTT 是多项式高级操作的基础。有了高效的卷积我们就可以在 O(n log n) 时间内实现多项式求逆给定多项式 A(x)求 B(x) 使得 A(x)B(x) ≡ 1 (mod x^n)。使用牛顿迭代法。多项式开根求 B(x) 使得 B(x)² ≡ A(x) (mod x^n)。多项式指数函数exp求 B(x) 使得 B(x) ≡ exp(A(x)) (mod x^n)。多项式对数函数ln求 B(x) 使得 B(x) ≡ ln(A(x)) (mod x^n)。这些操作构成了现代多项式算法的工具箱在生成函数、组合计数等问题中威力巨大。它们都依赖于 NTT 进行快速多项式乘法。9. 总结与个人体会实现一个正确且高效的 NTT就像搭建一个精密的机械钟表每一个齿轮步骤都必须严丝合缝。从理解原根与单位根的对应关系到实现迭代蝶形算法再到处理逆变换和边界条件每一步都需要清晰的逻辑和对细节的把握。我个人在最初实现 NTT 时曾花了整整一个晚上调试一个因为忘记逆变换缩放而导致的问题。从那以后我养成了一个习惯任何新的 NTT 实现首先用常数多项式和线性多项式测试。这两个简单测试能排除掉大部分低级错误。另一个深刻的体会是优化代码的前提是代码正确。不要一开始就追求极致的性能比如用位运算代替除法展开循环。先写出一个清晰、正确的版本确保逻辑无误。然后再逐步加入预处理单位根、减少取模等优化并且每加一个优化都要重新测试。这样当出现错误时你才能快速定位是哪个“优化”引入了问题。最后NTT 不仅仅是一个算法模板它更是一种思想——在看似不同的数学领域复数和有限域之间建立桥梁将一个问题转化为另一个更容易计算的问题。这种思想在算法设计中无处不在。当你下次遇到一个复杂的问题时不妨想想有没有一种“变换”能让它变得简单