1. 项目缘起从一道“卡脖子”的模数组合数说起几年前我在准备一场算法竞赛时遇到了一道让我记忆犹新的题目。题目本身描述很简单给定一个巨大的组合数 C(n, m)以及一个模数 P要求计算 C(n, m) mod P 的值。这听起来像是数论基础题我信心满满地写下了标准的预处理阶乘和逆元的代码。然而当我提交时系统返回了一个大大的“Wrong Answer”。仔细一看模数 P 的描述心里顿时凉了半截——这个 P 不是一个质数甚至不是一个质数的幂而是一个任意的正整数比如 10007 或者 999983 这种质数还好但题目给的可能是 1000、2016 甚至 1000000007 * 1000000009 这种合数。这就是经典的“组合数取模”问题在模数非质数时遇到的困境。我们熟悉的用费马小定理或扩展欧几里得求逆元的方法其前提是模数为质数这样才能保证在模意义下每个非零数都有乘法逆元。当模数是合数时许多数没有逆元整个基于阶乘和逆元的递推公式就失效了。我当时卡在这道题上很久直到后来系统学习了“扩展卢卡斯定理”Extended Lucas Theorem才豁然开朗。而 P4720正是洛谷上一道专门练习这个定理的经典模板题。今天我就结合自己踩坑和实战的经验把这个强大工具的原理、实现细节和避坑指南掰开揉碎了讲清楚。2. 核心问题拆解为什么普通卢卡斯定理不够用在深入扩展卢卡斯之前我们必须先理解普通卢卡斯定理Lucas Theorem的局限性这样才能明白我们到底要解决什么问题。2.1 卢卡斯定理的适用场景与限制普通卢卡斯定理表述为对于质数 p有 C(n, m) ≡ C(n mod p, m mod p) * C(n/p, m/p) (mod p)。这是一个递归公式能将大规模组合数计算分解为小规模组合数计算通常配合预处理小范围的阶乘和逆元来快速求解。它的效率很高代码也很简洁。然而它的核心限制就藏在前提里模数 p 必须是质数。这是因为在递归的底层我们需要计算 C(n‘, m’) mod p这通常通过公式 C n! / (m! * (n-m)!) 来计算而除法在模运算中需要转化为乘以其乘法逆元。逆元存在的充要条件就是该数与模数互质。当模数是质数 p 时只要分母的阶乘不被 p 整除其逆元就一定存在。但一旦模数 p 是合数分母的阶乘很可能与模数有公因子导致逆元不存在整个计算链就断裂了。2.2 合数模数带来的真正挑战当模数 P 是合数时直接计算 C(n, m) mod P 的难点可以归结为两点非互质导致的逆元缺失在计算 n! / (m! * (n-m)!) 时分母可能与模数 P 不互质因此无法直接求逆元进行模除。模数非质数中国剩余定理CRT成为桥梁解决这个问题的核心思路是将合数模数 P 质因数分解为 P p1^k1 * p2^k2 * ... * pt^kt。如果我们能分别求出 C(n, m) 对每个质数幂模数 pi^ki 的余数 ai即求解一系列同余方程x ≡ a1 (mod p1^k1) x ≡ a2 (mod p2^k2) ... x ≡ at (mod pt^kt) 那么根据中国剩余定理我们就可以唯一确定出 x 在模 P 意义下的值。所以问题的关键转化为如何计算 C(n, m) mod p^k其中 p 是质数k是正整数。这就是扩展卢卡斯定理要解决的核心子问题。普通卢卡斯定理处理的是 mod pk1而扩展卢卡斯将其推广到了 mod p^k。3. 扩展卢卡斯定理的核心原理剥离p因子与递归求解计算 C(n, m) mod p^k 不能直接用阶乘逆元因为分母可能包含因子 p导致与模数 p^k 不互质。扩展卢卡斯定理的精妙之处在于它通过一种“剥离”技巧将阶乘中所有 p 的因子分离出来单独处理。3.1 第一步将阶乘分解为“与p互质部分”和“p的幂次部分”定义函数F(n, p, pk)用于计算 n! 中所有与 p 互质的因子的乘积再对 pk (即 p^k) 取模。同时我们记录下 n! 中 p 这个质因子的总次数记为G(n, p)。以 n22, p3 为例计算 22! mod 3^2 22! 1 * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10 * 11 * 12 * 13 * 14 * 15 * 16 * 17 * 18 * 19 * 20 * 21 * 22 我们可以把它重写为 22! (124578101113141617192022) * (36912151821) 进一步把第二组每个数中的因子3提出来 22! (124578101113141617192022) * 3^7 * (1234567) 你会发现(124578101113141617192022) 这些数都与3互质而 (1234567) 正好是 floor(22/3) 7 的阶乘即 7!。于是我们得到一个递归定义 n! ≡ F(n, p, pk) * p^{G(n, p)} * (n/p)! (mod pk) 其中F(n, p, pk)计算了1到n中所有不被p整除的数的乘积模 pk。G(n, p) floor(n/p) floor(n/p^2) floor(n/p^3) ...即n!中质因子p的个数。(n/p)!是递归部分。对于F(n, p, pk)的计算也有技巧。因为模数是 pk而1到pk中与p互质的数会形成一个长度为 φ(pk)pk-p^{k-1} 的循环节。我们可以先计算一个完整循环节内所有与p互质的数的乘积模 pk记为prod。那么n 以内这样的完整循环节有n / pk个每个循环节的贡献是prod^{n/pk} mod pk。最后再乘以剩余的不完整部分即从floor(n/pk)*pk 1到 n 之间且与p互质的数的乘积。3.2 第二步计算组合数模 p^k有了上面的分解组合数可以表示为 C(n, m) n! / (m! * (n-m)!) 将其用 F 和 G 函数表示 C(n, m) [F(n) * p^{G(n)}] / [F(m) * p^{G(m)} * F(n-m) * p^{G(n-m)}] [F(n) / (F(m) * F(n-m))] * p^{G(n) - G(m) - G(n-m)}我们的目标是求 C(n, m) mod p^k。首先计算指数部分e G(n) - G(m) - G(n-m)。如果 e k说明组合数本身包含了至少 p^k 这个因子那么 C(n, m) mod p^k 0。如果 e k则继续。计算互质部分num F(n, p, pk) * inv(F(m, p, pk), pk) * inv(F(n-m, p, pk), pk) mod pk。这里inv(a, pk)表示 a 在模 pk 意义下的逆元。由于 F 函数计算的结果都是与 p 互质的所以它们对模数 pk 的逆元一定存在可以用扩展欧几里得算法求解。最终结果C(n, m) mod p^k num * p^e mod pk。3.3 第三步中国剩余定理CRT合成最终答案假设我们对合数模数 P 分解后得到了 t 个方程 x ≡ ans_i (mod pi^ki), i1, 2, ..., t。 其中 ans_i 就是我们用上述方法计算出来的 C(n, m) mod pi^ki。中国剩余定理的求解过程如下计算M P。对于每个 i计算Mi M / (pi^ki)。计算Mi在模pi^ki意义下的逆元inv_i因为 Mi 与 pi^ki 互质逆元存在。最终解为x Σ(ans_i * Mi * inv_i) mod M。这一步在算法实现中通常使用“增量法”合并同余方程每次合并两个方程逐步得到最终解比一次性计算所有逆元更易于编码。4. 手把手实现扩展卢卡斯模板理解了原理我们来看代码实现。我将结合 P4720 这道模板题的要求给出一个清晰、健壮且包含详细注释的 C 实现。代码会分为几个核心函数。4.1 基础工具函数快速幂与扩展欧几里得这些是数论算法的基石。// 快速幂计算 (base^exp) % mod long long qpow(long long base, long long exp, long long mod) { long long res 1 % mod; // 注意 mod1 的情况 base % mod; while (exp) { if (exp 1) res (res * base) % mod; base (base * base) % mod; exp 1; } return res; } // 扩展欧几里得算法求解 ax by gcd(a, b) // 返回 gcd(a, b)并通过引用返回 x, y long long exgcd(long long a, long long b, long long x, long long y) { if (b 0) { x 1; y 0; return a; } long long d exgcd(b, a % b, y, x); y - (a / b) * x; return d; } // 求 a 在模 mod 下的逆元前提是 gcd(a, mod) 1 long long inv(long long a, long long mod) { long long x, y; exgcd(a, mod, x, y); // 将逆元调整到 [0, mod) 范围内 return (x % mod mod) % mod; }4.2 核心函数 F计算剔除了p因子的阶乘模 p^k这个函数对应原理部分的F(n, p, pk)。/** * 计算 n! 中所有与质数 p 互质的因子的乘积再对 pk (p^k) 取模。 * param n 阶乘的上限 * param p 质数 * param pk p^k即当前处理的质数幂模数 * return n! 中与 p 互质部分的乘积模 pk */ long long factorial_prime(long long n, long long p, long long pk) { if (n 0) return 1; long long res 1; // 1. 处理完整循环节周期为 pk每个周期内与p互质的数乘积相同 // 计算一个周期内的乘积 long long cycle_prod 1; for (long long i 1; i pk; i) { if (i % p ! 0) { // 与p互质 cycle_prod (cycle_prod * i) % pk; } } // 共有 n/pk 个完整周期 res qpow(cycle_prod, n / pk, pk); // 2. 处理最后一个不完整的周期 for (long long i (n / pk) * pk 1; i n; i) { if (i % p ! 0) { res (res * (i % pk)) % pk; // i % pk 防止溢出且结果等价 } } // 3. 递归处理 (n/p)! 中与p互质的部分 // 因为 n! (1*2*...*n) (所有与p互质的数) * p * (所有与p互质的数) * ... // 递归部分正是 (n/p)! 中与p互质的部分 return res * factorial_prime(n / p, p, pk) % pk; }注意这里有一个非常关键的优化和易错点。在计算cycle_prod时我们是在模pk下计算1到pk之间与p互质的数的乘积。pk可能很大比如 5^10直接循环pk次在n很大时会被多次调用可能成为性能瓶颈。在实际的高性能模板中通常会预处理这个值。但为了代码清晰这里展示了最基本的逻辑。在真正解题时需要根据数据范围权衡是否预处理。4.3 核心函数 G计算 n! 中质因子 p 的个数这个函数用于计算指数e。/** * 计算 n! 中质因子 p 的个数。 * 公式G(n, p) floor(n/p) floor(n/p^2) floor(n/p^3) ... * param n 阶乘的上限 * param p 质数 * return n! 中质因子 p 的个数 */ long long count_prime_factor(long long n, long long p) { long long cnt 0; while (n) { cnt n / p; n / p; } return cnt; }4.4 核心函数 C_mod_pk计算组合数模质数幂这个函数整合前两步计算 C(n, m) mod p^k。/** * 计算组合数 C(n, m) 对质数幂 p^k 取模的结果。 * param n 组合数上标 * param m 组合数下标 * param p 质数 * param pk p^k * return C(n, m) mod pk */ long long combination_mod_prime_power(long long n, long long m, long long p, long long pk) { if (m n) return 0; if (m 0 || m n) return 1 % pk; // 1. 计算指数 e G(n) - G(m) - G(n-m) long long e count_prime_factor(n, p) - count_prime_factor(m, p) - count_prime_factor(n - m, p); if (e 0) { // 实际上e永远0这里判断是为了逻辑清晰也可以判断 if (e k) return 0; // 如果 e k (即pk中p的幂次)那么结果模pk为0。这里k可以通过pk和p计算得到。 // 简便写法如果 e 足够大使得 p^e 是 pk 的倍数则直接返回0。 // 更严谨的做法是计算 k log(pk) / log(p) 的整数部分然后比较 e 和 k。 // 下面我们采用另一种方式先计算互质部分如果 e 很大最后乘 p^e 时再取模。 } else { // 理论上不会出现负数因为组合数是整数 return 0; } // 2. 计算互质部分F(n) / (F(m) * F(n-m)) long long fn factorial_prime(n, p, pk); long long fm factorial_prime(m, p, pk); long long fnm factorial_prime(n - m, p, pk); long long num fn * inv(fm, pk) % pk * inv(fnm, pk) % pk; // 3. 乘以 p^e long long pe qpow(p, e, pk); // 注意这里是对 pk 取模因为最终结果是模 pk // 但是当 e k 时p^e mod pk 为 0所以这一步包含了 ek 时结果为0的情况。 return num * pe % pk; }踩坑点在计算num时一定要先对fm和fnm分别求逆元然后连乘取模。不能先计算fm * fnm % pk再求一次逆元因为乘法可能破坏互质性导致逆元不存在。必须保证每个与pk互质的数单独求逆。4.5 核心函数 exLucas主函数与CRT合并这是对外的接口处理合数模数 P。/** * 扩展卢卡斯定理主函数计算 C(n, m) mod PP 可为任意正整数。 * param n 组合数上标 * param m 组合数下标 * param P 模数任意正整数 * return C(n, m) mod P */ long long exLucas(long long n, long long m, long long P) { if (m n) return 0; if (m 0 || m n) return 1 % P; long long mod P; // 存储 (质数, 质数幂, 余数) 三元组 vectortuplelong long, long long, long long factors; // 1. 对模数 P 进行质因数分解 for (long long i 2; i * i mod; i) { if (mod % i 0) { long long pk 1; while (mod % i 0) { mod / i; pk * i; } // 计算 C(n, m) mod i^pk long long res combination_mod_prime_power(n, m, i, pk); factors.emplace_back(i, pk, res); } } if (mod 1) { // 处理剩余的大质数 factors.emplace_back(mod, mod, combination_mod_prime_power(n, m, mod, mod)); } // 2. 如果只有一个质因数直接返回结果 if (factors.size() 1) { return get2(factors[0]); } // 3. 使用中国剩余定理CRT合并所有同余方程 // 增量法合并x ≡ a1 (mod m1), x ≡ a2 (mod m2) // 合并为 x ≡ new_a (mod new_m)其中 new_m m1 * m2 long long a1 get2(factors[0]), m1 get1(factors[0]); for (size_t i 1; i factors.size(); i) { long long a2 get2(factors[i]), m2 get1(factors[i]); // 合并方程x a1 k1*m1 a2 k2*m2 // 即 k1*m1 - k2*m2 a2 - a1 // 令 g gcd(m1, m2)用扩展欧几里得求解 long long k1, k2; long long g exgcd(m1, m2, k1, k2); long long c a2 - a1; if (c % g ! 0) { // 理论上不会发生因为各模数两两互质 return -1; // 无解 } long long t m2 / g; // 调整 k1 为最小非负特解 k1 (k1 * (c / g) % t t) % t; // 新的余数和模数 long long new_a (a1 k1 * m1) % (m1 / g * m2); // 注意新模数是 lcm(m1, m2) m1/g*m2 long long new_m m1 / g * m2; a1 new_a; m1 new_m; } return (a1 % P P) % P; // 确保结果在 [0, P) 范围内 }5. 实战测试与性能优化要点将上述代码整合就可以通过 P4720 这道模板题了。输入 n, m, P调用exLucas(n, m, P)即可。但是直接使用上面的代码可能会在数据较大时超时我们需要关注几个性能瓶颈和优化点。5.1 性能瓶颈分析factorial_prime函数中的循环计算cycle_prod每次递归调用都会计算一次从1到pk的循环积。如果pk很大比如 10^6 级别且n也很大导致递归深度不浅这个开销是巨大的。递归调用factorial_prime递归本身有一定开销但更主要的是重复计算。factorial_prime(n/p, p, pk)会再次计算cycle_prod。质因数分解对 P 的分解是 O(√P) 的在 P 很大如 10^9时可以接受但也是常数开销。5.2 关键优化策略优化1预处理循环节乘积这是最重要的优化。对于给定的p和pkcycle_prod是一个定值。我们可以在计算combination_mod_prime_power之前先计算并存储它避免在递归中重复计算。// 在 combination_mod_prime_power 函数内部或外部预处理 long long cycle_prod 1; for (long long i 1; i pk; i) { if (i % p ! 0) { cycle_prod cycle_prod * i % pk; } } // 然后将 cycle_prod 作为参数传递给 factorial_prime或者设为全局/静态变量。 // 修改 factorial_prime 函数接收这个预计算好的 cycle_prod。 long long factorial_prime(long long n, long long p, long long pk, long long cycle_prod) { if (n 0) return 1; long long res qpow(cycle_prod, n / pk, pk); // ... 剩余部分不变 }优化2将递归改为迭代factorial_prime的递归形式清晰但可以改为迭代形式效率略高且避免了递归栈溢出的风险虽然此题一般不会。long long factorial_prime_iter(long long n, long long p, long long pk, long long cycle_prod) { long long res 1; while (n 0) { res res * qpow(cycle_prod, n / pk, pk) % pk; for (long long i (n / pk) * pk 1; i n; i) { if (i % p ! 0) { res res * (i % pk) % pk; } } n / p; // 关键对应递归中的 n/p } return res; }这个迭代版本模拟了递归过程每次循环处理当前n的互质部分和剩余部分然后将n更新为n/p直到n为 0。优化3使用更快的质因数分解对于巨大的 P可以使用 Pollard-Rho 算法进行质因数分解但这超出了模板题的一般范围。P4720 的数据范围下试除法足够。5.3 一个优化后的整合示例核心部分结合优化combination_mod_prime_power函数可以这样写long long combination_mod_prime_power(long long n, long long m, long long p, long long pk) { if (m n) return 0; // 预处理循环节乘积 long long cycle_prod 1; for (long long i 1; i pk; i) { if (i % p ! 0) cycle_prod cycle_prod * i % pk; } auto factorial [](long long x) - long long { long long res 1; long long tx x; while (tx 0) { res res * qpow(cycle_prod, tx / pk, pk) % pk; for (long long i (tx / pk) * pk 1; i tx; i) { if (i % p ! 0) res res * (i % pk) % pk; } tx / p; } return res; }; long long e count_prime_factor(n, p) - count_prime_factor(m, p) - count_prime_factor(n - m, p); // 如果 e 已经大于等于 k (即 pk 中 p 的幂次)可以提前返回0。 // 计算 k: 通过不断除以 p 得到 long long temp_pk pk, k 0; while (temp_pk % p 0) { temp_pk / p; k; } if (e k) return 0; long long fn factorial(n); long long fm factorial(m); long long fnm factorial(n - m); long long num fn * inv(fm, pk) % pk * inv(fnm, pk) % pk; long long pe qpow(p, e, pk); return num * pe % pk; }6. 边界条件与常见错误排查即使理解了原理和代码在实际编码和调试中依然会遇到一些隐蔽的坑。6.1 数据范围与溢出处理这是数论题最经典的坑。题目中 n, m 可能高达 10^18P 在 10^6 以内。qpow中的乘法溢出res * base或base * base可能超过long long范围约 9e18。即使对mod取模在乘法运算前就可能溢出了。必须使用快速乘龟速乘或__int128。// 使用 __int128 的快速幂推荐前提是编译器支持 long long qpow(long long base, long long exp, long long mod) { __int128 res 1 % mod; __int128 b base % mod; while (exp) { if (exp 1) res (res * b) % mod; b (b * b) % mod; exp 1; } return (long long)res; } // 或者在无法使用 __int128 时使用快速乘 long long mul_mod(long long a, long long b, long long mod) { long long res 0; a % mod; b % mod; while (b) { if (b 1) res (res a) % mod; a (a a) % mod; b 1; } return res; }递归/迭代中的变量范围在factorial_prime的循环for (long long i (n / pk) * pk 1; i n; i)中(n / pk) * pk的计算可能溢出。更安全的写法是long long start (n / pk) * pk; if (start 0) start 0;或者直接利用取模性质。6.2 特殊模数处理模数 P 1根据定义任何数模 1 都为 0。在代码开头应特判。质数幂 pk 可能为 1在质因数分解时如果 p^k 1这没有意义。实际上当 P 分解时pk 至少为 p。但若 P1已特判。CRT 合并中的模数在增量法合并同余方程时新的模数是m1 / g * m2即lcm(m1, m2)。务必注意计算顺序先除后乘避免中间结果溢出。可以使用__int128辅助计算。6.3 调试技巧当结果错误时可以按以下步骤隔离问题测试小数据用小的 n, m 和小的合数 P如 P6, 10进行测试与暴力计算的结果对比。分离测试子函数测试count_prime_factor验证 n! 中因子 p 的个数计算是否正确。测试factorial_prime选择小的 n, p, pk手动计算验证。测试combination_mod_prime_power针对单个质数幂模数测试。最后测试exLucas和 CRT 合并。验证质因数分解确保对 P 的分解是正确的。检查逆元计算确保在求逆元时inv函数传入的参数与模数互质。在combination_mod_prime_power中求逆元前可以加断言assert(gcd(fm, pk)1 gcd(fnm, pk)1)调试时。7. 扩展卢卡斯的其他应用与变体掌握了这个模板你不仅能解决 P4720还能处理一系列衍生问题。7.1 计算大组合数模任意数这是最直接的应用。在一些计数问题中模数可能不是质数比如998244353 * 1000000007这种扩展卢卡斯是唯一通用的方法。7.2 处理模数较小但n,m巨大的情况即使模数 P 是质数如果 n 和 m 巨大远超 P卢卡斯定理可以将问题规模缩小到 P 以内。而如果 P 不是质数扩展卢卡斯是唯一选择。它通过递归将 n! 分解巧妙地处理了 n 很大的情况。7.3 与多项式、生成函数结合在一些更复杂的组合恒等式证明或求和问题中需要处理模任意数的二项式系数。扩展卢卡斯提供的C(n, m) mod P的能力可以作为子程序嵌入到更大的算法框架中。7.4 局限性扩展卢卡斯定理的时间复杂度主要取决于模数 P 的质因数分解和每个质数幂的大小。设 P 分解为 ∏ pi^ki则时间复杂度约为 O(∑ (ki * pi log n))。当 P 包含大的质数幂如 2^30时计算cycle_prod的循环会非常慢。在这种情况下算法可能不再适用需要寻找其他数学方法或题目给定的特殊约束。在我自己的使用经验里扩展卢卡斯是一个“知道原理就能写但想写对、写快需要很多细节打磨”的算法。它不像快速幂或欧拉筛那样有几乎固定的短代码它的实现长度和细节处理恰恰体现了数论算法从理论到实践的复杂性。理解F(n, p, pk)那个递归式是第一步而处理好循环节、溢出、CRT合并以及各种边界条件才是能在比赛中稳定拿分的关键。建议在理解的基础上亲手实现并通过 P4720 这道模板题进行测试过程中遇到的每一个错误都会让你对模运算和递归分解有更深的认识。