扩展卢卡斯定理：非质数模数下组合数取模的算法实现与原理

1. 项目概述：当组合数遇上非质数模数

在算法竞赛和数论研究中，计算组合数 C(n, m) 对一个大整数 P 取模的结果，是一个经典且高频的需求。当模数 P 是一个质数时，我们有成熟的卢卡斯定理（Lucas Theorem）可以高效解决，其核心思想是将 n 和 m 按 P 进制分解，递归求解。然而，现实世界（尤其是题目里）往往没那么“理想”——模数 P 经常不是一个质数，它可能是一个合数，甚至是质数幂的乘积。这时，标准的卢卡斯定理就失效了。这正是“扩展卢卡斯定理”（exLucas）要解决的痛点：它允许我们在模数为任意正整数（不一定是质数）的情况下，计算组合数取模。

这个“模板”项目的核心，就是理解和实现 exLucas 算法。它不是一个单一的公式，而是一个精巧的算法框架，巧妙地将中国剩余定理、质因数分解、阶乘模质数幂的计算等多个数论工具串联起来。对于需要处理模数非质数的组合数问题，比如某些计数类动态规划、多项式系数计算，或者直接就是一道要求计算 C(n, m) mod P 的题目，掌握 exLucas 是必不可少的技能。

2. 算法核心思路拆解：化整为零，分而治之

exLucas 算法的核心思想是“分解与合并”。面对一个非质数模数 P，我们无法直接在其上定义模逆元（因为不是所有数都有逆元），所以需要转换思路。

2.1 第一步：模数质因数分解

假设我们要计算 C(n, m) mod P，其中 P 是任意正整数。 首先，对模数 P 进行质因数分解：P = p1^k1 * p2^k2 * ... * pt^kt，其中 pi 是互不相同的质数，ki 是正整数。

算法的目标转化为：分别求出C(n, m) mod pi^ki对于每一个 i 的结果，然后再利用中国剩余定理将这些结果合并，得到最终C(n, m) mod P。

为什么这么做？因为模一个质数幂pi^ki的环境，虽然比模质数复杂，但比模任意合数要规则得多。我们可以在每个pi^ki的“局部”环境下解决问题。

2.2 第二步：解决子问题 C(n, m) mod p^k

这是整个算法最难也是最核心的部分。对于某个质数 p 和指数 k，我们需要计算C(n, m) mod p^k。

组合数公式为：C(n, m) = n! / (m! * (n-m)!)。 在模p^k下，分母m!和(n-m)!可能含有因子 p，因此它们可能没有模p^k下的逆元（因为与模数不互素）。我们不能直接计算除法。

exLucas 的巧妙之处在于对阶乘进行“净化”处理：将阶乘 n! 写成两部分乘积的形式：n! = p^a * b。 其中，a是 n! 中质因子 p 的个数，b是 n! 中剔除了所有因子 p 之后剩下的部分，并且b与 p 互质。

例如，对于 p=2, k=3，计算 10! mod 8。 10! = 12345678910。 其中因子 2 出现在 2, 4=2^2, 6=23, 8=2^3, 10=25。 总共有 1+2+1+3+1 = 8 个因子 2，所以 a=8。 剔除所有因子 2 后，剩下的数相乘：1131537195。计算这个乘积 mod 8，得到 b。

这样，组合数就可以表示为：C(n, m) = [n! / (m! * (n-m)!)] = [p^{a_n} * b_n] / [p^{a_m} * b_m * p^{a_{n-m}} * b_{n-m}] = p^{a_n - a_m - a_{n-m}} * (b_n * inv(b_m) * inv(b_{n-m}))

其中，inv(x)表示 x 在模p^k下的逆元。由于 b_n, b_m, b_{n-m} 都与 p 互质，所以它们在模p^k下存在逆元，可以安全计算。

因此，问题进一步转化为两个子问题：

1. 计算指数e = a_n - a_m - a_{n-m}。
2. 计算“净化”后的部分b_n * inv(b_m) * inv(b_{n-m}) mod p^k。

如果e >= k，那么p^e mod p^k = 0，整个组合数模p^k就是 0。 否则，最终结果就是(p^e * (b部分的结果)) mod p^k。

2.3 第三步：递归计算净化阶乘 b

如何高效计算F(n) = n! 中剔除所有因子 p 后的乘积 mod p^k（即上面的 b）？ 这里采用递归思想。

观察 n!，我们可以把 1 到 n 的数按是否包含因子 p 来分组：

• 不包含因子 p 的数：1, 2, ..., p-1, p+1, ..., 2p-1, ...
• 包含至少一个因子 p 的数：p, 2p, 3p, ..., floor(n/p) * p。

对于包含因子 p 的数，我们可以提取出一个公因子 p，剩下的部分就变成了1, 2, 3, ..., floor(n/p)。而这正好是F(floor(n/p))要处理的问题！

因此，递归公式如下：F(n) = [乘积_{i=1, i%p!=0}^{p^k} i]^{floor(n / p^k)} * [乘积_{i=1}^{n mod p^k} i (if i%p!=0)] * F(floor(n/p)) mod p^k

公式解读：

1. 乘积_{i=1, i%p!=0}^{p^k} i：这是一个周期乘积。在模p^k的意义下，从 1 到p^k这p^k个数中，所有不与 p 互质的数（即 p 的倍数）我们暂时不考虑（它们被归入递归部分）。剩下的p^{k-1}*(p-1)个数构成一个“周期块”，它们的乘积记为pre。因为模运算的周期性，n!中完整的周期块有floor(n / p^k)个，所以贡献是pre^{floor(n / p^k)}。
2. 乘积_{i=1}^{n mod p^k} i (if i%p!=0)：这是最后一个不完整的周期块中，不与 p 互质的数的乘积。
3. F(floor(n/p))：这是递归部分，处理所有被提取出来的因子 p 之后剩下的那个整数序列的净化阶乘。

同时，计算指数a（即 n! 中因子 p 的个数）有一个著名的公式（勒让德定理）：a = floor(n/p) + floor(n/p^2) + floor(n/p^3) + ...这个计算可以在递归计算F(n)的过程中顺便完成。

2.4 第四步：中国剩余定理合并结果

经过第二步，我们得到了 t 组同余方程：x ≡ C(n, m) (mod p1^k1)x ≡ C(n, m) (mod p2^k2)...x ≡ C(n, m) (mod pt^kt)

由于p1^k1, p2^k2, ..., pt^kt两两互质（因为 pi 是不同质数），满足中国剩余定理的应用条件。我们可以使用 CRT 求出唯一的x mod P，这个 x 就是C(n, m) mod P。

CRT 的合并公式为： 设M = P，Mi = M / pi^ki，ti是Mi在模pi^ki下的逆元（因为 Mi 与 pi^ki 互质，逆元存在）。 则解为：x ≡ Σ(ai * Mi * ti) (mod M)，其中ai = C(n, m) mod pi^ki。

3. 算法实现细节与关键步骤

理解了思路，我们来看具体的实现。实现 exLucas 主要需要三个函数：快速幂、扩展欧几里得求逆元、以及核心的calc函数。

3.1 辅助函数：快速幂与扩展欧几里得

这些是基础数论工具。

// 快速幂，计算 base^exp % mod long long qpow(long long base, long long exp, long long mod) { long long res = 1; while (exp) { if (exp & 1) res = res * base % mod; base = base * base % mod; exp >>= 1; } return res; } // 扩展欧几里得，求解 ax + by = gcd(a, b)，返回 gcd，并通过引用返回 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); return (x % mod + mod) % mod; // 保证返回正数 }

3.2 核心函数：计算 F(n) 和指数 a

这个函数对应思路中的第二步，计算对于特定的(p, pk)（pk = p^k），C(n, m) mod pk的值。 我们实现一个函数lucas_pk(long long n, long long m, long long p, long long pk)。

// 计算净化阶乘 F(n) mod pk，同时返回指数 a（通过引用） long long factorial_p(long long n, long long p, long long pk, long long &a) { if (n == 0) { a = 0; return 1; } // 计算一个完整周期块的乘积 pre long long res = 1; for (long long i = 1; i <= pk; ++i) { if (i % p) { // 只取与 p 互质的数 res = res * i % pk; } } res = qpow(res, n / pk, pk); // 完整周期块的贡献 // 计算最后一个不完整周期块 for (long long i = 1; i <= n % pk; ++i) { if (i % p) { res = res * i % pk; } } // 递归计算 F(floor(n/p))，并累加指数 a long long next_res = factorial_p(n / p, p, pk, a); res = res * next_res % pk; a += n / p; // 勒让德公式，累加因子 p 的个数 return res; } // 计算 C(n, m) mod pk long long C_pk(long long n, long long m, long long p, long long pk) { if (m > n) return 0; long long a_n = 0, a_m = 0, a_nm = 0; // 计算三个净化阶乘及对应的指数 long long f_n = factorial_p(n, p, pk, a_n); long long f_m = factorial_p(m, p, pk, a_m); long long f_nm = factorial_p(n - m, p, pk, a_nm); long long e = a_n - a_m - a_nm; if (e >= 0) { // 计算净化部分：f_n * inv(f_m) * inv(f_nm) mod pk long long res = f_n * inv(f_m, pk) % pk * inv(f_nm, pk) % pk; // 乘上 p^e res = res * qpow(p, e, pk) % pk; return res; } else { // 如果指数为负，说明分母中 p 的因子比分子多，整个数不是整数，但在模 pk 下我们按 0 处理？ // 实际上，在组合数定义中 m <= n，且为整数，e 不可能为负。这里为了安全可以返回 0 或报错。 return 0; } }

3.3 主函数：exLucas 与 CRT 合并

现在实现最终的exLucas函数，它负责质因数分解 P，对每个质因子幂调用C_pk，最后用 CRT 合并。

// 扩展卢卡斯定理主函数：计算 C(n, m) mod P long long exLucas(long long n, long long m, long long P) { if (m > n) return 0; long long tmp = P; vector<pair<long long, long long>> factors; // 存储 (质数 p, 幂次 pk) // 质因数分解 P for (long long i = 2; i * i <= tmp; ++i) { if (tmp % i == 0) { long long pk = 1; while (tmp % i == 0) { tmp /= i; pk *= i; } factors.emplace_back(i, pk); } } if (tmp > 1) { factors.emplace_back(tmp, tmp); } // 分别计算 C(n, m) mod pk vector<long long> a(factors.size()), mod(factors.size()); for (size_t i = 0; i < factors.size(); ++i) { long long p = factors[i].first; long long pk = factors[i].second; a[i] = C_pk(n, m, p, pk); mod[i] = pk; } // 中国剩余定理合并 long long res = 0, M = P; for (size_t i = 0; i < factors.size(); ++i) { long long Mi = M / mod[i]; long long ti = inv(Mi, mod[i]); // 求 Mi 在模 mod[i] 下的逆元 res = (res + a[i] * Mi % M * ti % M) % M; } return res; }

4. 边界处理、优化与常见问题

4.1 边界情况与细节处理

1. m > n或m < 0：组合数定义为 0，应在函数入口处判断并返回 0。

2. P = 1：任何数模 1 都为 0，可以直接返回 0。

3. p^k可能很大：在计算周期块乘积pre时，pk可能达到1e9甚至更大（虽然通常题目会控制），循环for (i=1; i<=pk; ++i)会超时。这是实现中的一个关键优化点。优化方法：我们不需要真的计算1..pk的完整乘积。注意到我们只关心i % p != 0的i，并且是在模pk下计算。我们可以预处理一个数组pre[pk]，但pk太大时不行。实际上，由于我们最终要计算pre^{floor(n/pk)}，而floor(n/pk)通常很小（因为 n 有限），我们可以不显式计算出整个pre，而是在递归时直接计算(n % pk)!的部分，并通过递归关系处理完整周期。上面给出的factorial_p实现是一种简化描述，在pk较大时，计算res = qpow(res, n / pk, pk)中的res（即pre）会成为瓶颈。更高效的实现是避免直接计算pre，而是通过递归将n缩小。更优的factorial_p实现思路：

   long long factorial_p(long long n, long long p, long long pk) { if (n == 0) return 1; // 递归计算 F(n) = F(n/p) * (周期块乘积)^{n/pk} * (剩余块乘积) long long res = factorial_p(n / p, p, pk); // 递归处理提取p后的部分 // 计算周期块和剩余块 long long period_prod = 1; // 这里需要高效计算 1..pk 中非p倍数的乘积 mod pk long long rem_prod = 1; // 计算 1..(n%pk) 中非p倍数的乘积 mod pk // ... 计算 period_prod 和 rem_prod ... res = res * qpow(period_prod, n / pk, pk) % pk; res = res * rem_prod % pk; return res; }

   计算period_prod和rem_prod如果pk不大（比如 <= 1e6），可以预处理。如果pk很大，可能需要更数学的方法，或者题目通常会保证pk（即 p^k）不会太大，使得O(pk)的循环可接受。这是算法的时间复杂度瓶颈之一，通常认为p^k在1e6量级内是可接受的。

4. 求逆元inv(f_m, pk)：务必确保f_m与pk互质（即不含因子 p），这是我们设计算法时保证的。使用扩展欧几里得求逆元是安全的。

4.2 时间复杂度分析

假设模数 P 分解为Π pi^ki。

• 质因数分解 P：O(√P)，但 P 通常不大（如 1e9 以内）。
• 对于每个质因子幂(p, pk)：
  • factorial_p递归深度为O(log_p n)。
  • 每次递归需要计算周期块乘积（若pk不大，可O(pk)预处理；若直接循环，则每次O(pk)）和剩余块乘积O(n%pk)。
  • 粗略估计，处理一个(p, pk)的时间复杂度约为O(pk + log_p n)。
• CRT 合并：O(t)，t 是质因子个数，很小。

因此，总时间复杂度主要取决于最大的pk值。如果最大的pk在百万级别，算法可以在合理时间内运行。

4.3 常见问题与调试技巧

1. 结果错误为 0：

   • 检查是否m > n导致直接返回 0。
   • 检查C_pk函数中，当指数e >= k时，返回p^e * ... mod pk，如果e >= k，p^e mod pk确实为 0，这是正确的。但需确认e计算是否正确（factorial_p中的a累加是否正确）。
   • 检查模数pk是否在运算过程中因为溢出变成了负数或奇怪的值。确保使用long long并在乘法后及时取模。

2. 运行超时：

   • 最大的可能是pk太大，导致计算周期块乘积pre的循环O(pk)太慢。确认题目约束，pk是否真的可能很大（如 > 1e7）。如果很大，需要实现上文提到的优化版factorial_p，避免直接循环计算pre。
   • 递归深度过大？log_p n通常很小，比如 n=1e18, p=2，深度也就 60左右，不是问题。

3. 逆元计算失败：

   • 在C_pk中计算inv(f_m, pk)前，理论上f_m与pk互质。如果求逆失败（扩展欧几里得返回的 gcd 不为1），说明factorial_p函数有 bug，没有正确剔除因子 p。

4. CRT 合并结果错误：

   • 检查inv(Mi, mod[i])是否计算正确，即Mi与mod[i]（即pk）是否互质。由于Mi = P / pk，而pk是p^k，Mi包含其他质因子，肯定与pk互质，所以逆元存在。
   • 检查合并公式res = (res + a * Mi % M * ti % M) % M;中的取模是否正确，确保每次加法乘法后都取模M，防止溢出。

一个实用的调试方法：用小的、手算可验证的样例进行测试。 例如，计算C(5, 2) mod 6。

• 分解6 = 2 * 3。
• 计算C(5,2)=10。
• 10 mod 2 = 0。
• 10 mod 3 = 1。
• 解同余方程组x ≡ 0 (mod 2), x ≡ 1 (mod 3)。
• 解得x ≡ 4 (mod 6)。所以C(5,2) mod 6 = 4。 用你的 exLucas 程序计算，看结果是否为 4。

5. 实战应用与扩展思考

exLucas 算法虽然原理和实现略显复杂，但它解决了模数非质数时组合数计算的根本问题。在算法竞赛中，它通常以“模板题”的形式出现，要求你实现它来计算C(n, m) mod P。理解其每一步的数学原理，对于应对可能的变化至关重要。

可能的变化包括：

1. 模数 P 很大，但质因子幂p^k较小：这是 exLucas 发挥作用的典型场景。算法效率取决于最大的pk。
2. 需要计算多次组合数：可以对每个质因子幂pk预处理周期块乘积pre，这样在多次调用C_pk时，可以避免重复计算pre，提升效率。
3. 与其它数论定理结合：有时题目可能要求计算C(n, m) mod P，但 P 不是固定的，或者需要处理更复杂的求和式。exLucas 可以作为其中一个模块。

个人实现心得：

• 理解优于记忆：尝试自己推导一遍factorial_p的递归公式，比死记硬背代码更有用。理解了“提取因子p”和“周期块”的概念，就能应对各种变体。
• 重视边界：仔细处理n=0, m=0, p^k=1等情况，并确保在C_pk中正确处理指数e为 0 或大于等于k的情况。
• 测试驱动：实现后，务必用多个小数据（包括边界数据）和暴力计算（对于小的 n, m, P）进行对比测试，确保正确性。再找一些标准题目（如洛谷 P4720）的样例进行测试。
• 复杂度心里有数：明确算法的瓶颈在于pk的大小。如果题目中pk可能很大（比如由大质数组成），需要意识到可能存在的性能问题，或者考虑题目是否保证了pk不会太大。

最后，exLucas 是数论工具链中重要的一环。它将组合数计算、阶乘模运算、中国剩余定理、质因数分解等知识串联起来。掌握它，不仅能解决一类特定问题，更能加深你对模运算、同余方程和递归思想的理解。在遇到模数非质数的组合数问题时，你便可以有条不紊地将其分解为质因子幂上的子问题，然后逐个击破，最后合并得到答案。这种“化整为零”的思路，在算法设计中也是一种强大的策略。