1. 从组合数取模到Lucas定理一个经典问题的诞生在算法竞赛和数论编程中计算组合数 C(n, m) 取模一个质数 p 的值是一个高频出现的需求。比如计算从 n 个不同元素中选取 m 个的方案数并对一个大质数如 1e97取模。当 n 和 m 的数值较小时我们可以直接使用公式 C(n, m) n! / (m! * (n-m)!) 配合预处理阶乘和乘法逆元来 O(1) 计算。然而当 n 和 m 的数值非常大比如 1e18但模数 p 是一个相对较小的质数比如 1e5时直接计算阶乘会溢出甚至无法存储如此巨大的中间结果。这时我们就需要一种更聪明的工具——Lucas定理。我第一次遇到这个问题是在一个关于“网格路径计数”的题目里n 和 m 都达到了 1e15 级别模数 p 是 10007。直接算内存和时间都不允许。当时我就想肯定存在一种方法能将这个“庞然大物”分解成一系列我们能够处理的小问题。Lucas定理正是这把钥匙。它巧妙地将大数 n 和 m 转换到模数 p 的进制下将问题规模缩小到 p 以内从而使得我们可以用预处理阶乘这种常规手段来解决。理解并掌握 Lucas 定理是迈向解决更复杂模运算问题的重要一步。2. Lucas定理核心原理与严谨证明Lucas定理的核心思想是“降维打击”。它告诉我们当模数 p 是质数时计算 C(n, m) mod p 可以转化为计算一系列更小的组合数的乘积。2.1 定理陈述设 p 是一个质数n 和 m 是非负整数。将 n 和 m 分别用 p 进制表示n n_k * p^k n_{k-1} * p^{k-1} ... n_1 * p n_0m m_k * p^k m_{k-1} * p^{k-1} ... m_1 * p m_0 其中对于所有的 i有 0 ≤ n_i, m_i p。那么Lucas定理断言 C(n, m) ≡ Π_{i0}^{k} C(n_i, m_i) (mod p)简单来说就是把 n 和 m 写成 p 进制然后对应位上的数字分别求组合数最后把所有结果乘起来再对 p 取模。2.2 为什么它是对的一个生成函数的视角最经典的证明利用了二项式系数在模 p 下的性质以及生成函数。这里我分享一个更容易理解其“组合意义”的思路。考虑多项式 (1 x)^n 在模 p 下的展开。根据二项式定理(1 x)^n 的系数就是组合数 C(n, m)。现在我们把 n 写成 p 进制n a * p b其中 a n / p, b n % p。那么 (1 x)^n (1 x)^{a*p b} [(1 x)^p]^a * (1 x)^b这里有一个关键引理费马小定理的一个推论在模 p 下对于质数 p有 (1 x)^p ≡ 1 x^p (mod p)。这是因为除了 C(p, 0) 和 C(p, p) 为 1 外其他 C(p, k) 都包含因子 p在模 p 下为 0。所以上式在模 p 下等价于 (1 x)^n ≡ (1 x^p)^a * (1 x)^b (mod p)现在我们想求左边多项式中 x^m 的系数即 C(n, m) mod p。我们把 m 也写成 p 进制m c * p d。 观察右边(1 x^p)^a 展开后x 的指数都是 p 的倍数 (1 x)^b 展开后x 的指数小于 p。 因此要得到 x^{cp d}我们必须从 (1 x^p)^a 中取 x^{cp} 项其系数为 C(a, c)从 (1 x)^b 中取 x^d 项其系数为 C(b, d)。于是我们得到了 C(n, m) ≡ C(a, c) * C(b, d) (mod p)其中 n ap b, m cp d。这正好是 Lucas 定理在 p 进制下一位的情况。递归地对 a 和 c 应用这个过程即把它们当作新的 n 和 m就得到了完整的定理。这个证明过程清晰地展示了“按位分解”的思想来源。注意这个证明中隐含了一个重要条件即 p 必须是质数否则 (1 x)^p ≡ 1 x^p (mod p) 这个关键等式不成立。这也是经典 Lucas 定理只适用于质数模数的根本原因。2.3 一个具体例子感受降维威力假设我们要计算 C(1234, 789) mod 13。首先13 是质数满足条件。将 1234 和 789 转化为 13 进制。1234 ÷ 13 94 ... 12 (12 对应 C)94 ÷ 13 7 ... 37 ÷ 13 0 ... 7 所以 1234 的 13 进制是 (7, 3, 12)_{13}写作 73C这里用C代表12。789 ÷ 13 60 ... 960 ÷ 13 4 ... 84 ÷ 13 0 ... 4 所以 789 的 13 进制是 (4, 8, 9)_{13}写作 489。根据 Lucas 定理 C(1234, 789) mod 13 ≡ C(7, 4) * C(3, 8) * C(12, 9) mod 13计算每一位的组合数。这里有个细节如果某一位上 m_i n_i根据组合数的定义C(n_i, m_i) 0那么整个乘积就是 0。C(7, 4) 35 35 mod 13 9。C(3, 8)因为 8 3所以结果为 0。由于中间有一项为 0所以最终结果 C(1234, 789) mod 13 ≡ 0。看我们根本没有去计算那巨大的 C(1234, 789)只是算了几个很小的组合数就得到了结果。这就是 Lucas 定理的威力。在实际编程中我们会预处理出所有小于 p 的阶乘 fact[i] 和阶乘的逆元 inv_fact[i]这样计算任意 C(n_i, m_i) (n_i, m_i p) 都是 O(1) 的。3. Lucas定理的代码实现与细节打磨理解了原理实现起来就清晰了。代码的核心是递归或迭代地取出 n 和 m 在 p 进制下的每一位并计算对应的小组合数。3.1 预处理快速计算小组合数的基石由于 Lucas 定理最终归结为计算多个 C(n_i, m_i)其中 n_i, m_i p我们可以预处理出 0 到 p-1 的阶乘模 p 的值及其逆元。// 假设 p 是全局质数 long long fact[MAXP], inv_fact[MAXP]; void init(int p) { fact[0] 1; for (int i 1; i p; i) { fact[i] fact[i-1] * i % p; } // 计算阶乘逆元利用费马小定理 a^(p-2) ≡ a^(-1) (mod p) inv_fact[p-1] fast_pow(fact[p-1], p-2, p); // 快速幂 for (int i p-2; i 0; --i) { inv_fact[i] inv_fact[i1] * (i1) % p; } } // 计算 C(a, b) mod p, 其中 a, b p long long C_small(long long a, long long b, int p) { if (a b) return 0; return fact[a] * inv_fact[b] % p * inv_fact[a - b] % p; }这里fast_pow是标准的快速幂函数。预处理的时间复杂度是 O(p)对于 p 在 1e6 以内通常是可接受的。3.2 Lucas 函数的递归与迭代实现递归实现是最直观的完全对应定理的表述long long lucas(long long n, long long m, int p) { if (m 0) return 1; // C(n, 0) 1 // 递归核心C(n,m) mod p C(n%p, m%p) * lucas(n/p, m/p, p) mod p return C_small(n % p, m % p, p) * lucas(n / p, m / p, p) % p; }递归深度是 n 的 p 进制位数即 O(log_p n)效率很高。迭代实现避免了递归调用栈的开销思路一样用一个循环不断取出最低位long long lucas_iter(long long n, long long m, int p) { long long res 1; while (n 0 || m 0) { // 取出当前最低位 long long ni n % p; long long mi m % p; // 如果某一位上 m_i n_i组合数为0直接返回0 if (ni mi) return 0; // 计算当前位的组合数并乘入结果 res res * C_small(ni, mi, p) % p; // 移除最低位 n / p; m / p; } return res; }我个人更倾向于迭代版本它逻辑清晰且没有递归深度的限制虽然通常也不会超。3.3 实战中的注意事项与性能优化边界条件与快速返回在C_small函数或 Lucas 函数中一旦发现m n应立即返回 0。这是一个有效的剪枝能提前结束计算。模数 p 的范围预处理阶乘需要 O(p) 的时间和 O(p) 的空间。因此经典 Lucas 定理适用于p 是质数且 p 的大小在可接受范围内通常 ≤ 1e6的情况。如果 p 很大但仍是质数预处理会超时或超内存此时需要其他方法如分段计算。组合数计算的防溢出在C_small中我们依靠预处理好的阶乘和逆元进行计算。如果环境不允许预处理比如 p 较大计算 C(n_i, m_i) 时要注意使用long long中间结果并及时取模或者使用更稳健的组合数计算方法。p 必须为质数的验证在调用 Lucas 函数前务必确保传入的 p 是质数。可以写一个简单的素数测试函数如试除法到 sqrt(p)。如果 p 不是质数结果将是错误的。实操心得在竞赛中如果题目明确说明模数是质数如 1e97但 n, m 很大直接套 Lucas 定理是行不通的因为 1e97 太大了无法预处理。这种情况通常 n, m 不会超过模数直接用预处理阶乘即可。Lucas 定理的真正用武之地是模数 p 较小但 n, m 极大的场景。一定要分清应用场景。4. ExLucas当模数不是质数时我们如何破局现实世界和题目中模数常常不是质数。比如模数可能是 1000000000 (1e9)或者是一个合数如 2024。这时经典的 Lucas 定理就失效了。我们需要一个更强大的工具——扩展 Lucas 定理俗称 ExLucas。它的目标是解决C(n, m) mod p其中p 可以是任意正整数。4.1 核心思路中国剩余定理与质因数分解ExLucas 的核心思想是“化繁为简分而治之”。分解将模数 p 分解为质因数幂的形式p p1^k1 * p2^k2 * ... * pt^kt。求解对于每一个质因数幂 pi^ki单独计算 C(n, m) mod pi^ki。这是整个问题的难点所在。合并利用中国剩余定理将 t 个同余方程的解合并得到最终 C(n, m) mod p 的结果。所以问题的关键转化为如何计算 C(n, m) mod p^k其中 p 是质数k 是正整数4.2 攻克核心子问题C(n, m) mod p^k计算 C(n, m) n! / (m! * (n-m)!) mod p^k 的难点在于分母的逆元可能不存在因为分母中的因子可能与 p 不互质。我们不能直接求逆元。解决方法是将阶乘中所有 p 的因子剥离出来。 设函数F(n, p, pk)表示 n! 中除去所有 p 的因子后对 pk (即 p^k) 取模的值。同时设G(n, p)表示 n! 中 p 这个质因子的总次数。那么n! 可以表示为n! p^{G(n, p)} * F(n, p, pk) 于是组合数可以写为 C(n, m) [ p^{G(n,p)} * F(n,p,pk) ] / [ p^{G(m,p)} * F(m,p,pk) * p^{G(n-m,p)} * F(n-m,p,pk) ] p^{ [G(n,p) - G(m,p) - G(n-m,p)] } * [ F(n,p,pk) / ( F(m,p,pk) * F(n-m,p,pk) ) ]在模 pk 的意义下指数部分如果指数 ≥ k则整个式子模 pk 为 0。否则保留 p 的幂次。分数部分现在 F 函数的值是与 p 互质的因为我们把 p 的因子都提走了所以它在模 pk 下存在逆元可以正常计算。因此计算步骤为计算指数 e G(n,p) - G(m,p) - G(n-m,p)。若 e ≥ k则 C(n, m) mod p^k 0。计算 F(n,p,pk), F(m,p,pk), F(n-m,p,pk)。计算分数部分frac F(n,p,pk) * inv(F(m,p,pk), pk) * inv(F(n-m,p,pk), pk) mod pk。这里inv(a, pk)是求 a 在模 pk 下的逆元因为 a 与 p 互质逆元存在。最终结果C(n,m) mod p^k (p^e * frac) mod pk。4.3 实现关键函数 F 与 G计算 G(n, p)n! 中质因子 p 的个数这是一个经典问题公式为G(n, p) floor(n/p) floor(n/p^2) floor(n/p^3) ... 实现起来很简单long long get_g(long long n, int p) { long long cnt 0; while (n) { cnt n / p; n / p; } return cnt; }计算 F(n, p, pk)剥离 p 因子后的 n! mod pk这是 ExLucas 中最精妙的部分。我们不能直接计算 n! 再剥离因为 n 可能很大。观察 n! n! 1 * 2 * 3 * ... * n 我们可以把它按模 p 的余数分组所有 p 的倍数p, 2p, 3p, ..., floor(n/p)*p。这部分每个数都至少含有一个 p 因子。所有不是 p 的倍数的数。因此n! [ 1 * 2 * ... * (p-1) ] * [ (p1) * (p2) * ... * (2p-1) ] * ... * (最后一段非 p 倍数的数) * p^{floor(n/p)} * (floor(n/p))! 即n! (p-1)!^{floor(n/p)} * (最后一段不完整的周期) * p^{floor(n/p)} * (floor(n/p))!但是我们想要的是 F(n,p,pk)即提走所有 p 因子后的结果。注意在上面的式子中(floor(n/p))!这个部分里面可能还包含 p 的因子。所以我们需要递归地处理它 同时(p-1)! mod pk这一部分在一个完整的周期内长度为 p所有与 p 互质的数的乘积对 pk 取模是一个固定值。我们可以预处理这个值。由此得到递归式计算 F 的方法// 预处理计算 fact_no_p[i]表示从1到i中所有不是p的倍数的数的乘积模pk // 例如对于 p3, pk9 fact_no_p[5] (1*2*4*5) mod 9 40 mod 9 4 // 通常我们只需预处理到 pk 的长度。 long long F(long long n, int p, long long pk) { if (n 0) return 1; long long res 1; // 完整周期的贡献共有 n/pk 个完整周期每个周期乘积为 fact_no_p[pk] long long cycle_cnt n / pk; long long cycle_rem n % pk; res fast_pow(fact_no_p[pk], cycle_cnt, pk); // 完整周期的乘积 res res * fact_no_p[cycle_rem] % pk; // 最后一个不完整周期的贡献 // 递归处理 (n/p)! 中提走 p 因子后的部分 res res * F(n / p, p, pk) % pk; return res; }这里fact_no_p[pk]需要预处理。注意pk可能比较大比如 p^k但通常题目中 p^k 不会太大否则计算量剧增预处理是可行的。4.4 ExLucas 的完整组装现在我们可以组装出完整的 ExLucas 算法了。// 扩展欧几里得求逆元用于合并中国剩余定理 void exgcd(long long a, long long b, long long x, long long y) { if (b 0) { x 1; y 0; return; } exgcd(b, a % b, y, x); y - a / b * x; } long long inv(long long a, long long mod) { long long x, y; exgcd(a, mod, x, y); return (x % mod mod) % mod; } // 计算 C(n, m) mod p^k (p是质数) long long C_mod_pk(long long n, long long m, int p, long long pk) { if (m n) return 0; // 1. 计算p因子的指数 long long e get_g(n, p) - get_g(m, p) - get_g(n - m, p); if (e k) return 0; // 这里k是pk中p的指数需要作为参数传入或计算 // 2. 计算F函数值 long long fn F(n, p, pk); long long fm F(m, p, pk); long long fnm F(n - m, p, pk); // 3. 计算分数部分 long long frac fn * inv(fm, pk) % pk * inv(fnm, pk) % pk; // 4. 乘上p的幂 long long p_pow_e 1; for (int i 0; i e; i) p_pow_e p_pow_e * p % pk; return frac * p_pow_e % pk; } // 中国剩余定理合并 long long CRT(const vectorlong long a, const vectorlong long m) { long long M 1, res 0; for (long long mi : m) M * mi; for (int i 0; i a.size(); i) { long long Mi M / m[i]; long long inv_Mi inv(Mi, m[i]); res (res a[i] * Mi % M * inv_Mi) % M; } return res; } // ExLucas 主函数 long long exlucas(long long n, long long m, long long P) { if (m n) return 0; // 分解P vectorpairlong long, int factors; // 存储 (质因子p, 指数k) long long temp P; for (int i 2; (long long)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); vectorlong long a, mods; for (auto [p, k] : factors) { long long pk 1; for (int i 0; i k; i) pk * p; // 对于每个质因子幂计算 C(n, m) mod p^k long long ai C_mod_pk(n, m, p, pk); a.push_back(ai); mods.push_back(pk); } // 用中国剩余定理合并所有结果 return CRT(a, mods); }5. ExLucas实战复杂度、优化与避坑指南ExLucas 是一个强大的工具但它的实现比经典 Lucas 复杂得多也慢得多。理解其复杂度并掌握优化技巧至关重要。5.1 算法复杂度分析假设模数 P 分解为 t 个质因子幂最大的 p^k 记为max_pk。预处理对于每个质因子幂 p^k需要预处理fact_no_p数组长度为pk。总预处理复杂度为 O(Σ pk)最坏情况接近 O(P)但通常 pk 不会太大。计算 F 函数F(n, p, pk)是递归的递归深度为 O(log_p n)每层需要一次快速幂计算完整周期的贡献和一次乘法。单次计算 F 的复杂度约为 O(log_p n * log pk)。计算 C_mod_pk需要计算三次 F 函数和若干次简单运算。对于每个质因子幂计算复杂度主要取决于 F 函数。总体ExLucas 的复杂度与 P 的质因子分解、最大的 pk 大小以及 n 的大小有关。它是一个多项式时间算法但对于大的 P 和 n计算量依然可观。在竞赛中通常 P ≤ 1e6n, m ≤ 1e18 的场景下是可以接受的。5.2 常见优化策略预处理优化fact_no_p数组对于每个不同的(p, pk)都需要单独预处理。如果题目需要多次调用 ExLucas比如 T 组查询且模数 P 不变那么可以对 P 的每个质因子幂只预处理一次fact_no_p数组并存储起来避免重复计算。快速幂优化在F函数中计算fast_pow(fact_no_p[pk], cycle_cnt, pk)。当cycle_cnt很大时快速幂是必要的。可以检查如果cycle_cnt为 0 或 1直接返回结果避免快速幂调用。递归转迭代F函数本质上是递归的但可以写成迭代形式避免递归开销。不过递归深度为 log_p n通常不是瓶颈代码清晰更重要。提前返回在C_mod_pk中一旦计算出指数 e k可以立即返回 0。在get_g函数中如果发现 e 已经很大也可以提前终止计算。5.3 实战中的典型“坑”与解决方案坑1模数 P1组合数对 1 取模结果永远是 0。需要在 ExLucas 入口处特判。坑2n 或 m 为 0组合数 C(n,0)1, C(0,m)0 (m0)。在C_mod_pk和exlucas入口处都应处理这些边界情况。坑3p^k 较大时的预处理fact_no_p数组需要开到pk的大小。如果pk很大比如超过 1e7预处理可能会超内存。这时ExLucas 算法可能不再适用需要考虑其他方法或者题目数据可能保证了pk不会太大。在竞赛中如果模数 P 很大如 1e9但它是质数绝对不要用 ExLucas直接用普通逆元或 Lucas如果 n,mp。ExLucas 适用于 P 是合数且其质因子幂不大的情况。坑4中国剩余定理合并时的溢出在 CRT 函数中计算M * mi和Mi M / m[i]时如果模数 P 很大比如接近 1e18乘法可能会溢出long long。需要使用__int128或慢速乘龟速乘来处理大数乘法取模。// 使用 __int128 的 CRT 安全版本 long long CRT_safe(const vectorlong long a, const vectorlong long m) { __int128 M 1, res 0; for (long long mi : m) M * mi; for (int i 0; i a.size(); i) { __int128 Mi M / m[i]; __int128 inv_Mi inv(Mi, m[i]); // inv 函数需要支持 __int128 或 long long res (res (__int128)a[i] * Mi % M * inv_Mi) % M; } return (long long)res; }坑5逆元不存在在C_mod_pk中我们计算inv(fm, pk)和inv(fnm, pk)。由于F函数返回的值已经剔除了 p 因子所以它一定与pk互质吗是的因为F是剔除了所有 p 因子后剩余部分的乘积而pk的质因子只有 p所以F与pk互质逆元一定存在。这是算法正确性的关键保障。5.4 调试技巧ExLucas 代码较长容易写错。调试时可以从简单案例开始测试 P 为质数的情况与经典 Lucas 定理的结果对比。测试 P 为质数幂的情况如计算 C(5,2) mod 4, mod 8, mod 9 等可以手算验证。测试 P 为两个小质数乘积的情况如 P6计算一些小的组合数。使用暴力计算当 n,m 很小时来验证 ExLucas 的结果。写出 ExLucas 是一次很好的锻炼它融合了质因数分解、阶乘处理、递归、逆元、中国剩余定理等多个数论知识点。虽然代码有点长但一旦理解透彻就能应对模数为合数的组合数取模这一大类问题。在实际项目中如果模数是固定的可以预先分解并预处理使得每次查询的复杂度降低到 O(log n) 级别非常实用。