1. 从组合数计算到模运算Lucas定理的引入在算法竞赛和数论编程中计算组合数 C(n, m) 是一个基础且高频的需求。当 n 和 m 的数值较小时我们可以直接利用公式 C(n, m) n! / (m! * (n-m)!) 进行计算或者通过递推式 C(n, m) C(n-1, m-1) C(n-1, m) 来构建杨辉三角。然而当 n 和 m 非常大比如 10^18时这两种方法都会失效公式计算涉及超大整数的阶乘远超任何数据类型的表示范围递推则需要 O(n^2) 的时间和空间完全不可行。更常见且实际的需求是计算 C(n, m) mod p 的值其中 p 是一个素数。这就是 Lucas 定理要解决的核心问题。Lucas 定理提供了一个极其高效的方法能将大规模组合数取模问题分解为一系列小规模组合数取模问题的乘积。它的表述非常优雅对于素数 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) ≡ Π C(n_i, m_i) (mod p) 其中C(n_i, m_i) 是普通的小组合数当 m_i n_i 时定义 C(n_i, m_i) 0。我第一次接触这个定理时感觉像是发现了一个“降维打击”的武器。原本需要处理天文数字 n 的组合数现在只需要处理一系列最多为 p-1 的小数字 n_i 和 m_i 的组合数。计算复杂度从 O(n) 骤降到 O(log_p n)这对于 n 高达 10^18 而 p 为 10^5 量级的情况来说是质的飞跃。1.1 Lucas定理的证明思路与核心逻辑理解一个定理最好的方式就是搞懂它为什么成立。Lucas定理的证明巧妙利用了生成函数和模运算的性质。一个经典的证明思路如下考虑多项式 (1x)^n 在模 p 意义下的展开。一方面根据二项式定理(1x)^n Σ C(n, k) * x^k。 另一方面我们可以用p进制表示 n并利用模p运算的性质(1x)^n (1x)^{n_0} * ((1x)^p)^{n_1} * ... * ((1x)^{p^k})^{n_k}。这里有一个关键步骤对于素数 p有 (1x)^p ≡ 1 x^p (mod p)。这是因为根据二项式定理(1x)^p 展开后除了 C(p,0) 和 C(p,p) 项系数为1中间项系数 C(p, k) 都包含因子 p因此在模 p 意义下均为0。这个性质是费马小定理的一个多项式版本。于是((1x)^p)^{n_i} ≡ (1 x^p)^{n_i} (mod p)。最终(1x)^n 在模 p 下可以写成 (1x)^{n_0} * (1x^p)^{n_1} * ... * (1x^{p^k})^{n_k}。现在我们想从这个乘积的展开式中找出 x^m 项的系数。注意 m 也有p进制表示 m m_0 m_1p ... m_kp^k。要从上述乘积中得到 x^m唯一的方式是从第一个因子取 x^{m_0} 项系数为 C(n_0, m_0)从第二个因子取 (x^p)^{m_1} 项系数为 C(n_1, m_1)依此类推。将这些系数相乘就得到了 x^m 的总系数即 C(n, m)。由于每一步都是在模 p 意义下进行的所以最终有 C(n, m) ≡ Π C(n_i, m_i) (mod p)。这个证明过程清晰地揭示了定理的本质它将一个“整体”的模运算分解到了每个“数位”上独立进行这正是p进制表示带来的威力。注意在实现时我们通常直接对 n 和 m 循环除以 p并同步计算每一数位上的小组合数。当出现某一数位上 m_i n_i 时根据定义 C(n_i, m_i)0这意味着整个大组合数 C(n, m) 模 p 的结果就是 0。这是一个非常重要的优化和边界判断。1.2 标准Lucas算法的实现与细节理解了原理实现就是一个清晰的递归或迭代过程。以下是 Lucas 算法的典型实现步骤预处理阶乘与逆元由于我们需要反复计算 C(n_i, m_i) mod p而 n_i, m_i p我们可以预处理出 0 到 p-1 所有数的阶乘fac[i]模 p 的值以及对应的阶乘逆元invfac[i]。这样小组合数可以通过fac[n] * invfac[m] % p * invfac[n-m] % p在 O(1) 时间内得到。计算逆元通常使用费马小定理当 p 为素数时a^(p-2) % p即为 a 的逆元或扩展欧几里得算法。分解数位并计算循环进行当 n0 或 m0 时取当前数位ni n % p,mi m % p。判断如果mi ni直接返回 0。否则计算小组合数C(ni, mi)利用预处理好的阶乘和逆元。将结果乘入总答案ans中并对 p 取模。n 和 m 同时除以 pn / p; m / p;。返回最终结果。这里有一个实操心得预处理阶乘数组时务必把fac[0]和invfac[0]设为 1。因为在计算组合数时当 m0 或 nm 时会用到 0!。很多初学者在这里会出错。另一个常见问题是模数 p 可能很小比如 2, 3, 5但 n 和 m 极大。这时Lucas 定理的迭代次数log_p n会很多但每一步的计算量计算小组合数是 O(1) 的所以整体效率依然非常高。但是如果 p 非常大接近 n那么定理的分解优势就不明显了不过此时 n_i 和 m_i 就是 n 和 m 本身算法退化为直接计算一次组合数也是正确的。代码框架示例C风格// 预处理阶乘和逆元 (0!...p-1! mod p) vectorlong long fac(p), invfac(p); fac[0] 1; for (int i 1; i p; i) fac[i] fac[i-1] * i % p; invfac[p-1] pow_mod(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 auto C_small [](int a, int b) - long long { if (b a) return 0; return fac[a] * invfac[b] % p * invfac[a-b] % p; }; // Lucas 定理主函数 functionlong long(long long, long long) lucas [](long long n, long long m) - long long { if (m 0) return 1; return C_small(n % p, m % p) * lucas(n / p, m / p) % p; };2. 当模数非素数ExLucas扩展卢卡斯的挑战与方案Lucas 定理的强大建立在模数 p 是素数这一前提上。因为证明中用到了(1x)^p ≡ 1x^p (mod p)这一关键性质而这仅对素数 p 成立。在实际问题中尤其是密码学或某些特殊的计数问题里模数可能不是素数而是一个合数。例如可能需要计算 C(n, m) mod M其中 M 10007是素数可用Lucas但也可能是 M 1000000007是素数或者是 M 99991是素数但更棘手的是 M 1000, 2024, 或者 p^kp是素数k1如 2^10, 5^3 等。当模数不是素数时标准 Lucas 定理直接失效。我们不能直接求阶乘的逆元因为合数模下并非所有数都有乘法逆元。例如在模 4 意义下2 就没有逆元。这时我们就需要更强大的工具——ExLucas扩展卢卡斯定理。ExLucas 的核心思想是“分解与合并”运用了中国剩余定理CRT的思想。它的目标是将问题C(n, m) mod M分解为若干个模数为p_i^{k_i}的子问题其中M Π p_i^{k_i}是 M 的质因数分解。最后再将这些子问题的解用中国剩余定理合并回模 M 的解。2.1 问题拆解核心矛盾与解决思路所以ExLucas 要解决的根本问题是如何计算 C(n, m) mod p^k其中 p 是素数k 是正整数。计算 C(n, m) mod p^k 的难点在于分母m!和(n-m)!可能包含因子 p导致它们与模数 p^k 不互质从而无法直接求逆元。例如计算 C(10, 4) mod 2^3。C(10,4)210模8余2。但如果我们试图用阶乘公式10! / (4! * 6!)在模8下计算4! 24 模8为0根本没有逆元。解决思路是将阶乘中所有因子 p 剥离出来单独计算。具体来说对于任意一个数 x 的阶乘 x!我们可以把它写成x! p^a * b其中a是 x! 中质因子 p 的个数b是一个与 p 互质的整数即b是剔除了所有 p 因子后剩余部分的乘积模 p^k。那么组合数可以表示为C(n, m) n! / (m! * (n-m)!) p^{a_n - a_m - a_{n-m}} * (b_n / (b_m * b_{n-m}))在模 p^k 的意义下p^{a_n - a_m - a_{n-m}}这部分如果指数非负它就贡献了一个 p 的幂次如果指数为负理论上整个组合数应该是个分数但在模运算中我们通常只处理整数所以这种情况分母中p的幂次多于分子意味着结果不能被 p^k 整除但模 p^k 后可能不为0需要更精细的处理。实际上在计算模 p^k 时我们最终需要的是b_n * (b_m)^{-1} * (b_{n-m})^{-1} * p^{delta}模 p^k 的值其中 delta a_n - a_m - a_{n-m}而逆元运算(b)^{-1}是在与 p 互质的意义下进行的所以是可行的。因此问题的关键转化为两个子函数get_factor(x, p, pk)计算 x! 中剔除所有 p 因子后的部分b模 pk 的值以及因子 p 的个数a。C_mod_pk(n, m, p, k)利用上面的函数计算 C(n, m) mod p^k。2.2 计算剥离p因子后的阶乘get_factor 函数详解如何高效计算x!中剔除所有 p 因子后的乘积b模p^k呢我们不能直接计算 x! 再除因为 x 可能很大。这里需要一个巧妙的递归或循环方法。观察 x! 的乘积1 * 2 * 3 * ... * x。我们可以把这些数按模 p 的余数分组p 的倍数p, 2p, 3p, ..., floor(x/p) * p。把这些 p 提出来总共能提出floor(x/p)个 p。剩下的部分是1, 2, 3, ..., floor(x/p)的乘积这恰好是floor(x/p)!。于是这部分贡献了p^{floor(x/p)} * (floor(x/p)!)但floor(x/p)!里面可能还包含 p 的因子所以需要递归处理。非 p 的倍数这些数可以每p^k个一组因为模p^k下这些数的乘积具有周期性。具体来说在1..p^k这个周期内所有与 p 互质的数的乘积模p^k是一个常数记作prod。这样的完整周期有floor(x / p^k)个每个周期贡献prod^{floor(x / p^k)} mod p^k。最后不完整的周期剩下1..(x mod p^k)中与 p 互质的数需要单独计算乘积。此外还需要考虑每组中“非p倍数”的部分。更常见的优化写法是循环处理 x每次将 x! 分解为三部分的乘积第一部分所有不超过 x 且是 p 的倍数的数提出一个 p 后转化为(x/p)!的问题递归。第二部分在1..p^k这个块内所有与 p 互质的数的乘积记为pre[pk]可以预处理。这样的块有x / pk个。第三部分最后一个不完整的块1..(x % pk)中与 p 互质的数的乘积。同时因子 p 的个数a可以通过不断累加x / p,x / p^2, ... 来得到。实现细节与注意事项预处理pre[i]数组pre[i]表示1..i中所有与 p 互质的数的乘积模pk。这可以在 O(pk) 时间内完成。递归函数get_factor(x, p, pk)返回一个 pair(b, a)其中b是剔除p因子后的乘积模 pka是p因子的总个数。递归边界当x 0时返回(1, 0)。递归式(b, a) get_factor(x/p, p, pk)。然后当前层的b需要乘以pre[pk]^{x/pk} * pre[x % pk]并对 pk 取模。当前层的a需要加上x/p。这里有一个极易出错的地方计算pre[pk]^{x/pk}时指数x/pk可能很大需要使用快速幂算法。同时所有乘法运算都要随时对pk取模。代码思路示例// 预处理 pre 数组 pre[j] (1*2*...*j) mod pk, 其中跳过p的倍数 vectorlong long pre(pk 1); pre[0] 1; for (int i 1; i pk; i) { if (i % p 0) pre[i] pre[i-1]; // 跳过p的倍数 else pre[i] pre[i-1] * i % pk; } // 计算 x! 中剔除p因子后的部分b和p的指数a pairlong long, long long get_factor(long long x, long long p, long long pk) { if (x 0) return {1, 0}; auto [b_next, a_next] get_factor(x / p, p, pk); long long b b_next * pow_mod(pre[pk], x / pk, pk) % pk * pre[x % pk] % pk; long long a a_next x / p; return {b, a}; }3. 整合与求解C_mod_pk 与 CRT 合并有了get_factor函数计算C(n, m) mod p^k就变得直接了。3.1 计算单个质数幂模数下的组合数函数C_mod_pk(n, m, p, k)的步骤如下计算pk p^k。分别调用get_factor得到(b_n, a_n),(b_m, a_m),(b_nm, a_nm)对应 n!, m!, (n-m)!。计算指数差delta_a a_n - a_m - a_nm。如果delta_a k那么说明分子中 p 的幂次至少比分母多 k 次整个组合数能被p^k整除因此C(n, m) mod p^k 0。否则我们需要计算b b_n * inv(b_m, pk) * inv(b_nm, pk) % pk。注意这里的b_m和b_nm都是与 p 互质的因为我们在get_factor中已经剔除了所有 p 因子所以它们在模pk下存在逆元。求逆元可以用扩展欧几里得算法因为 gcd(b_m, pk) 1。最终结果ans b * pow_mod(p, delta_a, pk) % pk。因为delta_a可能很大但pk是p^k所以p^{delta_a}模pk可以直接计算实际上就是p^{delta_a}本身只要delta_a k或者如果delta_a不小可以用快速幂。这里有一个关键技巧在计算逆元时由于模数pk可能不是素数我们不能用费马小定理必须使用扩展欧几里得算法exgcd来求解。实现示例long long C_mod_pk(long long n, long long m, long long p, long long k) { long long pk 1; for (int i 0; i k; i) pk * p; // 计算 p^k if (pk 1) return 0; // 特殊情况处理 auto [bn, an] get_factor(n, p, pk); auto [bm, am] get_factor(m, p, pk); auto [bnm, anm] get_factor(n - m, p, pk); long long delta_a an - am - anm; if (delta_a k) return 0; // 结果能被 p^k 整除 long long b bn * inv_exgcd(bm, pk) % pk * inv_exgcd(bnm, pk) % pk; long long ans b * pow_mod(p, delta_a, pk) % pk; return ans; }3.2 使用中国剩余定理CRT合并答案现在我们已经能计算C(n, m) mod p_i^{k_i}对于每个质因数p_i^{k_i}的值记作a_i。我们的模数M Π p_i^{k_i}且这些p_i^{k_i}两两互质因为它们是不同质数的幂。根据中国剩余定理存在唯一解x mod M满足方程组x ≡ a_i (mod p_i^{k_i})对所有 i 成立。中国剩余定理的经典求解方法对于两两互质的模数是 设M_i M / m_i其中m_i p_i^{k_i}计算M_i在模m_i下的逆元t_i即M_i * t_i ≡ 1 (mod m_i)。则解为x ≡ Σ (a_i * M_i * t_i) (mod M)。在实现中我们通常依次合并两个同余方程。假设当前已经合并得到解x满足前几个方程模数为M_cur当前解为ans。现在要加入一个新的方程x ≡ a (mod m)其中m与M_cur互质。 我们需要找到一个新的解x满足x ≡ ans (mod M_cur)且x ≡ a (mod m)。 这等价于寻找整数k使得ans k * M_cur ≡ a (mod m)。 解这个关于k的线性同余方程k * M_cur ≡ a - ans (mod m)。 由于gcd(M_cur, m) 1M_cur在模m下有逆元可以用 exgcd 求出inv_M。 则k ≡ (a - ans) * inv_M (mod m)。 取k为最小非负解则新的解x ans k * M_cur新的模数M_new M_cur * m。实操中的注意事项在计算(a - ans)时由于是在模m下运算需要先对a - ans取模m防止负数。乘法运算(a - ans) * inv_M可能很大需要随时取模m。最终合并得到的ans可能很大但一定在[0, M_new)范围内。整个合并过程的时间复杂度是 O(l * log M)其中 l 是质因数的个数通常很小。合并代码示例// 使用扩展欧几里得求逆元 long long inv_exgcd(long long a, long long mod) { ... } // 合并同余方程组方程组表示为 vectorpairlong long, long long eq, 其中 .first 是模数 m, .second 是余数 a long long CRT(const vectorpairlong long, long long eq) { long long ans 0, M 1; for (auto [m, a] : eq) { // 求解 k * M ≡ a - ans (mod m) long long delta (a - ans % m m) % m; // 确保非负 long long gcd_val, inv_M, y; exgcd(M, m, gcd_val, inv_M, y); // 解 M * inv_M m * y gcd(M,m)1 inv_M (inv_M % m m) % m; // 得到 M 模 m 的逆元 long long k delta * inv_M % m; ans ans k * M; M M * m; ans % M; } return ans; }最后ExLucas 的主函数就是对模数 M 进行质因数分解对每个质因子幂调用C_mod_pk得到余数然后用 CRT 合并。4. ExLucas的实战技巧、优化与边界问题虽然 ExLucas 的原理和步骤已经清晰但在实际编码和竞赛中仍有不少细节和技巧需要掌握否则极易写出效率低下或错误的代码。4.1 预处理优化与时间复杂度分析ExLucas 的瓶颈主要在于get_factor函数中的递归和pre数组的计算。假设我们需要处理模数 M其一个质因子为 p对应的幂为 pk p^k。pre数组预处理需要计算长度为 pk 的数组复杂度 O(pk)。当 pk 很大时比如 p2, k30pk 超过10亿这是不可能的。但请注意在get_factor中我们只用到pre[pk]和pre[x % pk]。pre[pk]是一个固定值即1..pk中所有与 p 互质的数的乘积模 pk。这个值可以单独计算而不需要整个数组。计算pre[pk]可以用循环但更聪明的方法是注意到这个乘积具有周期性实际上pre[pk]在模 pk 下常常等于-1或某个特定值根据威尔逊定理的推广但为了通用性我们通常还是用循环计算但只算一次复杂度 O(pk)。对于pre[x % pk]x % pk 小于 pk我们可以实时计算也可以在递归时传入一个预处理的、长度不超过 pk 的pre数组如果 pk 不大。如果 pk 非常大比如 10^7那么 ExLucas 算法本身可能就不太适用了需要考虑其他方法或题目本身的特殊性。递归深度get_factor递归深度为O(log_p x)每次递归需要一次快速幂计算pre[pk]^{x/pk}和一些乘法。对于 n, m 高达 10^18p 很小的情况递归深度可能达到 60但每次递归的运算是常数时间快速幂是 O(log(x/pk))所以单次get_factor调用可以认为是 O(log n) 级别的。总体复杂度对于模数 M 的每个质因子 p^k需要计算三次get_factor和若干次逆元、乘法。设 M 的质因子个数为 ω(M)那么总复杂度大致为 O(ω(M) * (pk log n))。当 pk 很大时O(pk) 是主要开销。优化建议小模数预处理如果题目中模数 M 是固定的并且其质因子幂p^k都不大比如 pk 10^7可以预先为每个质因子幂计算好pre[pk]的值甚至预处理整个pre数组如果内存允许这样每次查询就是 O(log n) 的。避免重复计算在C_mod_pk中我们计算了 n, m, n-m 的get_factor。注意n-m可能小于 m但无法避免三次调用。使用迭代而非递归get_factor可以用 while 循环实现避免递归栈的开销但递归写法通常更清晰。4.2 常见错误与调试技巧实现 ExLucas 时以下几个坑点几乎每个初学者都会遇到逆元计算错误在模p^k下求逆元必须使用扩展欧几里得算法exgcd因为模数不是素数。绝对不要使用费马小定理pow_mod(a, pk-2, pk)这仅在 pk 是素数时成立。整数溢出中间计算涉及多次乘法即使最终结果在模数范围内中间过程a * b % pk也可能溢出 64 位整数如果 a 和 b 都接近 10^18pk 也很大。需要使用快速乘龟速乘或__int128来处理乘法取模。这是非常重要的一个技巧。// 快速乘 (防止溢出) long long mul_mod(long long a, long long b, long long mod) { long long res 0; a % mod; while (b) { if (b 1) res (res a) % mod; a (a a) % mod; b 1; } return res; } // 或者在支持的环境下直接使用 __int128 long long mul_mod(long long a, long long b, long long mod) { return (long long)((__int128)a * b % mod); }负数取模在 CRT 合并步骤计算delta (a - ans) % m时必须确保结果为非负。标准的写法是(a % m - ans % m m) % m。delta_a k的判断在C_mod_pk中如果delta_a k直接返回 0。这个判断必须做否则后续计算p^{delta_a}模p^k会得到 0但乘法顺序可能导致错误。更严谨地说当delta_a k时C(n, m)可被p^k整除模p^k的结果就是 0。get_factor中pre[pk]的计算pre[pk]是1..pk中所有与 p 互质的数的乘积模 pk。注意在计算这个乘积时如果遇到 p 的倍数就跳过。循环到 ipk 时i 就是 pk 本身它一定包含因子 p所以应该跳过。因此pre[pk]实际上是1..pk-1中与 p 互质的数的乘积模 pk。很多实现中会预处理pre[0..pk]但让pre[pk] pre[pk-1]这样在快速幂pow_mod(pre[pk], x/pk, pk)时才是正确的。4.3 典型问题场景与ExLucas的适用性ExLucas 虽然强大但并非万能。它的时间复杂度与模数 M 的质因子分解中最大的p^k密切相关。如果 M 有一个很大的质数 p那么 pk p算法效率很高因为pre数组长度是 p。但如果 M 包含一个像2^30这样的质数幂那么 pk 2^30 ≈ 10^9预处理pre数组或计算pre[pk]将非常耗时甚至不可能。因此在遇到需要计算大组合数模合数的问题时首先要分析模数 M 的特征M 是素数直接用普通 Lucas 定理简单高效。M 是无平方因子数square-free即 M 是若干不同素数的乘积。此时每个质因子幂p^k中的 k1pkp。ExLucas 对于每个质因子只需要处理模 p 的情况这时get_factor函数中的pre计算很简单因为 pkppre[p]就是 1到p-1所有数的乘积模 p根据威尔逊定理等于 -1 即 p-1。实际上这退化为了对每个质因子用普通 Lucas 定理求值再用 CRT 合并。效率很高。M 包含高次质数幂例如M 2^60。这时 ExLucas 的pk极大算法无法直接应用。这类问题通常需要利用组合数自身的数论性质或者题目有特殊限制如 n, m 较小可能需要其他方法。一个实用的建议在竞赛中如果模数 M 不是素数并且 M 的大小在 10^6 以内可以尝试用 ExLucas。如果 M 很大但 n, m 相对较小比如 10^7 以内可以考虑直接预处理阶乘和逆元但需要处理模数非素数时逆元不存在的问题这又回到了 ExLucas 要解决的问题。有时题目设计的模数会保证其质因子幂次不大使得 ExLucas 可行。最后分享一个我调试 ExLucas 代码时的技巧先写一个暴力计算小范围组合数取模的程序用于验证。然后针对小的合数模数如 6, 10, 12测试 ExLucas 的输出是否正确。再逐步测试包含质数幂的模数如 8, 9。确保小数据完全正确后再尝试大数据。由于算法涉及递归、快速幂、exgcd、CRT 等多个模块任何一个环节出错都可能导致结果不对耐心分模块测试是关键。