1. 项目概述为什么我们需要“扩展”卢卡斯定理在组合数学和算法竞赛的领域里计算大组合数取模是一个高频且棘手的问题。经典的卢卡斯定理Lucas Theorem为我们提供了一个优雅的解决方案但它有一个非常严格的限制模数p必须是一个质数。这个限制在现实问题中尤其是在密码学、数论和某些特定算法场景下显得过于苛刻。很多时候我们遇到的模数并非质数而是一个合数比如p 10007质数但更多时候可能是p 1000000007质数的变体或者干脆是p 1000、p 12这样的合数。当模数p是合数时经典的卢卡斯定理就失效了。这时“扩展卢卡斯定理”ExLucas便应运而生。它并不是一个像“扩展欧几里得”那样对原定理的直接推广而是一套系统性的算法思想用于解决模数为任意正整数特别是含有质数幂因子时组合数取模的计算问题。简单来说ExLucas 的核心思想是将合数模数分解质因数分别计算组合数对每个质数幂取模的结果最后用中国剩余定理CRT合并答案。我最初接触 ExLucas 是在解决一个涉及非质数模数的计数问题时当时卡了很久翻遍资料才发现这个“利器”。它不像快速幂那样直观也不像动态规划那样有套路更像是一个精巧的“工具箱”你需要理解其中每个工具质因数分解、逆元、阶乘处理、CRT的原理和组装方式。掌握了它你就能攻克一大类模数为合数的组合数难题。2. 核心思路拆解化整为零分而治之ExLucas 的算法流程可以清晰地分为几个步骤理解这个流程是掌握它的关键。我们先从宏观上看再深入每个细节。2.1 算法总览四步走战略假设我们要计算C(n, m) mod p其中p是任意正整数。质因数分解将模数p分解为若干个质数幂的乘积即p p1^k1 * p2^k2 * ... * pt^kt。例如p12分解为2^2 * 3^1。子问题求解对于每一个质数幂因子pi^ki单独计算C(n, m) mod (pi^ki)。这是整个算法最核心、最复杂的一步。因为pi^ki不是质数我们不能直接使用费马小定理求逆元需要特殊处理。中国剩余定理合并我们得到了t个同余方程C(n, m) ≡ a1 (mod p1^k1)C(n, m) ≡ a2 (mod p2^k2)...C(n, m) ≡ at (mod pt^kt)利用中国剩余定理CRT我们可以求出唯一满足所有方程的解x mod p这个x就是我们要的最终答案。结果整合将 CRT 计算出的x返回。整个思路体现了经典的“分治”思想将一个复杂的大问题模合数分解为若干个相对简单的子问题模质数幂分别解决后再合并。2.2 为什么不能直接用逆元这是新手最容易困惑的地方。在模质数p下我们可以用费马小定理a^(p-1) ≡ 1 (mod p)求出a的逆元a^(p-2)从而轻松计算n! / (m! (n-m)!) mod p。 但在模pi^ki下pi是质数pi^ki却不是质数。关键问题在于分母m!和(n-m)!可能与模数pi^ki不互质因为它们可能含有因子pi。而不互质的数在模运算下没有乘法逆元这意味着我们无法直接进行除法运算。注意这是 ExLucas 与普通 Lucas 最根本的区别。Lucas 定理成立的前提是模数为质数确保了分母与模数互质除非分母是模数的倍数但那种情况组合数为0。ExLucas 必须解决这个“不互质”的除法问题。3. 核心攻坚如何计算 C(n, m) mod (p^k)这是 ExLucas 算法的灵魂。设我们要计算C(n, m) mod (p^k)其中p是质数k是正整数。 组合数公式为C(n, m) n! / (m! * (n-m)!)由于分母可能与p^k不互质我们不能直接求逆元。解决思路是将阶乘中所有与p有关的因子即p的幂次和与p无关的部分分离开来。3.1 阶乘的“抽离”表示法我们定义函数F(n, p, k)用于计算n!在模p^k意义下的值但同时要将所有因子p提取出来。 具体地我们将n!写成如下形式n! p^e * f其中e是n!中质因子p的个数。f是n!去掉所有因子p之后剩下的部分对p^k取模的结果并且f与p互质。例如计算10! mod (2^3)。10! 3628800。 首先10!中因子2的个数e可以通过勒让德公式计算e floor(10/2)floor(10/4)floor(10/8)5218。 然后我们计算f即10! / 2^8对8取模。10! / 2^8 3628800 / 256 14175。14175 mod 8 7。 所以10!可以表示为2^8 * 7 (mod 8)。这里f7与2互质。3.2 递归计算 F(n, p, k)如何高效地计算这个f即去p因子后的乘积模p^k呢观察n!的构成n! 1 * 2 * 3 * ... * n我们可以把这些数按是否包含因子p来分组所有是p的倍数的数p, 2p, 3p, ..., floor(n/p)*p。将它们都提取一个p出来这部分贡献了p^(floor(n/p))和一个新的阶乘floor(n/p)!。所有不是p的倍数的数这些数在模p^k下会形成循环节。因为p^k通常不大我们可以预处理出一个周期内的乘积。由此得到递归公式F(n, p, k) [ (乘积周期) * F(floor(n/p), p, k) ] mod p^k具体计算步骤预处理周期乘积计算数组pre[i] (1 * 2 * ... * i) mod p^k其中i取1到p^k但跳过所有p的倍数。实际上我们只需要计算一个完整周期[1, p^k)内与p互质的数的乘积模p^k记为P。因为(a p^k) ≡ a (mod p^k)所以周期是p^k。递归计算n的完整周期个数cnt n / (p^k)最后一个不完整周期的长度rem n % (p^k)周期部分的贡献P^cnt mod p^k用快速幂计算。不完整周期部分的贡献pre[rem]即前rem个与p互质的数的乘积。递归部分F(n/p, p, k)因为提取了所有p的倍数中的因子p后剩下了一个(n/p)!需要处理。最终f (周期部分 * 不完整周期部分 * F(n/p, p, k)) mod p^k。计算指数 en!中因子p的个数e可以用勒让德公式直接计算e n/p n/p^2 n/p^3 ...直到p^i n。3.3 组合数的计算现在我们有了计算n! p^{e_n} * f_n mod p^k的能力。同理可以计算m! p^{e_m} * f_m和(n-m)! p^{e_{nm}} * f_{nm}。 那么C(n, m) n! / (m! * (n-m)!) p^{e_n - e_m - e_{nm}} * (f_n * inv(f_m) * inv(f_{nm}) mod p^k)这里inv(x)表示x在模p^k意义下的逆元。由于f_n, f_m, f_{nm}都与p互质所以它们在模p^k下确实存在逆元可以用扩展欧几里得算法求出。最终结果的处理如果e_n - e_m - e_{nm} k说明组合数分子中p的幂次比分母多至少k次那么整个组合数能被p^k整除所以C(n, m) mod p^k 0。否则结果为p^{e_n - e_m - e_{nm}} * (f_n * inv(f_m) * inv(f_{nm}) mod p^k) mod p^k。4. 完整算法实现与代码解析理解了原理我们来看一个完整的、适用于算法竞赛的实现。代码会包含质因数分解、计算F(n, p, k)、求逆元、快速幂、中国剩余定理等部分。4.1 辅助函数准备首先我们需要一些基础数论函数// 快速幂 (a^b % mod) long long qpow(long long a, long long b, long long mod) { long long res 1; while (b) { if (b 1) res res * a % mod; a a * a % mod; b 1; } return res; } // 扩展欧几里得算法求解 ax by gcd(a, b)返回 gcd(a,b) // 同时通过引用返回 x, y使得 ax by gcd(a,b) 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); return (x % mod mod) % mod; // 调整到 0~mod-1 范围内 }4.2 核心函数计算 n! 中去除 p 因子后的部分 mod p^k这个函数对应上文中的F(n, p, k)计算f部分。// 计算 n! 中剔除所有因子 p 后模 p^k 的结果 // 返回一个 pairlong long, long long first 是 f (乘积部分)second 是 e (p的指数) pairlong long, long long factorial_prime_power(long long n, long long p, long long pk) { if (n 0) return {1, 0}; // 递归计算n! (1*2*...*(p-1) * (p1)*...*(2p-1) * ... ) * p^(n/p) * (n/p)! // 周期部分每 pk 个数一个周期每个周期内与p互质的数的乘积相同 long long period_prod 1; // 一个周期的乘积 long long rem_prod 1; // 最后一个不完整周期的乘积 // 预处理一个周期的乘积跳过p的倍数 // 这里为了清晰直接循环计算。当 pk 较大时可以预处理。 // 注意实际竞赛中pk即 p^k通常不会太大否则计算量爆炸所以直接循环是可接受的。 for (long long i 1; i pk; i) { if (i % p ! 0) { period_prod period_prod * i % pk; } } long long cnt n / pk; // 完整周期个数 long long rem n % pk; // 最后一个不完整周期的长度 // 计算完整周期部分的贡献period_prod ^ cnt long long cycle_part qpow(period_prod, cnt, pk); // 计算不完整周期部分的贡献 for (long long i 1; i rem; i) { if (i % p ! 0) { rem_prod rem_prod * i % pk; } } // 递归计算 (n/p)! 的部分 pairlong long, long long rec factorial_prime_power(n / p, p, pk); // 合并结果f (cycle_part * rem_prod * rec.first) % pk long long f cycle_part * rem_prod % pk; f f * rec.first % pk; // 计算总指数 e floor(n/p) rec.second long long e n / p rec.second; return {f, e}; }4.3 计算 C(n, m) mod (p^k)// 计算 C(n, m) mod (p^k) p是质数 long long C_mod_prime_power(long long n, long long m, long long p, long long k) { if (m 0 || m n) return 0; long long pk 1; for (int i 0; i k; i) pk * p; // 计算 p^k注意可能溢出实际情况中p^k通常不大 // 计算 n!, m!, (n-m)! 的分解形式 auto fn factorial_prime_power(n, p, pk); auto fm factorial_prime_power(m, p, pk); auto fnm factorial_prime_power(n - m, p, pk); // 提取指数和乘积部分 long long e fn.second - fm.second - fnm.second; long long f_n fn.first; long long f_m fm.first; long long f_nm fnm.first; // 如果组合数中p的因子数量 k则结果为0 if (e k) return 0; // 计算 f_n * inv(f_m) * inv(f_nm) mod pk long long comb_f f_n * inv(f_m, pk) % pk; comb_f comb_f * inv(f_nm, pk) % pk; // 最终结果 p^e * comb_f mod pk long long res qpow(p, e, pk) * comb_f % pk; return res; }4.4 中国剩余定理 (CRT) 合并当我们对每个质数幂因子p_i^{k_i}都求出了a_i C(n, m) mod p_i^{k_i}后需要用 CRT 合并。// 中国剩余定理 (CRT) 合并一组同余方程 // equations: 向量每个元素是 pair(余数, 模数) // 所有模数必须两两互质 long long crt(const vectorpairlong long, long long equations) { long long M 1; // 所有模数的乘积 for (auto eq : equations) { M * eq.second; } long long res 0; for (auto eq : equations) { long long a eq.first; // 余数 long long m eq.second; // 模数 long long Mi M / m; // 除了当前模数外的乘积 long long inv_Mi inv(Mi, m); // Mi 在模 m 下的逆元 res (res a * Mi % M * inv_Mi % M) % M; } return res; }4.5 扩展卢卡斯主函数最后将所有部分组合起来。// 扩展卢卡斯定理主函数计算 C(n, m) mod p long long exlucas(long long n, long long m, long long p) { if (m 0 || m n) return 0; // 1. 质因数分解 p vectorpairlong long, int factors; // 存储 (质因子, 幂次) long long temp p; for (long long i 2; i * i temp; i) { if (temp % i 0) { int cnt 0; while (temp % i 0) { temp / i; cnt; } factors.emplace_back(i, cnt); } } if (temp 1) { factors.emplace_back(temp, 1); } // 2. 对每个质数幂因子分别计算 vectorpairlong long, long long equations; // 用于CRT的方程组 for (auto fac : factors) { long long prime fac.first; int power fac.second; long long pk 1; for (int i 0; i power; i) pk * prime; long long a C_mod_prime_power(n, m, prime, power); // C(n,m) mod p^k equations.emplace_back(a, pk); } // 3. 用中国剩余定理合并结果 if (equations.empty()) { // p1 的特殊情况 return 1 % p; // 实际上 p1 时任何数 mod 1 都是 0 } long long ans crt(equations); return ans; }5. 实战应用、常见问题与避坑指南理论很丰满现实很骨感。在实际编码和解题中你会遇到各种各样的问题。下面是我在多次使用 ExLucas 后总结的一些经验和坑点。5.1 典型应用场景算法竞赛这是 ExLucas 最主要的舞台。题目通常会明确给出一个合数模数p比如p10007是质数但有时就是p1000要求计算大组合数取模。直接套用上述模板即可。非质数模数下的计数问题在一些需要取模的计数 DP 或排列组合问题中如果模数不是质数且涉及除法求组合数ExLucas 是唯一选择。小范围验证有时在推导数学公式时可以用 ExLucas 对小模数进行暴力计算来验证公式的正确性。5.2 常见问题与排查技巧问题1时间复杂度与性能瓶颈ExLucas 的复杂度主要取决于模数p的质因数分解和每个质数幂p^k的大小。递归深度factorial_prime_power函数递归调用n/p深度为O(log_p n)。对于n高达1e18p很小的情况递归深度会很大但仍在可接受范围。周期计算在factorial_prime_power中我们循环pk次来计算周期积。这是最大的性能瓶颈如果p^k很大比如超过1e6这个循环将非常慢甚至不可行。优化预处理在调用exlucas前针对每个不同的(p, pk)预先计算好period_prod一个周期的乘积。因为对于同一组(p, pk)这个值是不变的。可以将它缓存起来避免重复计算。问题2数值溢出这是 C 实现中极易出错的地方。pk p^k的计算可能溢出long long。在竞赛中p和k通常不会太大使得p^k在1e18以内。但安全起见可以在循环中判断if (pk LLONG_MAX / p) break;或直接用__int128暂存。在qpow,factorial_prime_power的乘法运算中两个long long相乘可能溢出即使最后会取模。需要使用快速乘龟速乘或__int128来避免中间结果溢出。// 快速乘 (计算 a*b % mod防止溢出) long long qmul(long long a, long long b, long long mod) { long long res 0; while (b) { if (b 1) res (res a) % mod; a (a a) % mod; b 1; } return res; } // 然后在快速幂和乘法中用 qmul 代替直接的 *。问题3特殊边界条件m n或m 0组合数定义为 0需要在入口函数exlucas和C_mod_prime_power中处理。p 1模 1 结果恒为 0。我们的质因数分解循环会跳过导致equations为空。需要在 CRT 合并前判断直接返回0。n或m为 0C(n, 0) 1C(0, m) (m0) 0。递归函数factorial_prime_power(0, ...)能正确处理。问题4中国剩余定理CRT的适用条件我们实现的 CRT 要求所有模数p_i^{k_i}两两互质。这由p的质因数分解保证不同的质数幂必然互质。所以这个条件是天然满足的无需额外检查。5.3 调试与测试技巧当你觉得 ExLucas 代码写对了但结果不对时可以按以下步骤排查小数据暴力验证写一个暴力计算组合数不取模或用 Python 大整数和取模的脚本对小范围的n, m, p比如n,m20, p100进行对拍。这是最有效的调试方法。分步输出在factorial_prime_power、C_mod_prime_power函数中打印中间变量n, p, pk, f, e与手算结果对比。检查逆元单独测试inv函数确保在模pk下计算正确。检查质因数分解确保p被正确分解。检查 CRT 合并手动构造几个简单的同余方程组测试crt函数是否正确。5.4 一个完整的测试案例假设我们要计算C(100, 50) mod 121。质因数分解121 11^2。所以只有一个因子p11, k2, pk121。计算C(100, 50) mod 121。调用factorial_prime_power(100, 11, 121)计算100!。调用factorial_prime_power(50, 11, 121)计算50!。调用factorial_prime_power(50, 11, 121)计算(100-50)!。计算指数e和乘积部分f最终得到结果a1。方程组只有一个x ≡ a1 (mod 121)。CRT 合并后结果就是a1。你可以用 Python 验证math.comb(100, 50) % 121的结果应该与你的程序输出一致。6. 性能优化与扩展思考基础的 ExLucas 实现已经能解决大部分问题但在极端数据下可能需要优化。6.1 关键优化预处理周期乘积如前所述对于固定的(p, pk)period_prod是常量。我们可以在主函数exlucas中对每个分解出的(p, k)预先计算并存储period_prod[pk]。这样在递归调用factorial_prime_power时可以直接查表避免每次递归都进行O(pk)的循环。unordered_maplong long, long long period_prod_cache; long long get_period_prod(long long p, long long pk) { if (period_prod_cache.count(pk)) return period_prod_cache[pk]; long long res 1; for (long long i 1; i pk; i) { if (i % p ! 0) { res res * i % pk; } } period_prod_cache[pk] res; return res; } // 然后在 factorial_prime_power 函数中用 get_period_prod(p, pk) 代替循环计算 period_prod。6.2 应对大模数 p^k如果p^k非常大比如p很小但k很大使得p^k接近1e9预处理循环O(p^k)是不可接受的。此时需要更高级的优化例如利用威尔逊定理的推广和递归分治来快速计算一个周期内与p互质的数的乘积模p^k。但这已超出一般竞赛范围属于数论深水区。对于竞赛题目设计通常会保证p^k在一个可计算的范围内如p^k 1e6。6.3 与普通卢卡斯定理的选择如果模数p是质数一定要用普通的卢卡斯定理而不是 ExLucas。因为普通 Lucas 的时间复杂度是O(log_p n)且代码简单效率远高于 ExLucas。所以在解决问题时第一步应该是判断p是否为质数可以用 Miller-Rabin 素性测试。如果是质数就用 Lucas否则再用 ExLucas。bool is_prime(long long p); // 素性测试函数 long long lucas(long long n, long long m, long long p); // 普通卢卡斯定理实现 long long solve_combination(long long n, long long m, long long p) { if (is_prime(p)) { return lucas(n, m, p); } else { return exlucas(n, m, p); } }6.4 模数非质数但 n, m 较小时当n和m不大比如小于1000时即使模数是合数也有更简单的办法直接预处理组合数表计算时对合数模数取模。因为组合数计算过程中只涉及乘法和除法我们可以用扩展欧几里得算法对每个除法操作求逆元。但注意这要求每次除法时分母与当前模数p互质。由于p是合数这可能不总是成立。因此这种方法仅适用于你能确保运算过程中所有分母都与p互质的场景或者p的质因子很大n,m很小以至于不会触及这些因子。7. 总结与个人体会扩展卢卡斯定理是一个将多个基础数论知识质因数分解、阶乘性质、逆元、中国剩余定理综合应用的典型例子。它没有引入新的高深定理而是通过巧妙的组合解决了“模合数下求组合数”这一实际问题。我在实际使用中最大的体会是一定要理解其“分治”本质。不要试图去记忆那一长串递归公式而是记住“分解模数 - 解决质数幂子问题 - CRT合并”这个核心流程。对于子问题记住关键点是“把因子p抽离出来使得剩余部分与p互质从而可以求逆元”。另一个深刻的教训是关于溢出处理。在第一次实现时我因为long long乘法溢出调试了整整一个下午。后来养成了习惯在涉及大数乘法的模运算中要么用__int128要么就老老实实写上快速乘函数。最后ExLucas 在竞赛中属于“核武器”级别的算法代码较长运行时间也相对较慢。因此在比赛中要谨慎使用确保题目真的需要它模数是合数且n,m很大。如果模数是质数或者n,m很小都有更优的解法。把它作为你数论工具箱里的一件重型装备在合适的时机亮出来才能一击制胜。 扩展卢卡斯定理：解决模合数下大组合数计算的算法精解
1. 项目概述为什么我们需要“扩展”卢卡斯定理在组合数学和算法竞赛的领域里计算大组合数取模是一个高频且棘手的问题。经典的卢卡斯定理Lucas Theorem为我们提供了一个优雅的解决方案但它有一个非常严格的限制模数p必须是一个质数。这个限制在现实问题中尤其是在密码学、数论和某些特定算法场景下显得过于苛刻。很多时候我们遇到的模数并非质数而是一个合数比如p 10007质数但更多时候可能是p 1000000007质数的变体或者干脆是p 1000、p 12这样的合数。当模数p是合数时经典的卢卡斯定理就失效了。这时“扩展卢卡斯定理”ExLucas便应运而生。它并不是一个像“扩展欧几里得”那样对原定理的直接推广而是一套系统性的算法思想用于解决模数为任意正整数特别是含有质数幂因子时组合数取模的计算问题。简单来说ExLucas 的核心思想是将合数模数分解质因数分别计算组合数对每个质数幂取模的结果最后用中国剩余定理CRT合并答案。我最初接触 ExLucas 是在解决一个涉及非质数模数的计数问题时当时卡了很久翻遍资料才发现这个“利器”。它不像快速幂那样直观也不像动态规划那样有套路更像是一个精巧的“工具箱”你需要理解其中每个工具质因数分解、逆元、阶乘处理、CRT的原理和组装方式。掌握了它你就能攻克一大类模数为合数的组合数难题。2. 核心思路拆解化整为零分而治之ExLucas 的算法流程可以清晰地分为几个步骤理解这个流程是掌握它的关键。我们先从宏观上看再深入每个细节。2.1 算法总览四步走战略假设我们要计算C(n, m) mod p其中p是任意正整数。质因数分解将模数p分解为若干个质数幂的乘积即p p1^k1 * p2^k2 * ... * pt^kt。例如p12分解为2^2 * 3^1。子问题求解对于每一个质数幂因子pi^ki单独计算C(n, m) mod (pi^ki)。这是整个算法最核心、最复杂的一步。因为pi^ki不是质数我们不能直接使用费马小定理求逆元需要特殊处理。中国剩余定理合并我们得到了t个同余方程C(n, m) ≡ a1 (mod p1^k1)C(n, m) ≡ a2 (mod p2^k2)...C(n, m) ≡ at (mod pt^kt)利用中国剩余定理CRT我们可以求出唯一满足所有方程的解x mod p这个x就是我们要的最终答案。结果整合将 CRT 计算出的x返回。整个思路体现了经典的“分治”思想将一个复杂的大问题模合数分解为若干个相对简单的子问题模质数幂分别解决后再合并。2.2 为什么不能直接用逆元这是新手最容易困惑的地方。在模质数p下我们可以用费马小定理a^(p-1) ≡ 1 (mod p)求出a的逆元a^(p-2)从而轻松计算n! / (m! (n-m)!) mod p。 但在模pi^ki下pi是质数pi^ki却不是质数。关键问题在于分母m!和(n-m)!可能与模数pi^ki不互质因为它们可能含有因子pi。而不互质的数在模运算下没有乘法逆元这意味着我们无法直接进行除法运算。注意这是 ExLucas 与普通 Lucas 最根本的区别。Lucas 定理成立的前提是模数为质数确保了分母与模数互质除非分母是模数的倍数但那种情况组合数为0。ExLucas 必须解决这个“不互质”的除法问题。3. 核心攻坚如何计算 C(n, m) mod (p^k)这是 ExLucas 算法的灵魂。设我们要计算C(n, m) mod (p^k)其中p是质数k是正整数。 组合数公式为C(n, m) n! / (m! * (n-m)!)由于分母可能与p^k不互质我们不能直接求逆元。解决思路是将阶乘中所有与p有关的因子即p的幂次和与p无关的部分分离开来。3.1 阶乘的“抽离”表示法我们定义函数F(n, p, k)用于计算n!在模p^k意义下的值但同时要将所有因子p提取出来。 具体地我们将n!写成如下形式n! p^e * f其中e是n!中质因子p的个数。f是n!去掉所有因子p之后剩下的部分对p^k取模的结果并且f与p互质。例如计算10! mod (2^3)。10! 3628800。 首先10!中因子2的个数e可以通过勒让德公式计算e floor(10/2)floor(10/4)floor(10/8)5218。 然后我们计算f即10! / 2^8对8取模。10! / 2^8 3628800 / 256 14175。14175 mod 8 7。 所以10!可以表示为2^8 * 7 (mod 8)。这里f7与2互质。3.2 递归计算 F(n, p, k)如何高效地计算这个f即去p因子后的乘积模p^k呢观察n!的构成n! 1 * 2 * 3 * ... * n我们可以把这些数按是否包含因子p来分组所有是p的倍数的数p, 2p, 3p, ..., floor(n/p)*p。将它们都提取一个p出来这部分贡献了p^(floor(n/p))和一个新的阶乘floor(n/p)!。所有不是p的倍数的数这些数在模p^k下会形成循环节。因为p^k通常不大我们可以预处理出一个周期内的乘积。由此得到递归公式F(n, p, k) [ (乘积周期) * F(floor(n/p), p, k) ] mod p^k具体计算步骤预处理周期乘积计算数组pre[i] (1 * 2 * ... * i) mod p^k其中i取1到p^k但跳过所有p的倍数。实际上我们只需要计算一个完整周期[1, p^k)内与p互质的数的乘积模p^k记为P。因为(a p^k) ≡ a (mod p^k)所以周期是p^k。递归计算n的完整周期个数cnt n / (p^k)最后一个不完整周期的长度rem n % (p^k)周期部分的贡献P^cnt mod p^k用快速幂计算。不完整周期部分的贡献pre[rem]即前rem个与p互质的数的乘积。递归部分F(n/p, p, k)因为提取了所有p的倍数中的因子p后剩下了一个(n/p)!需要处理。最终f (周期部分 * 不完整周期部分 * F(n/p, p, k)) mod p^k。计算指数 en!中因子p的个数e可以用勒让德公式直接计算e n/p n/p^2 n/p^3 ...直到p^i n。3.3 组合数的计算现在我们有了计算n! p^{e_n} * f_n mod p^k的能力。同理可以计算m! p^{e_m} * f_m和(n-m)! p^{e_{nm}} * f_{nm}。 那么C(n, m) n! / (m! * (n-m)!) p^{e_n - e_m - e_{nm}} * (f_n * inv(f_m) * inv(f_{nm}) mod p^k)这里inv(x)表示x在模p^k意义下的逆元。由于f_n, f_m, f_{nm}都与p互质所以它们在模p^k下确实存在逆元可以用扩展欧几里得算法求出。最终结果的处理如果e_n - e_m - e_{nm} k说明组合数分子中p的幂次比分母多至少k次那么整个组合数能被p^k整除所以C(n, m) mod p^k 0。否则结果为p^{e_n - e_m - e_{nm}} * (f_n * inv(f_m) * inv(f_{nm}) mod p^k) mod p^k。4. 完整算法实现与代码解析理解了原理我们来看一个完整的、适用于算法竞赛的实现。代码会包含质因数分解、计算F(n, p, k)、求逆元、快速幂、中国剩余定理等部分。4.1 辅助函数准备首先我们需要一些基础数论函数// 快速幂 (a^b % mod) long long qpow(long long a, long long b, long long mod) { long long res 1; while (b) { if (b 1) res res * a % mod; a a * a % mod; b 1; } return res; } // 扩展欧几里得算法求解 ax by gcd(a, b)返回 gcd(a,b) // 同时通过引用返回 x, y使得 ax by gcd(a,b) 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); return (x % mod mod) % mod; // 调整到 0~mod-1 范围内 }4.2 核心函数计算 n! 中去除 p 因子后的部分 mod p^k这个函数对应上文中的F(n, p, k)计算f部分。// 计算 n! 中剔除所有因子 p 后模 p^k 的结果 // 返回一个 pairlong long, long long first 是 f (乘积部分)second 是 e (p的指数) pairlong long, long long factorial_prime_power(long long n, long long p, long long pk) { if (n 0) return {1, 0}; // 递归计算n! (1*2*...*(p-1) * (p1)*...*(2p-1) * ... ) * p^(n/p) * (n/p)! // 周期部分每 pk 个数一个周期每个周期内与p互质的数的乘积相同 long long period_prod 1; // 一个周期的乘积 long long rem_prod 1; // 最后一个不完整周期的乘积 // 预处理一个周期的乘积跳过p的倍数 // 这里为了清晰直接循环计算。当 pk 较大时可以预处理。 // 注意实际竞赛中pk即 p^k通常不会太大否则计算量爆炸所以直接循环是可接受的。 for (long long i 1; i pk; i) { if (i % p ! 0) { period_prod period_prod * i % pk; } } long long cnt n / pk; // 完整周期个数 long long rem n % pk; // 最后一个不完整周期的长度 // 计算完整周期部分的贡献period_prod ^ cnt long long cycle_part qpow(period_prod, cnt, pk); // 计算不完整周期部分的贡献 for (long long i 1; i rem; i) { if (i % p ! 0) { rem_prod rem_prod * i % pk; } } // 递归计算 (n/p)! 的部分 pairlong long, long long rec factorial_prime_power(n / p, p, pk); // 合并结果f (cycle_part * rem_prod * rec.first) % pk long long f cycle_part * rem_prod % pk; f f * rec.first % pk; // 计算总指数 e floor(n/p) rec.second long long e n / p rec.second; return {f, e}; }4.3 计算 C(n, m) mod (p^k)// 计算 C(n, m) mod (p^k) p是质数 long long C_mod_prime_power(long long n, long long m, long long p, long long k) { if (m 0 || m n) return 0; long long pk 1; for (int i 0; i k; i) pk * p; // 计算 p^k注意可能溢出实际情况中p^k通常不大 // 计算 n!, m!, (n-m)! 的分解形式 auto fn factorial_prime_power(n, p, pk); auto fm factorial_prime_power(m, p, pk); auto fnm factorial_prime_power(n - m, p, pk); // 提取指数和乘积部分 long long e fn.second - fm.second - fnm.second; long long f_n fn.first; long long f_m fm.first; long long f_nm fnm.first; // 如果组合数中p的因子数量 k则结果为0 if (e k) return 0; // 计算 f_n * inv(f_m) * inv(f_nm) mod pk long long comb_f f_n * inv(f_m, pk) % pk; comb_f comb_f * inv(f_nm, pk) % pk; // 最终结果 p^e * comb_f mod pk long long res qpow(p, e, pk) * comb_f % pk; return res; }4.4 中国剩余定理 (CRT) 合并当我们对每个质数幂因子p_i^{k_i}都求出了a_i C(n, m) mod p_i^{k_i}后需要用 CRT 合并。// 中国剩余定理 (CRT) 合并一组同余方程 // equations: 向量每个元素是 pair(余数, 模数) // 所有模数必须两两互质 long long crt(const vectorpairlong long, long long equations) { long long M 1; // 所有模数的乘积 for (auto eq : equations) { M * eq.second; } long long res 0; for (auto eq : equations) { long long a eq.first; // 余数 long long m eq.second; // 模数 long long Mi M / m; // 除了当前模数外的乘积 long long inv_Mi inv(Mi, m); // Mi 在模 m 下的逆元 res (res a * Mi % M * inv_Mi % M) % M; } return res; }4.5 扩展卢卡斯主函数最后将所有部分组合起来。// 扩展卢卡斯定理主函数计算 C(n, m) mod p long long exlucas(long long n, long long m, long long p) { if (m 0 || m n) return 0; // 1. 质因数分解 p vectorpairlong long, int factors; // 存储 (质因子, 幂次) long long temp p; for (long long i 2; i * i temp; i) { if (temp % i 0) { int cnt 0; while (temp % i 0) { temp / i; cnt; } factors.emplace_back(i, cnt); } } if (temp 1) { factors.emplace_back(temp, 1); } // 2. 对每个质数幂因子分别计算 vectorpairlong long, long long equations; // 用于CRT的方程组 for (auto fac : factors) { long long prime fac.first; int power fac.second; long long pk 1; for (int i 0; i power; i) pk * prime; long long a C_mod_prime_power(n, m, prime, power); // C(n,m) mod p^k equations.emplace_back(a, pk); } // 3. 用中国剩余定理合并结果 if (equations.empty()) { // p1 的特殊情况 return 1 % p; // 实际上 p1 时任何数 mod 1 都是 0 } long long ans crt(equations); return ans; }5. 实战应用、常见问题与避坑指南理论很丰满现实很骨感。在实际编码和解题中你会遇到各种各样的问题。下面是我在多次使用 ExLucas 后总结的一些经验和坑点。5.1 典型应用场景算法竞赛这是 ExLucas 最主要的舞台。题目通常会明确给出一个合数模数p比如p10007是质数但有时就是p1000要求计算大组合数取模。直接套用上述模板即可。非质数模数下的计数问题在一些需要取模的计数 DP 或排列组合问题中如果模数不是质数且涉及除法求组合数ExLucas 是唯一选择。小范围验证有时在推导数学公式时可以用 ExLucas 对小模数进行暴力计算来验证公式的正确性。5.2 常见问题与排查技巧问题1时间复杂度与性能瓶颈ExLucas 的复杂度主要取决于模数p的质因数分解和每个质数幂p^k的大小。递归深度factorial_prime_power函数递归调用n/p深度为O(log_p n)。对于n高达1e18p很小的情况递归深度会很大但仍在可接受范围。周期计算在factorial_prime_power中我们循环pk次来计算周期积。这是最大的性能瓶颈如果p^k很大比如超过1e6这个循环将非常慢甚至不可行。优化预处理在调用exlucas前针对每个不同的(p, pk)预先计算好period_prod一个周期的乘积。因为对于同一组(p, pk)这个值是不变的。可以将它缓存起来避免重复计算。问题2数值溢出这是 C 实现中极易出错的地方。pk p^k的计算可能溢出long long。在竞赛中p和k通常不会太大使得p^k在1e18以内。但安全起见可以在循环中判断if (pk LLONG_MAX / p) break;或直接用__int128暂存。在qpow,factorial_prime_power的乘法运算中两个long long相乘可能溢出即使最后会取模。需要使用快速乘龟速乘或__int128来避免中间结果溢出。// 快速乘 (计算 a*b % mod防止溢出) long long qmul(long long a, long long b, long long mod) { long long res 0; while (b) { if (b 1) res (res a) % mod; a (a a) % mod; b 1; } return res; } // 然后在快速幂和乘法中用 qmul 代替直接的 *。问题3特殊边界条件m n或m 0组合数定义为 0需要在入口函数exlucas和C_mod_prime_power中处理。p 1模 1 结果恒为 0。我们的质因数分解循环会跳过导致equations为空。需要在 CRT 合并前判断直接返回0。n或m为 0C(n, 0) 1C(0, m) (m0) 0。递归函数factorial_prime_power(0, ...)能正确处理。问题4中国剩余定理CRT的适用条件我们实现的 CRT 要求所有模数p_i^{k_i}两两互质。这由p的质因数分解保证不同的质数幂必然互质。所以这个条件是天然满足的无需额外检查。5.3 调试与测试技巧当你觉得 ExLucas 代码写对了但结果不对时可以按以下步骤排查小数据暴力验证写一个暴力计算组合数不取模或用 Python 大整数和取模的脚本对小范围的n, m, p比如n,m20, p100进行对拍。这是最有效的调试方法。分步输出在factorial_prime_power、C_mod_prime_power函数中打印中间变量n, p, pk, f, e与手算结果对比。检查逆元单独测试inv函数确保在模pk下计算正确。检查质因数分解确保p被正确分解。检查 CRT 合并手动构造几个简单的同余方程组测试crt函数是否正确。5.4 一个完整的测试案例假设我们要计算C(100, 50) mod 121。质因数分解121 11^2。所以只有一个因子p11, k2, pk121。计算C(100, 50) mod 121。调用factorial_prime_power(100, 11, 121)计算100!。调用factorial_prime_power(50, 11, 121)计算50!。调用factorial_prime_power(50, 11, 121)计算(100-50)!。计算指数e和乘积部分f最终得到结果a1。方程组只有一个x ≡ a1 (mod 121)。CRT 合并后结果就是a1。你可以用 Python 验证math.comb(100, 50) % 121的结果应该与你的程序输出一致。6. 性能优化与扩展思考基础的 ExLucas 实现已经能解决大部分问题但在极端数据下可能需要优化。6.1 关键优化预处理周期乘积如前所述对于固定的(p, pk)period_prod是常量。我们可以在主函数exlucas中对每个分解出的(p, k)预先计算并存储period_prod[pk]。这样在递归调用factorial_prime_power时可以直接查表避免每次递归都进行O(pk)的循环。unordered_maplong long, long long period_prod_cache; long long get_period_prod(long long p, long long pk) { if (period_prod_cache.count(pk)) return period_prod_cache[pk]; long long res 1; for (long long i 1; i pk; i) { if (i % p ! 0) { res res * i % pk; } } period_prod_cache[pk] res; return res; } // 然后在 factorial_prime_power 函数中用 get_period_prod(p, pk) 代替循环计算 period_prod。6.2 应对大模数 p^k如果p^k非常大比如p很小但k很大使得p^k接近1e9预处理循环O(p^k)是不可接受的。此时需要更高级的优化例如利用威尔逊定理的推广和递归分治来快速计算一个周期内与p互质的数的乘积模p^k。但这已超出一般竞赛范围属于数论深水区。对于竞赛题目设计通常会保证p^k在一个可计算的范围内如p^k 1e6。6.3 与普通卢卡斯定理的选择如果模数p是质数一定要用普通的卢卡斯定理而不是 ExLucas。因为普通 Lucas 的时间复杂度是O(log_p n)且代码简单效率远高于 ExLucas。所以在解决问题时第一步应该是判断p是否为质数可以用 Miller-Rabin 素性测试。如果是质数就用 Lucas否则再用 ExLucas。bool is_prime(long long p); // 素性测试函数 long long lucas(long long n, long long m, long long p); // 普通卢卡斯定理实现 long long solve_combination(long long n, long long m, long long p) { if (is_prime(p)) { return lucas(n, m, p); } else { return exlucas(n, m, p); } }6.4 模数非质数但 n, m 较小时当n和m不大比如小于1000时即使模数是合数也有更简单的办法直接预处理组合数表计算时对合数模数取模。因为组合数计算过程中只涉及乘法和除法我们可以用扩展欧几里得算法对每个除法操作求逆元。但注意这要求每次除法时分母与当前模数p互质。由于p是合数这可能不总是成立。因此这种方法仅适用于你能确保运算过程中所有分母都与p互质的场景或者p的质因子很大n,m很小以至于不会触及这些因子。7. 总结与个人体会扩展卢卡斯定理是一个将多个基础数论知识质因数分解、阶乘性质、逆元、中国剩余定理综合应用的典型例子。它没有引入新的高深定理而是通过巧妙的组合解决了“模合数下求组合数”这一实际问题。我在实际使用中最大的体会是一定要理解其“分治”本质。不要试图去记忆那一长串递归公式而是记住“分解模数 - 解决质数幂子问题 - CRT合并”这个核心流程。对于子问题记住关键点是“把因子p抽离出来使得剩余部分与p互质从而可以求逆元”。另一个深刻的教训是关于溢出处理。在第一次实现时我因为long long乘法溢出调试了整整一个下午。后来养成了习惯在涉及大数乘法的模运算中要么用__int128要么就老老实实写上快速乘函数。最后ExLucas 在竞赛中属于“核武器”级别的算法代码较长运行时间也相对较慢。因此在比赛中要谨慎使用确保题目真的需要它模数是合数且n,m很大。如果模数是质数或者n,m很小都有更优的解法。把它作为你数论工具箱里的一件重型装备在合适的时机亮出来才能一击制胜。