1. 从一道“毒瘤”题说起为什么我们需要 exLucas如果你在算法竞赛或者数论学习的路上走得足够远一定会遇到这样一类问题计算一个组合数 ( C_n^m \mod P )其中 ( n ) 和 ( m ) 可以非常大比如 ( 10^{18} )但模数 ( P ) 不一定是质数。更具体地说( P ) 可能是一个合数甚至是一个质数的幂次方如 ( p^k )。当你兴冲冲地准备用 Lucas 定理或者预处理阶乘逆元时会发现这些经典方法全都失效了。Lucas 定理要求模数是质数而预处理阶乘要求模数与待计算的数互质以便使用逆元。当 ( P ) 是一个合数特别是包含小质因子时阶乘中会出现大量与 ( P ) 不互质的项导致逆元不存在直接计算路径被堵死。这就是 exLucas扩展卢卡斯定理登场的场景。它不要求模数 ( P ) 是质数只要求我们能对 ( P ) 进行质因数分解。其核心思想非常巧妙既然整体不好算我们就“分而治之”。将模数 ( P ) 分解为若干个质数幂的乘积 ( P p_1^{k_1} \cdot p_2^{k_2} \cdots p_t^{k_t} )然后分别计算组合数 ( C_n^m ) 对每一个 ( p_i^{k_i} ) 取模的结果最后利用中国剩余定理将这些结果“组装”回模 ( P ) 下的答案。这个算法完美解决了模数为合数时组合数求值的问题是数论工具箱中一件处理边界情况的利器。虽然原理相对复杂但一旦掌握面对相关题目便能游刃有余。2. 核心思路拆解化整为零与各个击破exLucas 的整个流程可以清晰地分为三个层次理解这个层次是掌握其证明和实现的关键。2.1 第一层中国剩余定理降维这是整个算法的战略框架。给定模数 ( P )我们首先对其进行质因数分解( P \prod_{i1}^{t} p_i^{k_i} )。根据中国剩余定理如果我们能计算出方程组 [ \begin{cases} C_n^m \equiv a_1 \pmod{p_1^{k_1}} \ C_n^m \equiv a_2 \pmod{p_2^{k_2}} \ \vdots \ C_n^m \equiv a_t \pmod{p_t^{k_t}} \end{cases} ] 中的每一个 ( a_i )那么模 ( P ) 下的唯一解 ( x ) 可以通过标准 CRT 算法求得。因此原问题被简化为一个新的子问题如何计算 ( C_n^m \bmod p^k )这里 ( p ) 是质数( k \ge 1 )。我们只需要集中精力攻克这个子问题。注意在实际编码中CRT 的合并步骤需要计算模 ( P ) 意义下的逆元。由于每对 ( p_i^{k_i} ) 和 ( P / p_i^{k_i} ) 是互质的这是由质因数分解的唯一性保证的所以逆元一定存在可以用扩展欧几里得算法安全求出。2.2 第二层处理模质数幂 ( p^k )现在问题变成了计算 ( C_n^m \frac{n!}{m!(n-m)!} \mod p^k )。直接算的障碍在于分母的阶乘可能包含因子 ( p )导致其在模 ( p^k ) 下没有逆元。exLucas 的精髓在于它并不试图直接求分母的逆元而是选择将分子和分母中的质因子 ( p ) 全部“提取”出来分离计算。我们定义一个新的函数 ( f(n) ) 来表示剔除了所有质因子 ( p ) 后的 ( n! ) 在模 ( p^k ) 下的值。同时设 ( g(n) ) 为 ( n! ) 中所含质因子 ( p ) 的个数。那么阶乘可以表示为 [ n! p^{g(n)} \cdot f(n) ] 于是组合数可以重写为 [ C_n^m \frac{n!}{m!(n-m)!} \frac{p^{g(n)} \cdot f(n)}{p^{g(m)} \cdot f(m) \cdot p^{g(n-m)} \cdot f(n-m)} p^{g(n) - g(m) - g(n-m)} \cdot \frac{f(n)}{f(m) f(n-m)} ] 在这个表达式中指数部分 ( g(n) - g(m) - g(n-m) ) 是一个非负整数。如果它大于等于 ( k )那么整个组合数模 ( p^k ) 就是 0因为分子含有至少 ( k ) 个 ( p ) 因子。如果它小于 ( k )那么分数部分 ( \frac{f(n)}{f(m) f(n-m)} ) 中的所有项 ( f(\cdot) ) 都与 ( p ) 互质因为我们事先剔除了所有 ( p ) 因子因此它在模 ( p^k ) 下存在逆元可以正常计算。至此子问题进一步转化为两个更基础的问题如何快速计算 ( g(n) )即 ( n! ) 中质因子 ( p ) 的个数如何快速计算 ( f(n) )即剔除 ( p ) 因子后 ( n! ) 模 ( p^k ) 的值2.3 第三层递归求解 ( f(n) ) 与 ( g(n) )这是算法中最具技巧性的部分也是“扩展卢卡斯”得名的原因因为它用到了一个类似 Lucas 定理的递归结构。计算 ( g(n) )这是一个经典的数论公式也称作勒让德定理 [ g(n) \left\lfloor \frac{n}{p} \right\rfloor \left\lfloor \frac{n}{p^2} \right\rfloor \left\lfloor \frac{n}{p^3} \right\rfloor \cdots ] 其含义是1 到 n 中有 ( \lfloor n/p \rfloor ) 个数是 p 的倍数贡献至少一个 p 因子有 ( \lfloor n/p^2 \rfloor ) 个数是 ( p^2 ) 的倍数在刚才的基础上再贡献一个 p 因子以此类推。这个求和可以在 ( O(\log_p n) ) 时间内完成。计算 ( f(n) )考虑 ( n! 1 \times 2 \times 3 \times \cdots \times n )。我们可以把这些数按模 ( p ) 的余数进行分组所有不能被 ( p ) 整除的数。所有能被 ( p ) 整除的数。对于能被 ( p ) 整除的数我们可以提取出一个公因子 ( p )剩下的部分就变成了 ( 1, 2, ..., \lfloor n/p \rfloor )。而这正好是 ( \lfloor n/p \rfloor ! )。因此我们有递归关系 [ n! \left( \prod_{i1, p \nmid i}^{p^k} i \right)^{\lfloor n / p^k \rfloor} \times \left( \prod_{i1, p \nmid i}^{n \bmod p^k} i \right) \times \left( \lfloor n/p \rfloor ! \right) ] 这个式子需要仔细理解第一部分( \prod_{i1, p \nmid i}^{p^k} i ) 是一个长度为 ( p^k ) 的周期中所有与 ( p ) 互质的数的乘积模 ( p^k )。因为模 ( p^k ) 下的乘法具有周期性周期为 ( p^k )所以当 ( n ) 很大时这样的完整周期会出现 ( \lfloor n / p^k \rfloor ) 次。第二部分( \prod_{i1, p \nmid i}^{n \bmod p^k} i ) 处理最后一个不完整周期中与 ( p ) 互质的数的乘积。第三部分( \lfloor n/p \rfloor ! ) 来自于将所有 ( p ) 的倍数提取因子 ( p ) 后剩下的部分。注意这里递归计算 ( \lfloor n/p \rfloor ! ) 时我们同样需要剔除其中的 ( p ) 因子所以递归调用的是 ( f(\lfloor n/p \rfloor) )而不是普通的阶乘。然而我们定义 ( f(n) ) 是剔除了所有( p ) 因子后的 ( n! )。在第三部分 ( \lfloor n/p \rfloor ! ) 中仍然可能包含 ( p ) 因子。因此正确的递归式是 [ f(n) \left( \prod_{i1, p \nmid i}^{p^k} i \right)^{\lfloor n / p^k \rfloor} \times \left( \prod_{i1, p \nmid i}^{n \bmod p^k} i \right) \times f\left( \lfloor n/p \rfloor \right) \pmod{p^k} ] 同时对于 ( g(n) ) 的递归式也呼应了这一点 [ g(n) \lfloor n/p \rfloor g(\lfloor n/p \rfloor) ] 这个递归式的边界条件是( f(0) 1 ), ( g(0) 0 )。实操心得第一部分 ( \prod_{i1, p \nmid i}^{p^k} i \mod p^k ) 可以预处理出来记为pre[p^k]。因为对于固定的 ( p^k )这个值是一个常数。这能极大减少递归中重复的计算量。3. 算法实现与关键细节理解了上述三层思路后我们可以勾勒出 exLucas 算法的完整实现步骤。我将结合代码片段和详细注释解释每一个关键环节。3.1 第一步质因数分解模数 P这是整个算法的入口。我们需要将合数模数 ( P ) 分解为vectorpairlong long, int factors其中每个pair存储质因子p和其指数k。vectorpairlong long, int factorize(long long P) { vectorpairlong long, int res; for (long long i 2; i * i P; i) { if (P % i 0) { int cnt 0; while (P % i 0) { P / i; cnt; } res.emplace_back(i, cnt); // 存储 p^k 中的 p 和 k } } if (P 1) { // 处理剩余的大质数 res.emplace_back(P, 1); } return res; }3.2 第二步实现计算 C(n, m) mod p^k 的核心函数这个函数对应我们第二层和第三层的分析。它接收参数n, m, p, k返回C(n, m) mod p^k。// 快速幂取模 long long pow_mod(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; } // 扩展欧几里得求逆元要求 gcd(a, mod) 1 long long inv(long long a, long long mod) { long long x, y; exgcd(a, mod, x, y); // 假设 exgcd 已实现 return (x % mod mod) % mod; } // 计算 f(n) 和 g(n) pairlong long, long long factorial_prime_power(long long n, long long p, long long pk) { if (n 0) return {1, 0}; // 递归计算 long long res 1, cnt 0; // 处理完整周期pre[pk] 的 floor(n/pk) 次方 res res * pow_mod(pre[pk], n / pk, pk) % pk; // 处理不完整周期 for (long long i 1; i n % pk; i) { if (i % p ! 0) { res res * i % pk; } } auto rest factorial_prime_power(n / p, p, pk); res res * rest.first % pk; cnt n / p rest.second; return {res, cnt}; } // 计算 C(n, m) mod p^k long long C_prime_power(long long n, long long m, long long p, long long pk) { if (m n) return 0; 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 p_cnt fn.second - fm.second - fnm.second; if (p_cnt k) return 0; // 组合数模 p^k 为 0 long long res fn.first * inv(fm.first, pk) % pk * inv(fnm.first, pk) % pk; res res * pow_mod(p, p_cnt, pk) % pk; return res; }关键细节解析预处理pre[pk]在程序开始时对于每一个需要用到的pk即p^k我们需要预先计算pre[pk] prod_{i1, p∤i}^{pk} i mod pk。这个计算只需一次存储在全局数组或哈希表中。递归深度factorial_prime_power的递归调用n / p因此递归深度是 ( O(\log_p n) )效率很高。逆元计算在C_prime_power中我们对fm.first和fnm.first求逆元。由于f(x)的定义已经剔除了所有p因子所以fm.first和fnm.first都与p互质进而与pk互质逆元一定存在。3.3 第三步中国剩余定理合并结果现在我们有了计算C(n, m) mod p_i^{k_i}的能力最后一步就是用 CRT 将所有结果合并。// 中国剩余定理 (CRT) 合并方程组 x ≡ a_i (mod m_i)其中 m_i 两两互质 long long crt(const vectorlong long a, const vectorlong long m, long long P) { long long res 0; for (size_t i 0; i a.size(); i) { long long Mi P / m[i]; long long inv_Mi inv(Mi, m[i]); // 求 Mi 在模 m[i] 下的逆元 res (res a[i] * Mi % P * inv_Mi % P) % P; } return res; } // exLucas 主函数 long long exLucas(long long n, long long m, long long P) { if (m n) return 0; vectorpairlong long, int factors factorize(P); vectorlong long a, mods; for (auto [p, k] : factors) { long long pk 1; for (int i 0; i k; i) pk * p; // 预处理 pre[pk]如果尚未计算的话 if (pre.find(pk) pre.end()) { pre[pk] 1; for (long long i 1; i pk; i) { if (i % p ! 0) { pre[pk] pre[pk] * i % pk; } } } long long ai C_prime_power(n, m, p, pk); a.push_back(ai); mods.push_back(pk); } // 使用 CRT 合并所有结果 return crt(a, mods, P); }4. 复杂度分析、优化与边界处理一个完整的算法实现离不开对性能的考量和对特殊情况的处理。4.1 时间复杂度分析假设模数 ( P ) 分解为 ( t ) 个质因数幂最大的 ( p^k ) 为 ( M )。质因数分解( O(\sqrt{P}) )这是主要瓶颈之一但对于竞赛中常见的 ( P \le 10^9 ) 是可行的。预处理pre[pk]对于每个 ( p^k )需要 ( O(p^k) ) 时间。当 ( p ) 很小但 ( k ) 很大时例如 ( 2^{20} )这可能成为瓶颈。但通常 ( p^k ) 作为 ( P ) 的因子其值不会太大。计算C_prime_power递归计算factorial_prime_power的深度为 ( O(\log_p n) )每层递归需要 ( O(p^k) ) 的时间处理不完整周期如果pre[pk]已预处理则这部分是 ( O(1) )。因此对于单个 ( p^k )复杂度约为 ( O(p^k \log_p n) )。CRT合并( O(t) )可忽略。总复杂度大致为 ( O(\sqrt{P} \sum (p_i^{k_i} \log_{p_i} n)) )。在实战中最需要关注的是p^k的大小。如果 ( P ) 包含一个像 ( 2^{30} ) 这样的因子预处理pre[pk]将需要计算超过10亿个数的乘积这显然是不可接受的。4.2 针对大 ( p^k ) 的优化技巧当 ( p^k ) 较大时直接计算pre[pk]会超时。我们需要一个更聪明的办法来计算f(n)。观察递归式 [ f(n) \left( \prod_{i1, p \nmid i}^{p^k} i \right)^{\lfloor n / p^k \rfloor} \times \left( \prod_{i1, p \nmid i}^{n \bmod p^k} i \right) \times f\left( \lfloor n/p \rfloor \right) \pmod{p^k} ] 关键在于第一部分 ( \prod_{i1, p \nmid i}^{p^k} i )。在模 ( p^k ) 下这个乘积有一个重要性质它等于 ( -1 \mod p^k )当 ( p2 ) 或 ( k1 ) 时。这是因为在缩系与 ( p^k ) 互质的数构成的集合中每个数都有唯一的逆元且除了1和-1即 ( p^k-1 )之外其他数都是成对互为逆元出现的。对于 ( p2, k1 ) 的情况这个乘积等于1。因此我们可以优化当 ( p 2 ) 或 ( k 1 ) 时pre[pk] pk - 1。当 ( p 2 ) 且 ( k 1 ) 时pre[pk] 1。这个优化将预处理复杂度从 ( O(p^k) ) 降到了 ( O(1) )是 exLucas 算法能够处理较大模数的关键。修正后的factorial_prime_power函数中计算完整周期部分变为long long cycle_prod (p 2 k 1) ? 1 : (pk - 1); res res * pow_mod(cycle_prod, n / pk, pk) % pk;4.3 边界条件与常见陷阱m n或m 0组合数定义中此时值为0。应在函数入口处检查。P 1任何数模1都为0。可以特判。p^k可能溢出在计算pk p^k时尤其是p较小而k较大时可能会超出long long范围。在竞赛中通常保证 ( P \le 10^9 )所以pk也在可控范围内。但为了鲁棒性可以添加溢出检查。递归基factorial_prime_power(0, p, pk)必须返回{1, 0}。逆元不存在在基本算法流程中我们只在确保互质的情况下求逆元f(m)部分。但在实现 CRT 合并时求inv(Mi, m[i])也必须保证Mi与m[i]互质这由Mi P / m[i]和m[i]是P的质因数幂次保证是成立的。5. 实战演练与问题排查理论再完美也需要代码实现来验证。让我们通过一个具体例子并讨论调试中常见的问题。5.1 示例计算C(10, 3) mod 12分解模数( P 12 2^2 \times 3^1 )。所以我们需要计算C(10,3) mod 4和C(10,3) mod 3。计算 mod 4 (p2, k2, pk4)计算C(10,3) mod 4。n10, m3, n-m7。g(10) floor(10/2)floor(10/4)floor(10/8)5218。g(3)1, g(7)3。p_cnt 8-1-34。因为4 k2所以C(10,3) mod 4 0。验证C(10,3)120120 / 4 30余0。正确。计算 mod 3 (p3, k1, pk3)计算C(10,3) mod 3。g(10)3, g(3)1, g(7)2。p_cnt 3-1-20。计算f(10) mod 310!中剔除3的因子。f(10) (1*2*4*5*7*8*10) * f(3) mod 3。在模3下1*2*4*5*7*8*10 ≡ 1*2*1*2*1*2*1 ≡ 8 ≡ 2 mod 3。f(3) (1*2) * f(1) 2 mod 3。所以f(10) 2*24 ≡ 1 mod 3。同理f(3)2 mod 3,f(7)1*2*4*5*7 (剔除非3倍数) 1*2*1*2*1 ≡ 4 ≡ 1 mod 3。C 1 * inv(2,3) * inv(1,3) mod 3 1 * 2 * 1 2 mod 3。验证C(10,3)120120 mod 3 0等等这里出错了我们计算结果是2 mod 3但实际120 mod 3 0。问题出在哪检查g值g(10)floor(10/3)floor(10/9)314。g(3)1,g(7)floor(7/3)2。p_cnt4-1-21。啊哈这里p_cnt1而k1所以p_cnt k成立组合数模3^1应该为0我最初计算g(10)时漏掉了floor(10/9)项。修正后C(10,3) mod 3 0。CRT 合并我们有方程组x ≡ 0 (mod 4),x ≡ 0 (mod 3)。解显然是x ≡ 0 (mod 12)。最终结果C(10,3) mod 12 0。验证120 mod 12 0。正确。这个例子揭示了手动计算时容易出现的错误勒让德公式求g(n)时求和项必须持续到p^t n。5.2 常见问题与调试技巧结果错误特别是为0的情况首先检查g(n)的计算这是最容易出错的地方。确保递归或循环正确计算了所有floor(n / p^i)项。检查p_cnt k的判断如果p_cnt即g(n)-g(m)-g(n-m)大于等于k结果直接为0。很多错误源于漏判。验证逆元计算确保在计算f(n) / (f(m)*f(n-m))时分母的逆元存在即f(m)和f(n-m)与p互质。如果算法逻辑正确这应该自然满足。性能问题超时未使用-1 mod p^k优化这是最大的性能陷阱。务必对p2 or k1的情况将完整周期乘积设为pk-1。重复预处理对相同的(p, k)组合pre[pk]只需计算一次用哈希表或数组缓存起来。递归开销factorial_prime_power是尾递归可以很容易改为迭代形式以减少函数调用开销。数值溢出中间乘法溢出在计算res res * i % pk或类似操作时即使res和i都小于pk它们的乘积也可能超过long long范围如果pk接近10^9乘积可能约10^18。需要使用快速乘或__int128来处理。// 使用 __int128 的安全乘法取模 long long mul_mod(long long a, long long b, long long mod) { return (long long)((__int128)a * b % mod); }pk本身溢出在计算pk p^k时使用循环累乘并加入溢出判断。特殊模数P是质数此时直接使用 Lucas 定理或预处理阶乘逆元会更简单高效。可以在 exLucas 入口处添加判断如果P是质数则路由到更快的算法。P很大但质因子很少exLucas 表现良好。如果P很大且质因子很多比如 square-free number每个p^k p都很小但数量t很多CRT 合并步骤的复杂度是线性的可以接受。在实际比赛中实现 exLucas 需要格外细心。建议从一个简单、未优化的版本开始确保逻辑正确然后再逐步加入优化如-1优化、缓存。编写针对性的测试用例比如测试C(n,0),C(n,n),C(n,1)以及模数较小的情况并与暴力计算对于小的 n, m的结果进行对比是确保代码正确的有效方法。记住数论算法的调试往往需要你拿起纸笔重新演算一遍中间步骤。