多项式全家桶核心原理：牛顿迭代法统一求逆、开根、ln与exp

1. 项目概述：从“黑盒”到“白盒”的多项式运算工具箱

在算法竞赛和理论计算机科学领域，多项式运算早已不是新鲜话题。从基础的加减乘，到稍显复杂的求逆、开根，再到更高级的对数（ln）和指数（exp）运算，这些操作构成了一个被称为“多项式全家桶”的核心工具集。很多选手在初次接触时，往往直接套用模板，将其视为一个功能强大的“黑盒”——输入系数数组，调用函数，得到结果。然而，当题目条件发生变化，或者需要优化常数、处理边界情况时，仅仅会“用”是远远不够的。理解其背后的数学原理和实现细节，将这个“黑盒”彻底“白盒化”，才是从使用者迈向创造者的关键一步。

这篇文章，就是一次彻底的“白盒化”旅程。我们不满足于仅仅给出代码模板，而是要深入拆解多项式求逆、开根、ln、exp这四大核心操作的每一个步骤。我会结合自己多年打比赛和出题的经验，详细解释牛顿迭代法是如何在这些运算中统一应用的，FFT/NTT（快速傅里叶变换/数论变换）在其中扮演了什么角色，以及那些模板代码里看似神秘的参数和边界处理究竟为何如此设计。无论你是正在备赛的选手，还是对算法实现细节有浓厚兴趣的开发者，相信这篇超过5000字的深度解析，都能让你对“多项式全家桶”有一个全新的、透彻的认识。

2. 核心数学原理与统一框架：牛顿迭代法

在深入每个具体操作之前，我们必须先建立一个统一的视角。多项式求逆、开根、ln、exp，这些看似不同的运算，在算法实现层面，其实共享着一个强大的核心工具：牛顿迭代法。理解这一点，是理解整个“全家桶”的钥匙。

2.1 牛顿迭代法在多项式语境下的重塑

我们都知道牛顿迭代法用于求实数方程的根：从一个初始近似值x0开始，通过公式x_{n+1} = x_n - f(x_n)/f'(x_n)不断逼近真实根。在多项式运算中，我们将其进行了一次巧妙的“移植”。我们不再求一个实数的根，而是求一个“多项式函数”的“零点多项式”。

具体来说，我们把目标运算转化为求解一个关于多项式G(x)的方程F(G(x)) = 0。这里F是一个将多项式映射到多项式的函数。例如：

• 求逆：设A(x)是已知多项式，我们想求B(x)使得A(x) * B(x) ≡ 1 (mod x^n)。这可以转化为方程F(B) = 1/B - A ≡ 0，或者更常用的F(B) = A - 1/B ≡ 0。
• 开根：设A(x)是已知多项式，我们想求B(x)使得B(x)^2 ≡ A(x) (mod x^n)。这转化为方程F(B) = B^2 - A ≡ 0。

牛顿迭代公式相应地变为：G_{new}(x) ≡ G_{old}(x) - F(G_{old}(x)) / F'(G_{old}(x)) (mod x^m)。这里的导数F'是对G求导，mod x^m表示我们只关心前m项系数。关键在于，通过每次迭代，解的有效位数（即正确的项数）会翻倍。如果我们从mod x^1（即仅常数项正确）开始，那么迭代一次得到mod x^2正确的解，再迭代一次得到mod x^4，以此类推。这种“倍增”思想是算法效率的基石。

注意：这里的“导数”是形式导数，完全按照多项式求导法则进行，与实数微积分中的导数意义不同，但运算法则一致。例如，(x^n)' = n*x^{n-1}。

2.2 迭代的起点与收敛性分析

牛顿迭代需要一个初始值。对于多项式运算，这个初始值通常是常数项。

• 求逆：要求A(x)的常数项a0在模意义下可逆（即存在逆元）。我们取B_0(x) ≡ inv(a0) (mod x)，这里inv(a0)是a0的模逆元。
• 开根：要求A(x)的常数项a0是模意义下的二次剩余（即存在b0使得b0^2 ≡ a0）。我们取B_0(x) ≡ sqrt(a0) (mod x)，这里sqrt(a0)是a0的模平方根之一。
• ln和exp：ln要求多项式常数项为1，exp要求多项式常数项为0。它们的迭代通常依赖求逆操作，因此初始值也隐含在其中。

如果初始条件不满足，相应的运算在模意义下无解。这是实现时必须检查的第一步。

3. 核心操作一：多项式求逆详解

多项式求逆是全家桶中最基础也是最重要的操作，它是实现开根、ln、exp的基石。其定义是：对于多项式A(x)，求多项式B(x)，使得A(x) * B(x) ≡ 1 (mod x^n)。这里n是我们需要的项数。

3.1 牛顿迭代推导与实现步骤

我们从方程F(B) = A - 1/B ≡ 0出发。对F(B)关于B求导：F'(B) = 1/B^2。 代入牛顿迭代公式：B_{new} ≡ B_{old} - (A - 1/B_{old}) / (1/B_{old}^2) (mod x^{2m})≡ B_{old} - (A - 1/B_{old}) * B_{old}^2 (mod x^{2m})≡ B_{old} * (2 - A * B_{old}) (mod x^{2m})

这就是多项式求逆的核心迭代式：B ≡ B * (2 - A * B) (mod x^{2m})。

实现步骤（递归/倍增法）：

1. 边界条件：当n=1时，直接返回常数项a0的逆元。若a0无逆元，则整个运算无解。
2. 递归求解：设当前需要求解mod x^n下的逆。我们先递归求解mod x^{ceil(n/2)}下的逆，记为B0(x)。ceil表示向上取整，这是为了保证倍增性质。
3. 迭代升级：
   • 计算T(x) = A(x) * B0(x) (mod x^n)。注意，这里我们只需要前n项。
   • 计算R(x) = 2 - T(x) (mod x^n)。因为T(x)的前ceil(n/2)项应该是1（理论上），所以R(x)的前ceil(n/2)项为1，这保证了迭代的有效性。
   • 计算B(x) = B0(x) * R(x) (mod x^n)。根据公式，这就是新的、在mod x^n下更精确的逆。
4. 返回结果：返回B(x)。

在这个过程中，核心的乘法运算（A*B0,B0*R）都需要使用NTT来加速。这也是为什么多项式全家桶通常要求模数满足 NTT 条件（如 998244353，其原根为 3）。

3.2 关键细节与常数优化

• 长度与清零：进行 NTT 乘法时，长度必须扩展到大于等于2*n的最近 2 的幂。计算完成后，务必手动将n之后的系数清零（mod x^n的含义），否则在后续运算中会引入错误的高次项。
• 临时数组复用：为了减少内存分配和拷贝开销，通常会预分配几个大的临时数组，在函数间传递并复用。计算T(x)和R(x)时可以共用数组。
• 迭代与递归的选择：上述描述是递归形式，易于理解。在实际的高性能模板中，往往采用等价的循环倍增实现，以避免递归的函数调用开销。即从mod x^1开始，不断进行B = B * (2 - A * B) (mod x^{2m})，每次m翻倍，直到m >= n。

// 伪代码示意：多项式求逆 (循环倍增版本) void poly_inv(int *a, int *b, int n) { // b 是输出数组，初始为空 static int tmp[MAXN]; b[0] = qpow(a[0], MOD-2); // 初始值: mod x^1 for (int m = 1; m < n; m <<= 1) { // m 是当前已知正确的长度 int lim = m << 2; // 计算长度，为4*m以保证NTT精度 // 将 a 的前 2*m 项拷贝到 tmp_a // 将 b 的前 m 项拷贝到 tmp_b (实际上b就是当前的逆) ntt(tmp_a, lim, 1); ntt(tmp_b, lim, 1); for (int i=0; i<lim; ++i) tmp[i] = (2 - (ll)tmp_a[i]*tmp_b[i]%MOD + MOD) % MOD * tmp_b[i] % MOD; ntt(tmp, lim, -1); // 将 tmp 的前 2*m 项拷贝到 b，作为新的逆 for (int i=2*m; i<lim; ++i) b[i] = 0; // 清空高位 } for (int i=n; i<lim; ++i) b[i] = 0; // 最终只保留 n 项 }

4. 核心操作二：多项式开根解析

多项式开根，即求B(x)使得B(x)^2 ≡ A(x) (mod x^n)。有了求逆的基础，开根的实现就清晰多了。

4.1 基于牛顿迭代的推导

方程是F(B) = B^2 - A ≡ 0。求导得F'(B) = 2B。 代入牛顿迭代公式：B_{new} ≡ B_{old} - (B_{old}^2 - A) / (2B_{old}) (mod x^{2m})≡ (B_{old} + A / B_{old}) / 2 (mod x^{2m})≡ (B_{old} + A * inv(B_{old})) / 2 (mod x^{2m})

看，这里出现了多项式求逆inv(B_{old})。所以，开根的实现依赖于求逆。迭代式简化为：B ≡ (B + A * B^{-1}) / 2。

实现步骤：

1. 边界与初始值：当n=1时，返回常数项a0的模平方根（需预处理或使用 Cipolla 算法求解）。同样要求a0是二次剩余。
2. 递归/倍增求解：类似求逆，先递归求出mod x^{ceil(n/2)}下的根B0(x)。
3. 迭代升级：
   • 调用poly_inv计算B0(x)在mod x^n下的逆元I0(x)。
   • 计算T(x) = A(x) * I0(x) (mod x^n)。
   • 计算B(x) = (B0(x) + T(x)) * inv2 (mod x^n)。其中inv2是 2 的模逆元（在模意义下除以2等于乘以2的逆元）。
4. 返回结果：返回B(x)。

4.2 实现要点与边界处理

• 依赖求逆：这是开根运算的核心特点。在代码组织上，开根函数内部会调用求逆函数。
• 常数项处理：如果常数项a0不是二次剩余，在模意义下无法开根。通常题目会保证有解，但自己写代码时要留心。
• 除以2的处理：在模MOD下，除以2必须转换为乘以 (MOD+1)/2（即2的逆元）。这是模运算的基本要求，绝对不能直接做整数除法。
• 长度管理：在计算A * I0时，A只需要前n项，I0是B0的逆，长度为n。NTT 长度需要妥善设置。

// 伪代码示意：多项式开根 void poly_sqrt(int *a, int *b, int n) { static int tmp_inv[MAXN], tmp_a[MAXN]; if (n == 1) { b[0] = sqrt_mod(a[0]); return; } // sqrt_mod 求模平方根 poly_sqrt(a, b, (n+1)/2); // 递归求前 ceil(n/2) 项 poly_inv(b, tmp_inv, n); // 求当前 b 的逆 int lim = get_lim(n*2); // 准备 a 的前 n 项到 tmp_a ntt(tmp_a, lim, 1); ntt(tmp_inv, lim, 1); for (int i=0; i<lim; ++i) tmp_a[i] = (ll)tmp_a[i] * tmp_inv[i] % MOD; ntt(tmp_a, lim, -1); int inv2 = (MOD+1)/2; for (int i=0; i<n; ++i) b[i] = (ll)(b[i] + tmp_a[i]) * inv2 % MOD; for (int i=n; i<lim; ++i) b[i] = 0; }

5. 核心操作三：多项式对数函数（ln）

多项式ln的定义需要借助微积分。对于常数项为1的多项式A(x)，定义ln(A(x))为其形式幂级数展开。在实际计算中，我们使用导数积分法。

5.1 利用导数与积分的计算原理

公式：ln(A(x)) = ∫ (A'(x) / A(x)) dx这里A'(x)是A(x)的导数，∫是积分，/是多项式除法（即乘以逆元）。

原理简述：对ln(A(x))两边求导，得到(ln(A(x)))' = A'(x) / A(x)。然后两边积分，就得到上面的公式。注意，因为常数项为1，所以积分后的常数项为0，符合ln1=0。

计算步骤：

1. 检查条件：确保A(x)的常数项a0 ≡ 1 (mod MOD)。否则ln在形式幂级数意义上无良好定义。
2. 求导：计算A(x)的导数A'(x)。这是一个O(n)的操作：A'[i] = A[i+1] * (i+1) % MOD。
3. 求逆：计算A(x)的乘法逆元A_inv(x) (mod x^n)。
4. 卷积：计算P(x) = A'(x) * A_inv(x) (mod x^{n-1})。因为A'的次数是n-2，A_inv我们取前n-1项，卷积后我们只需要前n-1项。
5. 积分：对P(x)进行积分，得到结果B(x)。积分公式：B[i] = P[i-1] * inv(i) % MOD，其中inv(i)是i的模逆元，且B[0] = 0（因为ln1=0）。

可以看到，多项式ln的核心是求逆和卷积。

5.2 实现细节与预处理优化

• 逆元预处理：积分时需要用到1, 2, ..., n-1的模逆元。可以提前用线性方法inv[i] = MOD - MOD/i * inv[MOD%i] % MOD预处理出来，避免在函数内重复计算。
• 长度匹配：求逆时我们要求mod x^n的逆，用于后续与A'（长度为n-1）的卷积。卷积结果我们只取前n-1项用于积分。
• 常数项处理：输入必须保证A[0]==1。输出结果的B[0]固定为0。

// 伪代码示意：多项式 ln void poly_ln(int *a, int *b, int n) { static int tmp_a[MAXN], tmp_inv[MAXN]; // 步骤1: 求导 for (int i=0; i<n-1; ++i) tmp_a[i] = (ll)a[i+1] * (i+1) % MOD; // 步骤2: 求逆 poly_inv(a, tmp_inv, n); // 求 a 的逆，长度为 n // 步骤3: 卷积 int lim = get_lim(2*n); ntt(tmp_a, lim, 1); ntt(tmp_inv, lim, 1); for (int i=0; i<lim; ++i) tmp_a[i] = (ll)tmp_a[i] * tmp_inv[i] % MOD; ntt(tmp_a, lim, -1); // 此时 tmp_a 的前 n-1 项是 A'/A // 步骤4: 积分 b[0] = 0; for (int i=1; i<n; ++i) b[i] = (ll)tmp_a[i-1] * inv[i] % MOD; // inv[i] 已预处理 for (int i=n; i<lim; ++i) b[i] = 0; }

6. 核心操作四：多项式指数函数（exp）

多项式exp是ln的逆运算。求B(x)使得ln(B(x)) ≡ A(x) (mod x^n)，且B(x)常数项为1。这是全家桶中实现最复杂的一环，通常使用牛顿迭代法。

6.1 牛顿迭代法的再次应用

设F(B) = ln(B) - A ≡ 0。求导得F'(B) = 1/B。 代入牛顿迭代公式：B_{new} ≡ B_{old} - (ln(B_{old}) - A) / (1/B_{old}) (mod x^{2m})≡ B_{old} * (1 - ln(B_{old}) + A) (mod x^{2m})

化简后得到核心迭代式：B ≡ B * (1 - ln(B) + A)。注意，这里ln(B)是多项式对数运算。

实现步骤（倍增法）：

1. 边界与初始值：当n=1时，B(x) ≡ 1 (mod x)，因为exp(0)=1（要求A(x)常数项为0）。
2. 递归求解：先递归求出mod x^{ceil(n/2)}下的结果B0(x)。
3. 迭代升级：
   • 计算L(x) = ln(B0(x)) (mod x^n)。注意，这里我们调用poly_ln计算B0的对数，但B0只有前ceil(n/2)项是正确的，poly_ln函数内部会先对其补零到长度n再进行计算。理论上，L(x)的前ceil(n/2)项应等于A(x)的前ceil(n/2)项。
   • 计算D(x) = A(x) - L(x) (mod x^n)。因为B0是mod x^{ceil(n/2)}下的解，所以D(x)的前ceil(n/2)项应为0。
   • 将D(x)的常数项加1：D[0] = (D[0] + 1) % MOD。这对应着迭代式中的(1 - ln(B) + A)。
   • 计算B(x) = B0(x) * D(x) (mod x^n)。
4. 返回结果：返回B(x)。

6.2 复杂度分析与实现陷阱

• 主要开销：一次exp迭代中，包含了一次ln和两次多项式乘法（计算ln内部包含一次求逆和乘法，外部还有一次B0*D的乘法）。因此exp的常数是全家桶中最大的。
• 长度传递：这是最容易出错的地方。在递归调用poly_exp得到B0（长度为ceil(n/2)）后，我们需要将其作为poly_ln的输入。poly_ln要求输入长度是目标长度n，因此我们需要将B0的长度扩展到n（高位补零）。poly_ln输出的结果长度也是n。
• 常数项：必须保证输入A(x)的常数项为0，否则exp结果常数项不为1，与定义不符。
• 迭代的另一种形式：有些实现会将迭代式写为B ≡ B * (A + 1 - ln(B))，本质相同。关键是理解ln(B)的计算是基于当前近似解B0的。

// 伪代码示意：多项式 exp (简化版，展示流程) void poly_exp(int *a, int *b, int n) { static int tmp_ln[MAXN], tmp_d[MAXN]; if (n == 1) { b[0] = 1; return; } poly_exp(a, b, (n+1)/2); // 递归求 B0 // 现在 b 中存储的是长度为 (n+1)/2 的 B0 poly_ln(b, tmp_ln, n); // 计算 ln(B0)，结果长度 n // 计算 D = A - ln(B0) + 1 for (int i=0; i<n; ++i) { tmp_d[i] = (a[i] - tmp_ln[i] + MOD) % MOD; } tmp_d[0] = (tmp_d[0] + 1) % MOD; // 常数项+1 // 计算 B = B0 * D int lim = get_lim(2*n); // 将 b (B0) 补零到长度 lim，将 tmp_d 补零到长度 lim ntt(b, lim, 1); ntt(tmp_d, lim, 1); for (int i=0; i<lim; ++i) b[i] = (ll)b[i] * tmp_d[i] % MOD; ntt(b, lim, -1); for (int i=n; i<lim; ++i) b[i] = 0; // 只保留 n 项 }

7. 实战应用与组合技巧

掌握了这四个基本操作，我们就拥有了强大的多项式处理能力。它们很少单独使用，更多的是组合起来解决复杂问题。

7.1 典型问题建模

1. 生成函数与计数：这是多项式全家桶最经典的应用场景。例如，求某个组合对象的生成函数，其运算可能涉及乘法、求逆（求生成函数的倒数对应某种反演）、exp（比如集合的 exp 对应无序组合）等。

   • 例题：有标号连通图计数。设G(x)是所有有标号图的生成函数，C(x)是所有有标号连通图的生成函数。根据指数生成函数原理，G = exp(C)。已知G（容易计算），则C = ln(G)。这就直接化为了一个多项式ln问题。

2. 多项式复合与快速幂：计算A(x)^k mod x^n。当k很大时，我们可以利用ln和exp：A^k = exp(k * ln(A))。前提是A(x)常数项不为0（通常为1）。这比做k-1次多项式乘法要快得多（O(n log n)vsO(k n log n)）。

3. 多项式三角函数：利用欧拉公式，sin(A(x))和cos(A(x))可以通过exp(i*A(x))和exp(-i*A(x))来表示，其中i是模意义下的单位根（如果模数支持，如 998244353，其i = 86583718）。这又归结到了exp运算。

7.2 实现中的组合调用与优化

在实际的模板代码中，这些函数是相互调用的：

• sqrt调用inv。
• ln调用inv和求导积分。
• exp调用ln。 这意味着，一个exp操作内部，会递归调用ln，而ln又会调用inv。因此，exp的常数非常大，在时间紧张的题目中要谨慎使用。

优化技巧：

• 内存池化：为 NTT 和临时计算预分配全局数组，避免频繁new/delete或vector扩容。
• 逆元预处理：提前预处理1到n的逆元，供ln的积分和求逆中的常数使用。
• 封装与复用：将 NTT 操作、数组拷贝、清零等封装成函数，确保代码清晰且不易出错。
• 长度计算优化：实现一个get_lim函数，根据所需长度快速计算最小的 2 的幂，用于 NTT。

8. 常见问题、调试技巧与心得

即使理解了原理，实现一个健壮高效的多项式全家桶也充满挑战。下面分享一些我踩过的坑和调试经验。

8.1 常见问题速查表

8.2 调试心得与性能压榨

1. 单元测试是王道：不要直接拿复杂题目测试。写一个test()函数，用n=4的小多项式，手动计算出求逆、开根、ln、exp 的预期结果（可以用 Python 的sympy辅助），然后与你的模板输出逐项对比。这是定位问题最有效的方法。

2. 打印中间变量：在怀疑出错的函数里，比如poly_inv的每次迭代后，打印出B数组的前若干项。观察它是否如理论所述快速收敛到正确值。

3. 关注常数项：很多错误都源于常数项处理不当。求逆要求常数项可逆，开根要求常数项是二次剩余，ln要求常数项为1，exp要求常数项为0。这些检查不仅能避免错误，也能帮你快速定位问题阶段。

4. 长度，长度，还是长度：多项式模板 90% 的 bug 都和长度有关。mod x^n意味着只保留前n项。进行 NTT 乘法时，长度必须是2的幂且足够容纳结果（至少deg(A)+deg(B)+1）。每次操作后，都要清晰地知道当前多项式的“有效长度”是多少，并清除无效的高位数据。

5. 空间与时间的权衡：为了极致优化，模板代码往往看起来“脏乱差”，充满了全局数组和指针操作。在竞赛中这是必要的。但在学习和调试阶段，可以先用vector<int>实现一个清晰易懂的版本，确保逻辑正确后，再将其优化为静态数组版本。理解永远比代码风格更重要。

实现一个完全正确且高效的多项式全家桶，就像组装一台精密的机械表。每一个齿轮（函数）都必须严丝合缝，每一次传动（数据传递）都必须精准无误。这个过程充满挑战，但一旦完成，你会发现面对许多复杂的生成函数问题，你手中多了一把万能钥匙。从“黑盒”调用到“白盒”掌控，这种对底层原理的深刻理解，是提升算法能力道路上最坚实的阶梯。