组合数取模进阶：从Lucas定理到ExLucas的算法精解

1. 从一道经典面试题说起组合数取模的“陷阱”如果你在准备算法面试或者刷过一些在线评测平台的数论题目大概率遇到过这样一道题计算组合数 C(n, m) 对一个大质数 P比如 1e97取模的结果。很多人的第一反应是这还不简单直接用公式 C(n, m) n! / (m! * (n-m)!)然后分别计算阶乘再用费马小定理求逆元不就行了这个思路在 P 是质数且足够大大于 n时确实是标准且高效的解法。但问题往往就藏在这个“但是”里。如果题目条件变成了计算 C(n, m) mod p其中 p 不一定是一个大质数它可能是一个合数甚至可能小到和 n、m 同数量级。这时你精心准备的阶乘逆元模板就会瞬间失效。为什么核心矛盾在于“逆元”的存在性。在模运算中一个数 a 关于模 p 的逆元 a⁻¹是满足 a * a⁻¹ ≡ 1 (mod p) 的整数。根据数论知识a 在模 p 下有逆元的充要条件是 a 与 p 互质即 gcd(a, p) 1。在我们计算 n! / (m! * (n-m)!) mod p 时我们实际上需要计算 n! * inv(m!) * inv((n-m)!) mod p。如果 p 是一个质数并且 p n, m, (n-m)那么 n!, m!, (n-m)! 都不可能含有质因子 p因此它们都与 p 互质逆元一定存在费马小定理a^(p-2) mod p或扩展欧几里得算法都可以求出逆元。然而一旦 p 不是质数或者 p 虽然质但小于等于 n情况就复杂了。例如计算 C(5, 2) mod 3。5! 120 2! 2 3! 6。120 / (26) 10 10 mod 3 1。但如果试图用逆元法在模 3 下2 的逆元是 2因为 224≡1 mod 36 的逆元呢6 和 3 不互质gcd(6,3)36 在模 3 下没有逆元因为 6 ≡ 0 (mod 3)任何数乘以 0 都是 0不可能等于 1。这时直接套用公式就会因为试图计算一个不存在的逆元而导致错误。这就是Lucas 定理和扩展 Lucas 定理ExLucas登场的背景。它们就是为了系统性地解决“模数 p 不是质数或者 p 较小”时组合数取模的计算问题。简单来说Lucas 定理是解决“p 是质数但可能很小”的利器而 ExLucas 则是解决“p 是任意正整数合数”的终极方案。理解它们不仅是解决特定题目的需要更是深入理解模运算、中国剩余定理CRT等数论核心思想的绝佳实践。2. Lucas 定理当模数是质数时的“降维打击”我们先来啃第一块硬骨头Lucas 定理。它的应用场景非常明确模数 p 是一个质数。至于 p 和 n, m 的大小关系定理本身并不要求 p n这正是它强大之处。2.1 定理陈述与核心思想Lucas 定理的表述非常简洁 对于质数 p 有 C(n, m) ≡ Π C(n_i, m_i) (mod p) 其中n_i 和 m_i 分别是 n 和 m 在 p 进制下的各个位数。 也就是说如果我们把 n 和 m 写成 p 进制的形式 n n_k * p^k n_{k-1} * p^{k-1} ... n_1 * p n_0 m m_k * p^k m_{k-1} * p^{k-1} ... m_1 * p m_0 那么组合数 C(n, m) 模 p 的值等于所有这些 p 进制位上对应的小组合数 C(n_i, m_i) 模 p 的乘积。这个定理的精妙之处在于“降维”。它将一个可能很大的 n 和 m 的组合数计算分解成了若干个很小的、在 p 进制下的“局部”组合数计算。因为 n_i 和 m_i 都是小于 p 的数p进制数的性质所以计算 C(n_i, m_i) mod p 就回到了我们熟悉的、可以用阶乘逆元安全计算的情形因为此时 n_i, m_i p 保证了阶乘与 p 互质。一个简单的例子计算 C(12, 5) mod 7。 首先p7 是质数。将 12 和 5 转化为 7 进制 12 (十进制) 1 * 7^1 5 * 7^0 (1, 5)_7 5 (十进制) 0 * 7^1 5 * 7^0 (0, 5)_7 根据 Lucas 定理 C(12, 5) mod 7 ≡ C(1, 0) * C(5, 5) mod 7 计算小组合数C(1,0)1 C(5,5)1。 所以结果为 1 * 1 1。 我们可以验证C(12,5)792 792 mod 7 792 - 113*7 792 - 791 1。结果正确。再看一个例子其中 m 的某一位大于 n 的对应位计算 C(10, 4) mod 5。 p5。转换进制 10 2 * 5^1 0 * 5^0 (2, 0)_5 4 0 * 5^1 4 * 5^0 (0, 4)_5 根据定理C(10,4) mod 5 ≡ C(2,0) * C(0,4) mod 5。 C(2,0)1 但 C(0,4) 是多少从组合数学定义从 0 个元素中选 4 个这是不可能的所以 C(0,4)0。 因此结果为 1 * 0 0。 验证C(10,4)210 210 mod 5 0。这也揭示了 Lucas 定理的一个关键点如果在 p 进制下m 的某一位数值 m_i 大于 n 的对应位 n_i那么 C(n_i, m_i) 0从而导致整个乘积为 0。这意味着 C(n, m) 能被 p 整除。2.2 算法实现与代码解析理解了原理实现起来就清晰了。算法步骤如下边界条件如果 m n 组合数为 0。循环分解当 n 0 或 m 0 时循环获取 n 和 m 在模 p 下的最后一位即 n % p 和 m % p。局部计算如果某次循环中m_i n_i直接返回 0。累乘否则用预处理的阶乘和逆元表计算 C(n_i, m_i) mod p并与结果相乘再取模。降维n / p, m / p 继续处理更高位。预处理为了快速计算 C(n_i, m_i)我们需要预处理出 0 到 p-1 的阶乘数组 fac[] 和对应的逆元数组 invFac[]。这里有一个关键的预处理优化因为 p 是质数且 n_i p我们可以用 O(p) 的时间预处理所有阶乘和阶乘逆元。计算逆元时通常采用费马小定理a^(p-2) mod p利用快速幂。更高效的方法是线性求逆元先求出 fac[p-1] 的逆元然后倒推回去。invFac[i] invFac[i1] * (i1) % p下面是 Lucas 算法的 C 实现核心代码#include iostream using namespace std; typedef long long ll; // 快速幂用于求逆元 (当p是质数时 inv(a) a^(p-2) mod p) ll qpow(ll a, ll b, ll p) { ll res 1; while (b) { if (b 1) res res * a % p; a a * a % p; b 1; } return res; } // 预处理阶乘和阶乘逆元范围 [0, p-1] ll fac[100005], invFac[100005]; void init(ll p) { fac[0] 1; for (int i 1; i p; i) { fac[i] fac[i-1] * i % p; } // 费马小定理求最大阶乘的逆元 invFac[p-1] qpow(fac[p-1], p-2, p); // 线性倒推阶乘逆元 for (int i p-2; i 0; --i) { invFac[i] invFac[i1] * (i1) % p; } } // 计算小组合数 C(a, b) mod p, 其中 a, b p ll smallC(ll a, ll b, ll p) { if (b a) return 0; // C(a, b) a! / (b! * (a-b)!) return fac[a] * invFac[b] % p * invFac[a-b] % p; } // Lucas 定理主函数 ll lucas(ll n, ll m, ll p) { if (m 0) return 1; // C(n, 0) 1 // 递归或循环形式将问题分解为p进制下的每一位 return smallC(n % p, m % p, p) * lucas(n / p, m / p, p) % p; } int main() { ll n, m, p; // 假设 p 是质数 cin n m p; init(p); // 预处理模 p 下的阶乘表 cout lucas(n, m, p) endl; return 0; }注意上述代码中lucas函数采用了递归形式非常直观地反映了定理的分解思想。在递归调用lucas(n/p, m/p, p)时它继续处理更高 p 进制位。递归的终止条件是m 0返回 1因为 C(anything, 0) 1。你也可以用循环来实现本质是一样的。2.3 时间复杂度分析与适用场景Lucas 算法的时间复杂度主要取决于两部分预处理O(p)用于构建阶乘和逆元表。这是算法的主要开销。递归/循环计算O(log_p n)因为我们需要把 n 和 m 逐位除以 p直到它们变为 0。因此总时间复杂度为 O(p log_p n)。这决定了 Lucas 定理的适用边界模数 p 不能太大。通常p 在 10^5 到 10^6 量级时预处理 O(p) 是可行的。如果 p 很大比如接近 10^9预处理 O(p) 的时空开销都无法承受此时 Lucas 定理就不再适用。幸运的是当 p 很大且为质数时通常也满足 p n 我们可以直接使用普通的阶乘逆元法复杂度 O(n)而无需 Lucas 定理。所以Lucas 定理的典型应用场景是p 是质数但 p 相对较小例如 1e5 以内而 n 和 m 可以非常大甚至超过 1e18。它巧妙地将大数运算转化为了小数运算。3. 扩展 Lucas 定理ExLucas应对合数模数的通用解法现在我们来面对更一般、也更复杂的情况模数 p 不是一个质数而是一个任意的正整数。这就是ExLucasExtended Lucas要解决的问题。虽然名字里带“Lucas”但它的核心思想与 Lucas 定理截然不同更像是中国剩余定理CRT和质因数分解的一个经典应用。3.1 问题拆解中国剩余定理的视角ExLucas 的核心思路基于这样一个事实如果我们可以解决模数为质数幂p^k的组合数取模问题那么对于任意合数模数我们都可以通过质因数分解和中国剩余定理CRT将问题组合起来。具体来说假设我们要计算 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}。这是整个 ExLucas 算法中最困难的部分。结果合并利用中国剩余定理CRT将 t 个同余方程的解合并得到最终关于模 P 的解。所以ExLucas 的关键和难点就落在了第 2 步如何计算 C(n, m) mod p^kp是质数k是正整数。3.2 攻克核心子问题C(n, m) mod p^k计算 C(n, m) mod p^k 不能直接使用阶乘逆元因为分母中的阶乘可能与 p 不互质导致逆元不存在。ExLucas 采用了一种“剥离质因子 p”的策略。我们回顾组合数公式C(n, m) n! / (m! * (n-m)!) 计算这个式子模 p^k 的值等价于计算 [n!]{p^k} / ([m!]{p^k} * [(n-m)!]_{p^k}) mod p^k 但除法在模 p^k 下不一定可行。我们需要将阶乘中所有质因子 p 提取出来使得剩下的部分与 p 互质从而可以求逆元。定义函数F(n, p, pk) 它计算 n! 中除去所有质因子 p 后剩余部分模 pk 的值。同时我们还需要知道 n! 中质因子 p 的个数记为G(n, p)。那么n! 可以表示为n! p^{G(n, p)} * F(n, p, pk) * (一些与p互质的因子)。 同理m! 和 (n-m)! 也可以这样表示。因此C(n, m) n! / (m! * (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))]现在指数部分 G(n,p) - G(m,p) - G(n-m,p) 是一个非负整数根据组合数的整数性。如果这个指数大于等于 k那么整个 C(n, m) 模 p^k 就是 0因为被 p^k 整除了。否则我们只需要计算剩余部分[F(n, p, pk) / (F(m, p, pk) * F(n-m, p, pk))] mod pk。关键来了现在分子分母中的 F 函数值都是与 p 互质的因为我们把所有的因子 p 都提走了所以在模 pk 下分母的逆元存在我们可以用扩展欧几里得算法求出逆元。所以问题转化为两个子函数G(n, p)计算 n! 中质因子 p 的个数。有一个经典的公式G(n, p) floor(n/p) floor(n/p^2) floor(n/p^3) ...。例如n10, p2 G floor(10/2)floor(10/4)floor(10/8)5218。这很好理解1到10中2的倍数有5个贡献至少1个24的倍数有2个额外再贡献1个28的倍数有1个再额外贡献1个2。F(n, p, pk)计算 n! 中移除所有 p 因子后模 pk 的值。这是 ExLucas 实现中最繁琐的部分。3.3 F(n, p, pk) 的计算递归与周期性的利用计算 F(n, p, pk) 不能暴力从 1 乘到 n因为 n 可能很大。我们需要利用模 pk 运算的周期性。考虑 n! 1 * 2 * 3 * ... * n。 我们把所有 p 的倍数都挑出来p, 2p, 3p, ..., floor(n/p)*p。 那么 n! 可以写成两部分乘积 n! (1 * 2 * ... * (p-1) * (p1) * ... ) * (p * 2p * 3p * ... * floor(n/p)*p) 第一部分是不含因子 p 的数的乘积第二部分是所有含因子 p 的数的乘积。对于第二部分我们可以提取出公因子 p p * 2p * 3p * ... * floor(n/p)*p p^{floor(n/p)} * (1 * 2 * 3 * ... * floor(n/p)) p^{floor(n/p)} * floor(n/p)!看到了吗第二部分里又出现了一个阶乘 floor(n/p)!。这启发我们使用递归。更精确地我们考虑模 pk 下的一个“周期”。因为模 pk 一个数与 p 互质其逆元存在。我们可以把 1 到 n 的数分组第一组1, 2, ..., pk。这是一个完整的周期。第二组pk1, pk2, ..., 2pk。以此类推最后可能有一个不完整的周期。在每个完整的周期内所有与 p 互质的数的乘积模 pk 是一个常数记为prod。我们可以预处理这个prod。例如对于 p3, k2, pk9 周期 1~9 中与 3 互质的数是 {1,2,4,5,7,8}它们的乘积 12457*82240 2240 mod 9 8。所以 prod 8。那么F(n, p, pk) 可以递归计算 F(n, p, pk) [ (prod)^{floor(n/pk)} * F(n mod pk, p, pk) ] * F(floor(n/p), p, pk) mod pk解释(prod)^{floor(n/pk)}计算所有完整周期内与 p 互质部分乘积的贡献。F(n mod pk, p, pk)计算最后一个不完整周期内与 p 互质部分乘积的贡献。F(floor(n/p), p, pk)这是递归部分处理从所有数中提取出因子 p 后剩下的那个阶乘 floor(n/p)!。注意这里递归调用 F 参数变成了 floor(n/p) 因为我们要处理的是p * 2p * ... * floor(n/p)*p中提取 p 后剩下的floor(n/p)!。递归的边界条件是 F(0, p, pk) 1。3.4 ExLucas 算法完整实现步骤与代码框架结合以上所有分析我们可以勾勒出 ExLucas 的完整算法流程主函数exlucas(n, m, P):如果 m n 返回 0。对模数 P 进行质因数分解得到向量vectorpairll, ll factors 每个元素是 (p, k) 表示质因子 p 和它的幂次 k。初始化两个数组a[]和m[] 用于存储中国剩余定理的方程。对于每个 (p, k) a. 计算pk pow(p, k)。 b. 调用子函数C_mod_pk(n, m, p, k, pk) 得到结果res 即 C(n, m) mod pk。 c. 将res存入a[i] 将pk存入m[i]。调用中国剩余定理函数crt(a, m) 合并所有同余方程得到最终答案C(n, m) mod P。子函数C_mod_pk(n, m, p, k, pk):计算指数e G(n,p) - G(m,p) - G(n-m,p)。如果e k 直接返回 0。计算a F(n, p, pk) mod pk。计算b F(m, p, pk) mod pk。计算c F(n-m, p, pk) mod pk。计算res a * inv(b, pk) % pk * inv(c, pk) % pk。其中inv(x, pk)使用扩展欧几里得算法计算因为 x 与 p 互质所以逆元存在。计算res res * pow_mod(p, e, pk) % pk。这里pow_mod是快速幂。返回res。子函数G(n, p): 使用循环计算while (n) { cnt n/p; n / p; }。子函数F(n, p, pk):如果 n 0 返回 1。预处理周期乘积prod计算 1 到 pk 之间所有与 p 互质的数的乘积模 pk。设res pow_mod(prod, n/pk, pk)。计算剩余部分rem n % pk 然后计算 1 到 rem 之间所有与 p 互质的数的乘积模 pk记为left。递归计算res res * left % pk * F(n/p, p, pk) % pk。返回res。由于代码较长这里给出一个高度概括的 C 框架并强调几个极易出错的实现细节#include iostream #include vector using namespace std; typedef long long ll; // 快速幂 ll qpow(ll a, ll b, ll p) { ... } // 扩展欧几里得求逆元 (ax by gcd(a,b)) ll exgcd(ll a, ll b, ll x, ll y) { ... } ll inv(ll a, ll p) { // 求 a 在模 p 下的逆元gcd(a,p)1 ll x, y; exgcd(a, p, x, y); return (x % p p) % p; } // 计算 n! 中质因子 p 的个数 ll G(ll n, ll p) { ll cnt 0; while (n) { cnt n / p; n / p; } return cnt; } // 计算 F(n, p, pk) ll F(ll n, ll p, ll pk) { if (n 0) return 1; // 预处理周期乘积 prod (1..pk中与p互质的数的积 % pk) // 这里需要预处理可以全局缓存或每次计算。假设我们有一个函数 get_prod(p, pk) ll prod get_prod(p, pk); // 需要实现 ll res qpow(prod, n / pk, pk); ll rem n % pk; // 计算剩余部分 left (1..rem中与p互质的数的积 % pk) ll left 1; for (ll i 1; i rem; i) { if (i % p ! 0) { left left * i % pk; } } res res * left % pk; // 递归处理 res res * F(n / p, p, pk) % pk; return res; } // 计算 C(n, m) mod p^k ll C_mod_pk(ll n, ll m, ll p, ll k, ll pk) { if (m n) return 0; ll e G(n, p) - G(m, p) - G(n - m, p); if (e k) return 0; // 被 p^k 整除 ll a F(n, p, pk); ll b F(m, p, pk); ll c F(n - m, p, pk); ll res a * inv(b, pk) % pk * inv(c, pk) % pk; res res * qpow(p, e, pk) % pk; return res; } // 质因数分解 vectorpairll, ll factorize(ll P) { ... } // 中国剩余定理 (模数两两互质) ll crt(const vectorll a, const vectorll m) { ... } // ExLucas 主函数 ll exlucas(ll n, ll m, ll P) { if (m n) return 0; vectorpairll, ll factors factorize(P); vectorll rems, mods; for (auto [p, k] : factors) { ll pk 1; for (int i 0; i k; i) pk * p; ll r C_mod_pk(n, m, p, k, pk); rems.push_back(r); mods.push_back(pk); } // 使用中国剩余定理合并 return crt(rems, mods); }关键实现细节与踩坑点get_prod(p, pk)的预处理与缓存对于固定的 (p, pk) 周期乘积prod是常数。在多次查询 ExLucas 时应该缓存这个值避免重复计算。计算prod时需要遍历 1 到 pk 跳过 p 的倍数相乘并取模。当 pk 较大时比如超过 1e6这个预处理可能成为瓶颈这也是 ExLucas 的主要性能限制之一。递归深度F(n, p, pk)的递归深度是 O(log_p n) 对于极大的 n如 1e18是可以接受的。中国剩余定理的实现需要保证所有pk两两互质这是由质因数分解保证的。CRT 的合并过程可能涉及大数乘法需要注意使用__int128或类似技术防止中间结果溢出。性能瓶颈ExLucas 的复杂度主要受限于模数 P 的最大质因子幂pk。如果 P 有一个质因子 p 的幂次 k 很高导致pk很大那么预处理prod和计算F函数中的循环部分计算left开销会很大。在实践中如果pk超过 1e7 算法可能会非常慢。4. 实战对比与场景选择Lucas vs. ExLucas vs. 普通逆元法理解了原理和实现我们最后来梳理一下面对不同的模数 P 应该如何选择正确的算法。这往往是比赛或面试中决定效率甚至能否AC的关键。场景特征推荐算法理由与注意事项P 是质数且 P n, m普通阶乘逆元法这是最简单、最快速的情况。直接预处理 1~n 的阶乘和逆元O(n) 预处理O(1) 查询。复杂度最低。P 是质数但 P 可能 ≤ n, mLucas 定理Lucas 定理的经典场景。它能处理 P 小于 n 的情况复杂度 O(P log_P n)。注意 P 不能太大通常 P ≤ 1e6否则预处理 O(P) 会超时或超内存。P 是合数扩展 Lucas 定理 (ExLucas)唯一通用的解法。通过质因数分解和CRT将问题转化为多个模质数幂的子问题。性能取决于 P 的最大质因子幂p^k。如果p^k很大1e7算法会变慢。P 不是质数但 n, m 非常小 (≤ 50)暴力计算或动态规划对于非常小的 n, m 可以直接用定义计算或者用杨辉三角递推 C[n][m] C[n-1][m-1] C[n-1][m] 结果对 P 取模。简单可靠。需要多次查询不同的 n, m 模数 P 固定根据 P 的类型选择上述方法并做好预处理如果使用普通逆元法或 Lucas预处理阶乘和逆元表一次即可。如果使用 ExLucas 可以对每个质因子幂p^k预处理其周期乘积prod 能大幅提升多次查询的效率。个人踩坑经验分享永远先判断模数 P 的类型。拿到题目第一步不是默写模板而是分析 P 的范围和性质。很多题目会刻意设置 P1e97大质数来引导你使用普通逆元法或者设置 P10007小质数来考察 Lucas 又或者设置 P999911658合数分解后是 2 * 3 * 4679 * 35617来考察 ExLucas。ExLucas 的调试技巧。ExLucas 代码复杂容易写错。调试时可以先用小数据暴力计算组合数验证。重点检查三个函数G质因子个数、F提取p后的阶乘、C_mod_pk。可以单独为F函数写测试比如验证 F(10, 2, 8) 是否等于 (10! 中除去所有2后模8的值)。注意数据范围与溢出。在计算阶乘、乘积、快速幂时即使对pk取模中间结果也可能溢出long long。在关键乘法处使用__int128临时转换或手动实现慢速乘mul_mod是稳妥的做法尤其是在 CRT 合并步骤。时间复杂度心里有数。如果题目中 n, m 高达 1e18 P1009质数那么 Lucas 定理的 O(P) 预处理是可行的。但如果 P1e97 n,m1e5 那么就应该用普通逆元法O(n)。如果 P 是一个有较大质因子幂的合数如 10^9 以内ExLucas 可能会超时需要审视是否有其他数学性质可以利用。最后无论是 Lucas 还是 ExLucas它们都是工具。理解其背后的数论原理——模运算、逆元、质因数分解、中国剩余定理——远比记住模板代码更重要。当你透彻理解了为什么在模合数下直接求逆元会失败以及 ExLucas 是如何通过“剥离质因子”来解决这个问题的你就能在面对千变万化的组合数取模问题时灵活地选用或组合这些工具甚至推导出属于自己的优化方法。