多项式全家桶核心原理:牛顿迭代法统一求逆、开根、ln与exp 1. 项目概述从“黑盒”到“白盒”的多项式运算工具箱在算法竞赛和理论计算机科学领域多项式运算早已不是新鲜话题。从基础的加减乘到稍显复杂的求逆、开根再到更高级的对数ln和指数exp运算这些操作构成了一个被称为“多项式全家桶”的核心工具集。很多选手在初次接触时往往直接套用模板将其视为一个功能强大的“黑盒”——输入系数数组调用函数得到结果。然而当题目条件发生变化或者需要优化常数、处理边界情况时仅仅会“用”是远远不够的。理解其背后的数学原理和实现细节将这个“黑盒”彻底“白盒化”才是从使用者迈向创造者的关键一步。这篇文章就是一次彻底的“白盒化”旅程。我们不满足于仅仅给出代码模板而是要深入拆解多项式求逆、开根、ln、exp这四大核心操作的每一个步骤。我会结合自己多年打比赛和出题的经验详细解释牛顿迭代法是如何在这些运算中统一应用的FFT/NTT快速傅里叶变换/数论变换在其中扮演了什么角色以及那些模板代码里看似神秘的参数和边界处理究竟为何如此设计。无论你是正在备赛的选手还是对算法实现细节有浓厚兴趣的开发者相信这篇超过5000字的深度解析都能让你对“多项式全家桶”有一个全新的、透彻的认识。2. 核心数学原理与统一框架牛顿迭代法在深入每个具体操作之前我们必须先建立一个统一的视角。多项式求逆、开根、ln、exp这些看似不同的运算在算法实现层面其实共享着一个强大的核心工具牛顿迭代法。理解这一点是理解整个“全家桶”的钥匙。2.1 牛顿迭代法在多项式语境下的重塑我们都知道牛顿迭代法用于求实数方程的根从一个初始近似值x0开始通过公式x_{n1} 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和expln要求多项式常数项为1exp要求多项式常数项为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})。实现步骤递归/倍增法边界条件当n1时直接返回常数项a0的逆元。若a0无逆元则整个运算无解。递归求解设当前需要求解mod x^n下的逆。我们先递归求解mod x^{ceil(n/2)}下的逆记为B0(x)。ceil表示向上取整这是为了保证倍增性质。迭代升级计算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下更精确的逆。返回结果返回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 i0; ilim; 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 i2*m; ilim; i) b[i] 0; // 清空高位 } for (int in; ilim; 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。实现步骤边界与初始值当n1时返回常数项a0的模平方根需预处理或使用 Cipolla 算法求解。同样要求a0是二次剩余。递归/倍增求解类似求逆先递归求出mod x^{ceil(n/2)}下的根B0(x)。迭代升级调用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的逆元。返回结果返回B(x)。4.2 实现要点与边界处理依赖求逆这是开根运算的核心特点。在代码组织上开根函数内部会调用求逆函数。常数项处理如果常数项a0不是二次剩余在模意义下无法开根。通常题目会保证有解但自己写代码时要留心。除以2的处理在模MOD下除以2必须转换为乘以 (MOD1)/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, (n1)/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 i0; ilim; i) tmp_a[i] (ll)tmp_a[i] * tmp_inv[i] % MOD; ntt(tmp_a, lim, -1); int inv2 (MOD1)/2; for (int i0; in; i) b[i] (ll)(b[i] tmp_a[i]) * inv2 % MOD; for (int in; ilim; 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符合ln10。计算步骤检查条件确保A(x)的常数项a0 ≡ 1 (mod MOD)。否则ln在形式幂级数意义上无良好定义。求导计算A(x)的导数A(x)。这是一个O(n)的操作A[i] A[i1] * (i1) % MOD。求逆计算A(x)的乘法逆元A_inv(x) (mod x^n)。卷积计算P(x) A(x) * A_inv(x) (mod x^{n-1})。因为A的次数是n-2A_inv我们取前n-1项卷积后我们只需要前n-1项。积分对P(x)进行积分得到结果B(x)。积分公式B[i] P[i-1] * inv(i) % MOD其中inv(i)是i的模逆元且B[0] 0因为ln10。可以看到多项式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 i0; in-1; i) tmp_a[i] (ll)a[i1] * (i1) % 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 i0; ilim; 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 i1; in; i) b[i] (ll)tmp_a[i-1] * inv[i] % MOD; // inv[i] 已预处理 for (int in; ilim; 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)是多项式对数运算。实现步骤倍增法边界与初始值当n1时B(x) ≡ 1 (mod x)因为exp(0)1要求A(x)常数项为0。递归求解先递归求出mod x^{ceil(n/2)}下的结果B0(x)。迭代升级计算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)的常数项加1D[0] (D[0] 1) % MOD。这对应着迭代式中的(1 - ln(B) A)。计算B(x) B0(x) * D(x) (mod x^n)。返回结果返回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, (n1)/2); // 递归求 B0 // 现在 b 中存储的是长度为 (n1)/2 的 B0 poly_ln(b, tmp_ln, n); // 计算 ln(B0)结果长度 n // 计算 D A - ln(B0) 1 for (int i0; in; 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 i0; ilim; i) b[i] (ll)b[i] * tmp_d[i] % MOD; ntt(b, lim, -1); for (int in; ilim; i) b[i] 0; // 只保留 n 项 }7. 实战应用与组合技巧掌握了这四个基本操作我们就拥有了强大的多项式处理能力。它们很少单独使用更多的是组合起来解决复杂问题。7.1 典型问题建模生成函数与计数这是多项式全家桶最经典的应用场景。例如求某个组合对象的生成函数其运算可能涉及乘法、求逆求生成函数的倒数对应某种反演、exp比如集合的 exp 对应无序组合等。例题有标号连通图计数。设G(x)是所有有标号图的生成函数C(x)是所有有标号连通图的生成函数。根据指数生成函数原理G exp(C)。已知G容易计算则C ln(G)。这就直接化为了一个多项式ln问题。多项式复合与快速幂计算A(x)^k mod x^n。当k很大时我们可以利用ln和expA^k exp(k * ln(A))。前提是A(x)常数项不为0通常为1。这比做k-1次多项式乘法要快得多O(n log n)vsO(k n log n)。多项式三角函数利用欧拉公式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 常见问题速查表问题现象可能原因排查方法结果全为0或明显错误1. NTT 的len或lim计算错误。2. 忘记在 NTT 前后进行位逆序置换。3. 模数MOD或原根G写错。1. 打印每次 NTT 调用时的长度lim。2. 检查 NTT 的rev数组是否正确初始化。3. 用简单数据如{1, 1}测试 NTT 正逆变换。求逆或开根结果前几项对后面错1. 迭代后没有正确清零高位系数mod x^n操作未执行。2. 递归/倍增边界处理错误长度传递混乱。1. 在每次迭代或乘法后手动将n之后的系数置0。2. 仔细检查递归函数中传入的长度和需要的长度是否匹配。ln或exp结果爆炸或溢出1.ln的输入多项式常数项不为1。2.exp的输入多项式常数项不为0。3. 积分时使用的逆元inv[i]计算错误。1. 在ln开头检查a[0] 1。2. 在exp开头检查a[0] 0。3. 验证预处理逆元数组的正确性。答案与暴力计算对不上1. 题目模数不是 NTT 友好模数需要三模 NTT 或 MTT。2. 多项式长度超过 NTT 能处理的范围lim太大。3. 运算顺序或公式推导有误。1. 确认模数如1e97需用三模 NTT。2. 用小的n如4进行单元测试打印每一步的中间结果与手算对比。8.2 调试心得与性能压榨单元测试是王道不要直接拿复杂题目测试。写一个test()函数用n4的小多项式手动计算出求逆、开根、ln、exp 的预期结果可以用 Python 的sympy辅助然后与你的模板输出逐项对比。这是定位问题最有效的方法。打印中间变量在怀疑出错的函数里比如poly_inv的每次迭代后打印出B数组的前若干项。观察它是否如理论所述快速收敛到正确值。关注常数项很多错误都源于常数项处理不当。求逆要求常数项可逆开根要求常数项是二次剩余ln要求常数项为1exp要求常数项为0。这些检查不仅能避免错误也能帮你快速定位问题阶段。长度长度还是长度多项式模板 90% 的 bug 都和长度有关。mod x^n意味着只保留前n项。进行 NTT 乘法时长度必须是2的幂且足够容纳结果至少deg(A)deg(B)1。每次操作后都要清晰地知道当前多项式的“有效长度”是多少并清除无效的高位数据。空间与时间的权衡为了极致优化模板代码往往看起来“脏乱差”充满了全局数组和指针操作。在竞赛中这是必要的。但在学习和调试阶段可以先用vectorint实现一个清晰易懂的版本确保逻辑正确后再将其优化为静态数组版本。理解永远比代码风格更重要。实现一个完全正确且高效的多项式全家桶就像组装一台精密的机械表。每一个齿轮函数都必须严丝合缝每一次传动数据传递都必须精准无误。这个过程充满挑战但一旦完成你会发现面对许多复杂的生成函数问题你手中多了一把万能钥匙。从“黑盒”调用到“白盒”掌控这种对底层原理的深刻理解是提升算法能力道路上最坚实的阶梯。