尧图网站设计 尧图网站设计YAOTU DESIGN
ARTICLE DETAIL

资讯详情

深耕网站设计与一线实操的经验洞察。

多项式对数函数(ln)算法详解:从公式推导到NTT实现与调试

多项式对数函数(ln)算法详解:从公式推导到NTT实现与调试 1. 从一道模板题说起多项式对数函数ln到底是什么如果你在洛谷、Codeforces或者任何一个算法竞赛社区混迹过一段时间大概率会刷到过“P4725 【模板】多项式对数函数多项式 ln”这道题。它就像算法竞赛选手在多项式领域的一个“成人礼”标志着从只会做加减乘除的“小学生”进阶到开始触及多项式更深刻运算的“中学生”。但很多人在第一次接触时都会懵多项式还能取对数这玩意儿有什么用难道是把1 2x 3x^2丢进计算器按ln键吗显然不是。这里的“多项式对数函数”是一个形式幂级数上的形式运算。它解决的核心问题是给定一个常数项为1的多项式或形式幂级数A(x)求另一个多项式B(x)使得在形式幂级数的意义下exp(B(x)) A(x)。这里 exp 是指数函数。换句话说B(x) 就是 A(x) 的“形式对数”。这个运算在组合数学、生成函数、多项式算法中有着极其重要的地位。比如当你用生成函数刻画一个组合结构时对其取 ln 往往对应着将连通分量拆解出来在多项式牛顿迭代求解中ln 和 exp 是一对关键的基础算子。这道模板题之所以经典是因为它完美地将多项式求导、积分、求逆、乘法这几个基础操作串联了起来形成了一个完整的算法链条。网上能找到的题解和代码很多但大多只给出了“怎么做”的步骤和代码对于“为什么这么做”、“每一步背后的数学原理是什么”、“实现时有哪些一踩就炸的坑”却语焉不详。我这篇文章就想结合我多次实现和调试的经验把这些隐藏在水面下的东西彻底讲透。我们不止要会套模板更要理解这个模板的每一颗螺丝钉是怎么拧上去的。2. 核心公式推导为什么求ln变成了求导、求逆和积分几乎所有教程都会直接甩给你这个公式 若 A(x) 1 a_1 x a_2 x^2 ...且 A(0)1则ln(A(x)) ∫ [A(x) / A(x)] dx这个公式是整套算法的基石。我们来一步步拆解它看它到底是怎么来的。2.1 从形式微分的定义出发首先我们得认同对形式幂级数也可以定义“导数”。这很直观对于多项式 A(x) ∑_{i0}^{n} a_i x^i其形式导数 A(x) 就是 ∑_{i1}^{n} i * a_i x^{i-1}。就是把每一项的指数拿下来当系数然后指数减一。现在考虑我们想求的 B(x) ln(A(x))。这里 ln 是一个形式运算。我们对这个等式两边同时关于 x 求形式导数利用链式法则左边B(x) 右边d/dx [ln(A(x))] A(x) / A(x) 这里直接类比了实数域上 ln(f(x)) 的导数为 f(x)/f(x)于是我们得到了一个关键等式B(x) A(x) / A(x)。2.2 从微分到积分得到了 B(x) 的导数那么 B(x) 本身自然就是对其积分B(x) ∫ B(x) dx ∫ [A(x) / A(x)] dx注意这里积分会有一个积分常数 C。因为是不定积分。那么 C 是多少我们需要利用初始条件A(0) 1。我们希望 B(x) 也是一个形式幂级数并且通常定义 ln(1) 0。所以 B(0) ln(A(0)) ln(1) 0。另一方面我们对 ∫ [A(x) / A(x)] dx 求出的结果其常数项就是积分产生的常数 C。为了让 B(0) 0我们必须令 C 0。所以最终公式里我们直接写为定积分形式从 0 积到 x或者理解为取不定积分后忽略常数项因为常数项为0。在实现时我们做不定积分然后手动将结果的常数项设为0即可。2.3 公式的可行性分析这个公式将 ln 运算转化为了三个我们已知能做的操作求导 (A(x))O(n) 复杂度极其简单。求逆 (1 / A(x))这里需要计算 A(x) 的乘法逆元。这需要用到多项式求逆算法通常使用牛顿迭代法复杂度 O(n log n)。积分 (∫ ... dx)O(n) 复杂度是求导的逆过程同样简单。所以整个多项式 ln 的算法复杂度就卡在了多项式求逆这一步为 O(n log n)。这也就是为什么多项式求逆是多项式全家桶里更基础的一个模板。注意这个公式成立有一个绝对的前提A(x) 的常数项必须为 1。为什么 从数学上看ln(A(x)) 要想展开成形式幂级数必须在 x0 处有定义且 ln(A(0)) 需要是一个有限值我们取0。如果 A(0)0ln(0) 无定义如果 A(0) 是其他非1常数 c那么 ln(A(x)) ln(c) ln(1 (A(x)-c)/c)。这里 ln(c) 是一个实数常数但我们的多项式是在某个模数如998244353的有限域上运算的ln(c) 在这个域里可能没有定义除非 c 是模数的原根相关。为了保证运算纯粹在模意义下进行且结果是一个多项式常数项为0最方便且通用的约定就是要求 A(0)1。这样 ln(1)0一切都很干净。3. 算法步骤拆解与零基础实现指南理解了公式我们来把算法步骤彻底细化。假设我们有多项式 A(x)其次数为 n-1通常我们处理长度为 n 的数组下标 0 到 n-1 对应次数 0 到 n-1 的系数且满足 A[0] 1。3.1 第一步计算 A(x) 的导数 A(x)这一步是热身。设 A(x) a0 a1x a2x^2 ... a_{n-1}x^{n-1}。 那么 A(x) a1 2a2x 3a3*x^2 ... (n-1)*a_{n-1}*x^{n-2}。在代码中这就是一个简单的循环// 假设系数存储在数组 a 中长度为 n for (int i 1; i n; i) { da[i-1] 1LL * a[i] * i % mod; // da 存储导数系数 } // da 的有效长度变为 n-1注意边界导数的次数比原多项式低一次。3.2 第二步计算 A(x) 的乘法逆元 B(x) 1 / A(x)这是整个算法的核心和性能瓶颈。我们需要求一个多项式 B(x)使得 A(x) * B(x) ≡ 1 (mod x^n)。这里mod x^n的意思是我们只关心乘积的前 n 项0 到 n-1 次更高次的项可以忽略。求逆通常使用牛顿迭代法。其思想是假设我们已经求出了在模 x^{ceil(m/2)} 意义下的逆元 B_0(x)如何快速得到在模 x^m 意义下的逆元 B(x)推导过程略涉及泰勒展开结论是迭代公式为B(x) ≡ B_0(x) * (2 - A(x) * B_0(x)) (mod x^m)实际操作时我们采用递归或迭代倍增的方式初始条件当 n1 时A(x) 只有一个常数项 a0。由前提 a0 1所以在模 x^1 意义下其逆元就是 1。假设我们已经求出在模 x^{ceil(n/2)} 意义下的逆元 B_0(x)。目标计算模 x^n 意义下的逆元 B(x)。根据公式我们需要计算T(x) A(x) * B_0(x) (mod x^n)// 注意这里模数要提升到 x^nT(x) (2 - T(x)) (mod x^n)// 对 T(x) 的每一项做 2 - t_i 运算B(x) T(x) * B_0(x) (mod x^n)这个过程需要多项式乘法。利用 NTT快速数论变换可以将乘法优化到 O(n log n)。由于这部分是独立模板代码较长。其关键点在于每次迭代时对于 A(x) 我们只需要前 m 项当前目标长度对于 B_0(x) 我们知道它在模 x^{m/2} 下是精确的。计算A(x) * B_0(x)时结果长度会增长但我们只取前 m 项。然后进行2 - T(x)的系数运算最后再乘一次 B_0(x) 并取前 m 项。3.3 第三步计算 C(x) A(x) * B(x)现在我们有了导数da长度为 n-1和逆元b长度为 n。将它们相乘。注意da的次数是 n-2b的次数是 n-1它们的乘积次数最高为 (n-2)(n-1)2n-3。但我们最终只需要前 n-1 项因为下一步积分后我们要得到 n 项结果。所以我们可以只计算到长度至少为 n-1 的卷积。设dc da * b。我们取dc的前 n-1 项。注意dc[0]对应的是A(x)*B(x)的常数项。3.4 第四步对 C(x) 积分得到最终结果积分是导数的逆运算。如果C(x) c0 c1*x c2*x^2 ...那么它的积分∫ C(x) dx C c0*x (c1/2)*x^2 (c2/3)*x^3 ...其中 C 是积分常数。我们已经知道结果的常数项必须为 0。所以我们计算 对于 i 从 0 到 n-2res[i1] dc[i] * inv(i1) % mod其中inv(i1)是 i1 在模 mod 下的乘法逆元需要预处理。 而res[0] 0。这样得到的res就是ln(A(x))的前 n 项系数。3.5 完整流程图示与复杂度分析输入: A(x), 满足 A[0]1, 次数界 n 输出: B(x) ln(A(x)) mod x^n 1. 求导: DA(x) derivative(A(x)) // O(n) 2. 求逆: IA(x) inverse(A(x), n) // O(n log n) 使用牛顿迭代NTT 3. 乘法: C(x) DA(x) * IA(x) mod x^{n-1} // O(n log n) NTT乘法结果取前n-1项 4. 积分: B(x) integral(C(x)) // O(n) 常数项设为0 5. 返回 B(x)总时间复杂度由两次 O(n log n) 的操作主导即求逆和乘法。空间上需要一些临时数组进行变换和计算。4. 实战代码剖析从模块构建到边界处理光说不练假把式。下面我结合一个典型的基于 NTT模数 998244353原根为 3的实现来逐块解析代码并指出那些容易写错、调试到崩溃的细节。4.1 基础工具函数快速幂与逆元const int mod 998244353, g 3; // 原根 int qpow(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; }qpow是标准快速幂。inv函数通常直接调用qpow(a, mod-2)但频繁调用时建议预处理 1~n 的逆元。4.2 核心NTT 与多项式乘法这是所有多项式操作的基础。代码较长但结构固定。关键点在于rev数组的蝴蝶变换以及三层循环的迭代实现。这里我强调几个易错点长度与界限NTT 要求长度是 2 的幂。每次进行多项式乘法前必须计算lim 1, bit 0; while (lim n m) lim 1, bit;。然后初始化 rev 数组。结果清零对于长度为lim的数组一定要确保lim范围内的数据是有效的或者在使用前清空。特别是多次调用时旧数据可能残留。逆变换后的缩放NTT 逆变换后每一项需要乘以lim的逆元。int invlim qpow(lim, mod-2); for (int i0; ilim; i) a[i]1LL*a[i]*invlim%mod;4.3 多项式求逆的实现细节这是最难写对的部分。我给出一个相对清晰的迭代版本框架void poly_inv(int *a, int *b, int n) { // 计算 b(x)使得 a(x)*b(x) ≡ 1 (mod x^n) static int tmp[N]; // 临时数组需要足够大如4倍n b[0] qpow(a[0], mod-2); // 初始条件常数项逆元 for (int len 2; (len 1) n; len 1) { // 当前目标是求出模 x^len 下的逆元 int lim len 1; // 乘法需要长度 // 将 a 的前 len 项拷贝到 tmp并做 NTT for (int i 0; i len; i) tmp[i] i n ? a[i] : 0; for (int i len; i lim; i) tmp[i] b[i] 0; ntt(tmp, lim, 1); ntt(b, lim, 1); // 根据公式 B B0 * (2 - A * B0) 计算 for (int i 0; i lim; i) { b[i] 1LL * b[i] * (2 - 1LL * tmp[i] * b[i] % mod mod) % mod; } ntt(b, lim, -1); // 重要b 中 len 之后的项是无意义的必须清零防止影响下一轮 for (int i len; i lim; i) b[i] 0; } // 最后确保 b 只有前 n 项有效后面清零如果传入的n不是2的幂 for (int i n; i lim; i) b[i] 0; }踩坑实录1清零清零清零这是多项式题最经典的错误。在牛顿迭代的每一轮结束后b数组在len之后的系数必须手动设为0。因为 NTT 逆变换后这些位置可能留有上一轮或计算过程中的垃圾值。下一轮循环时我们会把整个b数组包括后面的垃圾值做 NTT这些垃圾值会污染整个频域导致结果完全错误。这个 bug 非常隐蔽因为小数据时可能因为长度不够碰不到垃圾内存而侥幸正确大数据一定挂。4.4 多项式 ln 的完整实现集成了求导、求逆、乘法和积分。void poly_derivative(int *a, int *da, int n) { for (int i 1; i n; i) da[i-1] 1LL * a[i] * i % mod; da[n-1] 0; // 导数长度减一最后一位可置0 } void poly_integral(int *a, int *ia, int n) { ia[0] 0; // 常数项为0 // 预处理1~n的逆元 inv[i] for (int i 1; i n; i) ia[i] 1LL * a[i-1] * inv[i] % mod; } void poly_ln(int *a, int *res, int n) { // 前提检查a[0] 必须为 1 assert(a[0] 1); static int da[N], ia[N], tmp[N]; // 1. 求导 poly_derivative(a, da, n); // da 长度为 n-1 // 2. 求逆 poly_inv(a, ia, n); // ia 是 A(x) 的逆长度为 n // 3. 乘法da * ia int lim 1, bit 0; while (lim (n-1) n) lim 1, bit; // 结果需要前 n-1 项 for (int i 0; i lim; i) tmp[i] (i n-1) ? da[i] : 0; for (int i 0; i lim; i) { // 注意ia 只有前 n 项有效后面在 poly_inv 中已清零 // 但这里为了安全可以只拷贝前 n 项后面置0 if (i n) ia[i] ia[i]; else if (i lim) ia[i] 0; // 确保 ia 在 lim 长度内有效 } // 这里需要调用一个标准的 NTT 乘法函数输入 tmp 和 ia结果存回 tmp ntt_mul(tmp, ia, lim); // 假设这个函数处理了 NTT 变换、点乘、逆变换和缩放 // 现在 tmp 的前 n-1 项是 da * ia 的结果 // 4. 积分 poly_integral(tmp, res, n); // 积分后长度变为 n }踩坑实录2长度对齐与数组越界在调用poly_inv时我们传入的长度是n它会计算模x^n的逆。但注意poly_inv内部可能按2的幂分配内存。我们传给poly_ln的数组a其有效长度就是n。但在求导后da有效长度是n-1。进行乘法da * ia时da长度n-1ia长度n卷积结果长度至少需要(n-1)(n-1)2n-2才能保证前n-1项精确。我们设置的lim必须大于等于这个值。同时要确保传入 NTT 乘法的数组在lim长度内都有定义要么是有效系数要么是0。任何未初始化的值都会导致错误。5. 调试技巧与常见问题排查就算你完全理解了算法第一遍代码也几乎不可能一次 AC。以下是我总结的排查清单5.1 结果完全不对/随机数检查 NTT 的正确性这是根源。写一个简单的测试比如计算(12x) * (13x)看结果是不是1 5x 6x^2。确保正变换、点乘、逆变换、缩放每一步都正确。检查数组清零如 4.3 节所述在poly_inv的每一轮迭代后必须清零b数组len之后的部分。在poly_ln中调用 NTT 前确保tmp和ia在lim范围内的数据是干净的。检查长度计算lim是否足够大while (lim n m) lim 1中的n和m是否正确对于da * ian n-1da的长度m nia的长度所以lim需要至少2n-1的下一个2的幂。5.2 结果前几项对后面错检查求逆的边界poly_inv函数是否保证了结果严格只有前n项有效在倍增过程中我们计算的是模x^len的逆但len可能大于n。函数最后需要把n之后的系数清零。检查积分用的逆元表inv[i]是否预处理正确inv[i]是i在模mod下的逆元通常用线性递推inv[i] mod - 1LL * (mod/i) * inv[mod%i] % mod来求。确保inv[1] 1。5.3 常数项不为0检查输入确认输入多项式a[0]是否真的为 1。题目可能不保证需要自己先判断或处理。检查积分函数poly_integral是否将res[0]设为了 05.4 性能问题TLENTT 的蝴蝶变换rev数组是否预处理每次乘法都重新计算会超时。不必要的拷贝在poly_inv和poly_ln中尽量减少大数组的memcpy操作。使用指针和就地计算。乘法优化对于da * ia我们只需要前n-1项。可以使用“半在线卷积”的思路进行优化但模板题通常不需要标准的 NTT 乘法即可通过。5.5 一个实用的调试方法对拍写一个暴力版本的poly_ln用于小数据范围比如 n 10的验证。 暴力版本可以模拟形式幂级数的运算先预处理逆元然后根据定义ln(A(x)) ∑_{k1} (-1)^{k-1} * (A(x)-1)^k / k。因为 A(x)-1 的常数项为0所以这个级数在模 x^n 意义下是有限的只需要算到 kn-1 即可。 用这个暴力程序去验证你的 NTT 优化版本在小数据n5,6,7...下的结果是否一致。这是定位问题最有效的方式。6. 从模板到应用ln 在生成函数中的意义搞懂了实现我们再来聊聊它到底有什么用这样下次遇到问题你才能想到用它。6.1 组合意义的连接集合与连通分量这是 ln 最经典的应用。假设我们有一个组合类 A比如所有的图其指数生成函数EGF为 A(x)。那么A(x)的 expB(x) exp(A(x))通常代表了由 A 中的“连通”对象任意组合而成的“所有”对象。例如A(x) 是连通图的 EGF那么 B(x) 就是所有图的 EGF。反过来A(x)的 lnC(x) ln(B(x))就代表了从“所有”对象中提取出“连通”分量。例如已知所有图的 EGF B(x)那么 ln(B(x)) 就是连通图的 EGF。在很多计数问题中我们更容易求出所有方案的生成函数而想要得到连通方案的生成函数就需要对其取 ln。6.2 多项式牛顿迭代中的角色牛顿迭代是求解多项式方程F(G(x)) 0的强大工具。例如求exp、求sqrt开根、求复合逆函数等。 在推导这些迭代式时ln和exp经常作为一对互逆的运算出现。例如求G(x) exp(F(x))可以转化为方程ln(G(x)) - F(x) 0然后应用牛顿迭代。此时poly_ln就成了迭代过程中必须调用的子程序。6.3 形式微分与形式积分的工具ln的公式本身完美结合了求导、求逆和积分。这使得它成为学习多项式形式运算的一个优秀案例。掌握了它你就掌握了处理形式幂级数的一整套基本工具链。7. 总结与扩展思考实现一个poly_ln就像搭乐高。你需要准备好“求导”、“求逆”其内部又需要“NTT乘法”、“积分”这几个基础模块然后按照公式∫ (A / A) dx把它们正确地拼接起来。其中多项式求逆是最复杂、最容易出错的一环务必理解其牛顿迭代的倍增思想并牢记迭代后清零的纪律。在竞赛中poly_ln很少单独出题它往往是更大问题的一块拼图。比如你需要先对某个生成函数取 ln进行一些操作再取 exp。因此将它写对、写熟封装成一个可靠的函数是进军更高级多项式算法如指数函数、三角函数、快速幂、复合逆的必经之路。最后关于常数项不为1的情况理论上可以通过提取公因式解决若 A(0) c ≠ 0则 ln(A(x)) ln(c) ln(A(x)/c)。但 ln(c) 在模意义下需要离散对数来求解这超出了普通多项式模板的范围。所以模板题和常见应用都默认常数项为1。写多项式代码是对耐心和细心的双重考验。一个符号的错误、一次忘记的清零都可能导致调试数小时。但一旦你彻底征服了它那种对复杂算法了如指掌、对每一行代码都充满自信的感觉是无与伦比的。希望这篇超详细的拆解能帮你少走些弯路真正把这块硬骨头啃下来。
返回列表