1. P4238 在考什么从“截断乘法”到“形式幂级数的除法”1.1 模 x^n 的乘法不是整多项式乘法P4238 这道题题面写得很学术又是“多项式乘法逆”又是“给定 (A(x)) 求 (B(x))”第一次接触的人往往会懵两个多项式相乘直接乘不就完了怎么还要逆元关键在于“模 (x^n)”这五个字。两个多项式正常相乘结果长度会变长但是题目只关心次数低于 n 的那些项即从常数项到 (x^{n-1}) 的系数次数大于等于 n 的项全部当作“截断扔掉”的高位噪声。这个操作有点像我平时做整数除法时只取不够一位除数的“余数”只不过这里的模对象变成了 (x^n)。于是“乘法逆元”的含义也随之改变了在模 (x^n) 意义下如果存在多项式 (B(x))使得[ A(x)\cdot B(x) \equiv 1 \pmod{x^n} ]成立那么我们就把 (B(x)) 叫作 (A(x)) 的乘法逆元。这里要注意它不是普通意义上两个多项式相除得到的“商”而是一种在截断环境里的特殊关系。简单说我们只需要保证乘积的低 n 项里常数项是 1其余项全是 0就算完成目标。刚学的时候我一度想不通为什么要定义这种看起来“不算完整”的乘法后来看多了多项式算子才明白很多算法需要的是“每算一步只往前推进一段精度”并不需要完整的无穷级数。比如后面要讲到的多项式对数、指数、开方全部依赖这种“截断到前 n 项”的运算环境。逆元就是第一个把这种环境用起来的经典模板。1.2 逆元是很多高级板子的地基除了题目本身是模板题之外多项式乘法逆最大的意义在于它被反复当作子过程使用多项式除法/取模需要用到除数的逆元多项式 (\ln) 要用到积分和逆元多项式 (\exp) 的牛顿迭代里又要反复调用逆元多项式开方也会借助逆元。可以这么理解只要能快速做完“一次 NTT 卷积”再配合“求逆”很多看起来复杂的高阶运算就能拼出来。这道题表面上只输出一组系数实际上是在给你搭一个工具后面同类板子几乎都从这里长出来。因此不要把 P4238 当成一道背代码题来刷它更值得做的是“推导一遍 对拍一遍 自己亲手把 NTT 调通”。调通之后你再看其他多项式模板会发现骨架都差不多区别只是迭代式长什么样。1.3 可解性的前提常数项不能是 0求逆有一个天然限制在模 (x^n) 意义下(A(x)) 的逆元存在当且仅当常数项 (a_0) 在模 998244353 意义下有逆元也就是 (a_0 \not\equiv 0)。道理很简单因为我们要让乘积的常数项等于 1而乘积的常数项只由 (a_0) 和 (b_0) 决定。若 (a_00)那么不管 (b_0) 取什么常数项都是 0永远凑不出 1。所以做题时第一步一定是先判断边界如果 (n1)答案就是 (a_0^{-1})。如果在实际工程里调用这个板子也应该在上层做好这个检查否则 NTT 函数进去了之后会发现怎么算都得不到正确结果反而浪费时间。这道模板题会保证逆元存在但我还是建议在自制板子里加一行防御性的特判。2. 为什么用 NTT998244353 这个模数的特殊结构2.1 模数是从哪儿“长”出来的题面给出的模数是 998244353也就是常说的 NTT 友好模数。很多人直接背代码却不理解为什么偏偏是这个数。看 (998244353-1998244352)对它做质因数分解你会发现[ 998244353 - 1 119 \times 2^{23} ]这意味着模数减一里面含有高达 23 次方的因子 2。而 NTT 要求变换长度必须是 2 的幂长度 (L) 最大不能超过模数里因子 2 的次数。998244353 支持到 (2^{23})约 838 万长度对于通常 (n \le 10^5) 的多项式题完全够用。此外998244353 的原根是 3。原根在 NTT 里扮演的角色相当于单位根在 DFT 里的角色。简单说我们能用 3 的幂次构造出一组长度恰好为 (2^k) 的“单位根”然后在这些根上做蝴蝶变换实现快速卷积。如果模数不满足“减一含大 2 的幂 有原根”这两个条件普通 FFT 就得用浮点数然后还要处理精度和取模问题那代码会痛苦得多。2.2 迭代版本的 NTT位逆置换 蝴蝶变换虽然 NTT 有很多写法但竞赛里最常见的还是迭代式。核心流程分三步先把系数按位逆序排列然后自底向上做长度 2、4、8……的蝴蝶变换最后如果是逆变换还要乘上长度的逆元。位逆置换可以不用递归用一个很经典的“按位爬台阶”写法for (int i 1, j 0; i n; i) { int bit n 1; for (; j bit; bit 1) j ^ bit; j ^ bit; if (i j) swap(a[i], a[j]); }第一次看这段代码会觉得像魔法其实它做的是把二进制下标倒过来例如长度为 8 时原下标 3 换成二进制 011反过来是 110也就是 6。这样换之后后续的分治过程就能用循环直接写不用递归也不用显式建一棵分治树。蝴蝶变换那步才真正消耗时间for (int len 2; len n; len 1) { int wlen ksm(3, (MOD - 1) / len); if (invert) wlen ksm(wlen, MOD - 2); for (int i 0; i n; i len) { int w 1; for (int j i; j i len / 2; j) { int u a[j]; int v 1LL * a[j len / 2] * w % MOD; a[j] (u v) % MOD; a[j len / 2] (u - v MOD) % MOD; w 1LL * w * wlen % MOD; } } }这组循环的思想是每一层把序列分成若干块每一块里做“一分为二”的合并。(uv) 是低频部分(u-v) 是高频部分用当前旋转因子 (w) 调整后就能把分治的两个小 DFT 合成一个大 DFT。逆变换最后乘上 (n^{-1}) 即可。2.3 一个可以直接抄的 NTT 核心我把整个 NTT 写成函数方便后面逆元复用#include bits/stdc.h using namespace std; typedef long long ll; const int MOD 998244353; const int G 3; int ksm(int a, int b) { int res 1; while (b) { if (b 1) res 1LL * res * a % MOD; a 1LL * a * a % MOD; b 1; } return res; } void ntt(vectorint a, bool invert) { int n (int)a.size(); for (int i 1, j 0; i n; i) { int bit n 1; for (; j bit; bit 1) j ^ bit; j ^ bit; if (i j) swap(a[i], a[j]); } for (int len 2; len n; len 1) { int wlen ksm(G, (MOD - 1) / len); if (invert) wlen ksm(wlen, MOD - 2); for (int i 0; i n; i len) { int w 1; for (int j i; j i len / 2; j) { int u a[j]; int v 1LL * a[j len / 2] * w % MOD; a[j] (u v) % MOD; a[j len / 2] (u - v MOD) % MOD; w 1LL * w * wlen % MOD; } } } if (invert) { int invn ksm(n, MOD - 2); for (int x : a) x 1LL * x * invn % MOD; } }这段代码适合先理解再动手改。真正比赛时限紧的时候把vector换成全局数组把ksm换成预处理的单位根表常数会明显下降这部分我后面单独说。3. 牛顿迭代为什么每次只“翻倍”而不是一步到位3.1 朴素做法和倍增想法的差距如果你第一次看到这题直接想到的朴素做法可能是解一个 (n) 元一次方程组。因为已知[ \sum_{j0}^{i} a_j b_{i-j} \begin{cases} 1 i0 \ 0 1 \le i n \end{cases} ]所以从 (b_0) 开始每一项都被前面已经算出的 (b_j) 决定。这样确实能 (O(n^2)) 算出答案但它不是我们要的模板算法。经典的 (O(n\log n)) 思路来源于“倍增”也叫二进制迭代或牛顿迭代。核心想法很简单如果我已经知道 (B_m(x)) 满足模 (x^m) 下的逆元关系那就通过一个固定的式子推出新的 (B_{2m}(x))让精度从 m 项直接翻倍到 2m 项。反复操作最终达到 n 项。为什么非要从 1 项、2 项、4 项这样长起来因为 NTT 的复杂度是 (O(L\log L))如果每一步只算当前需要的长度总时间会是一个等比数列求和最后收敛到 (O(n\log n))。如果一开始就硬算整个 n 项复杂度至少是 (O(n^2))那还不如别用 NTT。3.2 迭代公式的两种推导假设已有 (B_m) 满足[ A \cdot B_m \equiv 1 \pmod{x^m} ]我们要找到 (B_{2m})让[ A \cdot B_{2m} \equiv 1 \pmod{x^{2m}} ]最直观的代数推导是这样令误差[ E 1 - A \cdot B_m ]因为 (B_m) 已经是模 (x^m) 下的逆元所以 (E) 必定被 (x^m) 整除也就是说它的前 m 项全是 0。接下来构造[ B_{2m} B_m , (2 - A \cdot B_m) ]把它代进去[ \begin{aligned} A \cdot B_{2m} A B_m (2 - A B_m) \ 2 A B_m - (A B_m)^2 \end{aligned} ]已知 (A B_m 1 - E)所以[ A \cdot B_{2m} 2(1-E) - (1-E)^2 1 - E^2 ]由于 (E) 被 (x^m) 整除(E^2) 一定被 (x^{2m}) 整除。换句话说(A \cdot B_{2m}) 和 1 的差是 (x^{2m}) 的倍数这正是模 (x^{2m}) 意义下的逆元关系。另外一条推法是从牛顿迭代的角度看。函数 (F(B)AB-1)牛顿法求一个“根” (B) 让 (F(B)0)。直接套牛顿迭代[ B_{new} B - \frac{F(B)}{F(B)} ]但这里除法代价不小。绕一下用当前近似逆元 (B) 去代替 (F(B)^{-1})得到[ B_{new} B - (AB - 1) B 2B - A B^2 B(2 - AB) ]正好就是同一个式子。因此这个模板题背后的数学结构本质上就是牛顿法在多项式环上的应用。3.3 拿 (1x) 验证一遍公式理论推导完我习惯先手算一个简单例子压压惊。设[ A(x) 1 x ]它在模 (x^4) 下的逆元应该是 (1 - x x^2 - x^3)因为[ (1x)(1-xx^2-x^3) 1 - x^4 ]现在看迭代过程对不对。初始只有一项(B_1 1)。第一轮从长度 1 翻到 2[ A B_1 (1x)\cdot 1 1x ]所以[ 2 - A B_1 2 - (1x) 1 - x ]于是[ B_2 B_1(2 - AB_1) 1-x ]这就得到了正确的前两项 (1-x)。第二轮从长度 2 翻到 4。计算[ A B_2 (1x)(1-x)1-x^2 ]注意这里的 (1-x^2) 是完整乘积但我在计算 (2-AB_2) 时只需保留前 4 项也就是 (2-(1-x^2)1x^2)。接着[ B_4 B_2 (2-AB_2) (1-x)(1x^2) 1-xx^2-x^3 ]正好和直接展开完全一样。这个例子虽然简单却能很清楚地看到“前 m 项保持正确、多出的 m 项被修正”的过程理解了它代码里的循环结构就该浮出水面了。4. 手写代码给每个 NTT 卷积安排好位置4.1 每一轮到底要算几个卷积有了迭代公式剩下的问题怎么用 NTT 实现。每一轮从长度 len 翻到 curcur 一般是 2len最后一轮可能是 n要做的事情其实是两个卷积第一步计算 (A_{\text{cur}} \cdot B_{\text{len}})。这里的 (A_{\text{cur}}) 表示取原多项式前 cur 个系数。只算这个卷积的前 cur 项就够了但 NTT 的卷积长度不能只开 cur原因我下一节详细讲。第二步构造 (2 - A B) 的系数也就是把第一步卷积结果取负常数项加 2。再把这个新多项式与 (B_{\text{len}}) 做卷积。第二次卷积也只需要前 cur 项。两个卷积中间有一个很省时间的细节因为第二次卷积还是和同一个 (B_{\text{len}}) 相乘所以 (B_{\text{len}}) 的 NTT 结果可以直接复用。第一次算 (A \cdot B_{\text{len}}) 时我们已经对 (B_{\text{len}}) 做过一次正变换第二次接着用那份点值就行不用再变换一遍。4.2 非递归的倍增实现我的习惯是写非递归版本这样不容易在递归边界上出错也方便临时换数组。代码框架如下vectorint poly_inv(const vectorint a, int n) { // 初始长度为 1逆元的常数项就是 a0 的逆元 vectorint b(1, ksm(a[0], MOD - 2)); for (int len 1; len n; len 1) { int cur min(len 1, n); // 本轮要推进到的长度 int size 1; // NTT 长度必须能装下完整线性卷积 while (size cur len - 1) size 1; vectorint f(size, 0), g(size, 0), h(size, 0); for (int i 0; i cur i (int)a.size(); i) f[i] a[i]; for (int i 0; i len; i) g[i] b[i]; ntt(f, false); ntt(g, false); // h A[cur] * B[len] for (int i 0; i size; i) h[i] 1LL * f[i] * g[i] % MOD; ntt(h, true); // 只保留前 cur 项然后取负常数项加 2 h.resize(cur); for (int i 0; i cur; i) h[i] (MOD - h[i]) % MOD; h[0] (h[0] 2) % MOD; // h (2 - A*B) * B[len] h.resize(size, 0); ntt(h, false); for (int i 0; i size; i) h[i] 1LL * h[i] * g[i] % MOD; ntt(h, true); // 新逆元就是 h 的前 cur 项 b.assign(h.begin(), h.begin() cur); } b.resize(n); return b; }有个容易困惑的细节为什么第二次卷积后要直接b.assign(...)而不是只把新增的半段接上去实际上按迭代公式前 len 项的确应该和旧 (B_{\text{len}}) 完全相同理论上是不会变的。但直接整体赋值更省心也避免了“旧数组后半段还是脏数据”这类隐患。只要代码里没有写错整体赋值不会带来任何额外开销。4.3 完整可提交版本串起来之后主函数很简单int main() { int n; scanf(%d, n); vectorint a(n); for (int i 0; i n; i) scanf(%d, a[i]); vectorint ans poly_inv(a, n); for (int i 0; i n; i) printf(%d , ans[i]); printf(\n); return 0; }输入一个 (n)然后是 (a_0, a_1, \dots, a_{n-1})输出就是逆元数组。这段代码在正确性上没有难点主要是每一步都要记得取模和 resize。我第一次写的时候问题恰恰出在“h.resize(cur)”之后忘记再补回 size 个 0第二次ntt(h, false)直接对长度为 cur 的数组做变换长度对不上点值域NTT 结果自然全错。5. 那些必须踩一次才会避开的坑5.1 NTT 长度不足引发的混叠看似能算实则污染这是本题最阴间的坑没有之一。我的理解是这样的一次正常的多项式乘法如果两个多项式分别有 (L_1) 项和 (L_2) 项那么完整乘积最多有 (L_1L_2-1) 项。而 NTT 的变换长度 (size) 必须不小于这个数得到的才是完整线性卷积。有些题解为了省时间会把 NTT 长度只开到当前目标长度 cur得到一个“循环卷积”。循环卷积会有高位项卷回低位[ z_k c_k c_{ksize} c_{k2size} \cdots ]本来不需要的高位项会混到低位里污染正确答案。尤其做多项式逆元的时候每一轮只看前 cur 项高位被截断的部分恰好是逆元的高阶项混回来之后整段数据全烂。所以我在代码里坚持让size至少是cur len - 1的 2 的幂。虽然这样会偶尔多开一些长度比如 (cur4, len2) 时curlen-15所以 size 开到 8看起来有点浪费但换来的是正确性值。5.2 单位根反复求幂会卡常数上面给的ntt里每一层蝴蝶变换都会调用ksm求一次wlen。以 (n10^5) 级别的题目来说代码能过但如果将来做更长的多项式、更极限的数据每次ksm的开销会非常明显。常用的优化方法有两种一是预处理单位根表把每个长度对应的 (wlen) 先算好存起来。二是更彻底的做法预处理所有分治层需要的旋转因子甚至把w的递推也换成查表。但预处理代码写起来复杂得多建议先跑通朴素 NTT再根据时限决定要不要优化。还有一个更小的优化逆变换的 (wlen) 可以直接用原根高次幂的逆或者用ksm(wlen, MOD - 2)求。需要注意的是不能在逆变换后少乘n的逆元。我见过不少人 NTT 正变换没问题逆变换忘了最后统一乘invn结果所有系数都比真实值大了若干倍。5.3 边界条件n1、常数项、以及 resize 的时机边界条件里有三个点特别容易炸第一(n1) 时循环根本不会进入直接返回长度为 1 的数组[a0^{-1}]这是对的。但要保证a[0]不为 0否则ksm返回 0题目却没答案。第二h.resize(cur)之后如果你马上就做取负操作没问题但如果接着又h.resize(size, 0)注意这个调用会把新增位置补 0但是不会清除前 cur 个位置所以不要以为“resize 会清空整个数组”。第三读入时如果题目给的 (n) 是 1而数组长度是 1poly_inv里for (int i 0; i cur i (int)a.size(); i)这种写法要保证不越界。很多人默认传进来的数组长度就是 n但实际调用时可能会传一个更长的空余数组所以我还是写上了取地址的长度限制。另外b.assign(h.begin(), h.begin() cur)这个操作要求h至少有 cur 个元素。第二次 NTT 之后h保留了整个 size 长度的数组所以不会出问题。如果为了省内存把h再缩到 cur那反而会丢数据别这么做。6. 怎么验证你的逆元真的对对拍才是硬道理6.1 先跑特殊数据直接用 (1x)写完板子第一件事不是拿样例测而是找这种“一眼能看出答案”的数据。比如输入4 1 1 0 0因为 (1x) 在模 (x^4) 下的逆元是 (1-xx^2-x^3)输出应该是1 998244352 1 998244352这里的 998244352 就是 (-1) 在模 998244353 法则下的表示。如果这组数据过了说明整条迭代链路的框架基本正确。6.2 用朴素 (O(n^2)) 对拍光过一组不够我一般会写一个暴力版本进行随机对拍。暴力版不需要 NTT直接用递推关系[ b_0 a_0^{-1},\quad b_i -a_0^{-1} \sum_{j1}^{i} a_j b_{i-j} ]写成代码就是vectorint brute_inv(const vectorint a, int n) { vectorint b(n); b[0] ksm(a[0], MOD - 2); for (int i 1; i n; i) { int sum 0; for (int j 1; j i; j) { sum (sum 1LL * a[j] * b[i - j]) % MOD; } b[i] 1LL * (MOD - sum) * b[0] % MOD; } return b; }然后随机生成长度不等的数组长度覆盖 1、2、3、4、5、8、9、16、17、100 这类不规整数字。比较poly_inv和brute_inv的每一项。随机对拍能很快暴露长度不足、resize 错位、取模非法等问题。这里有个经验对拍时别忘了测试“n 不是 2 的幂”的情况。因为倍增循环里最后一轮的cur min(len 1, n)不等于 2lensize的计算逻辑、b.assign的长度都会走到特殊分支。很多实现能过整齐的 2 的幂数据一到 n100000 就出问题就是因为没有覆盖这种“半截”长度。6.3 时间和空间观察如果题目数据达到了 (n10^5) 级别这段 vector 版本通常能跑但不会特别快。原因在于每轮迭代要开 3 个 vector还要做 5 次 NTT 变换常数偏大。如果你卡时间可以改成全部使用全局数组预先分配好最大长度NTT 里把ksm(G, ...)改成预处理表第二次卷积时继续复用g的点值结果这些都能省下不少时间。空间上主要就是每轮三个 size 长度的数组。最坏那一轮 size 大约是2n向上取整到 2 的幂所以数组长度约 (4n) 的量级也是正常的。如果 (n10^5)开 (4\times 10^5) 长度的int数组完全不用担心内存。我做这题最大的感悟是多项式模板题不是“把公式背下来”而是“把每一轮为什么这么开长度、为什么能复用 NTT 点值理解清楚”。一旦你把 P4238 的迭代式吃透后面学到多项式 (\exp)、(\ln)、开方的时候框架基本可以直接平移过去只是迭代公式换了一下。