1. 从一道“超纲”的组合数取模题说起几年前我在准备一场算法竞赛时遇到了一道让我卡壳很久的题目。题目本身描述很简单给定三个巨大的整数n,m,p要求计算组合数C(n, m) mod p的值。这看起来是组合数学的入门题我熟练地写下了预处理阶乘和逆元的代码信心满满地提交。结果返回了一个大大的“Wrong Answer”。我检查了无数遍代码逻辑确信无误。直到我重新审视数据范围才发现p并不是一个质数而是一个任意的正整数比如p 10007一个合数。那一刻我才明白常规的费马小定理求逆元、预处理阶乘的方法其前提是模数p必须为质数。当p是合数时许多数在模p下可能没有乘法逆元整个算法的基础就崩塌了。这个问题在算法竞赛和密码学等领域其实非常典型。我们常常需要计算大组合数对一个非质数取模的结果。比如在计算一些概率问题时模数可能是一个方便的大整数如10^97是质数但10^96就不是又或者在一些加密协议中模数直接就是两个大质数的乘积。这时扩展卢卡斯定理就成了解决这类问题的“瑞士军刀”。它并不是一个全新的、独立的定理而是经典卢卡斯定理在模数为合数特别是质数幂情况下的自然延伸和组合应用。今天我就结合自己多次实现和调试的经验把扩展卢卡斯定理的核心思想、推导过程、实现细节以及那些容易踩的坑系统地梳理一遍。2. 问题拆解为什么常规方法会失效要理解扩展卢卡斯定理必须先弄清楚我们面临的障碍是什么。为什么模数是质数时一切顺利换成合数就出问题2.1 质数模下的“舒适区”逆元与阶乘当模数p是一个质数时数论为我们提供了强大的工具。根据费马小定理对于任意不被p整除的整数a都有a^(p-1) ≡ 1 (mod p)。这意味着a * a^(p-2) ≡ 1 (mod p)所以a^(p-2)就是a在模p下的乘法逆元。计算组合数C(n, m) n! / (m! * (n-m)!)时我们可以预处理出1!到n!模p的值然后通过乘以分母的逆元来完成“除法”运算。整个过程高效且优雅。2.2 合数模下的“雷区”逆元不存在一旦p是合数情况急转直下。一个数a在模p下存在乘法逆元的充要条件是gcd(a, p) 1即a与p互质。在计算组合数的分母m!或(n-m)!时它们很可能包含p的质因子。例如p6m3那么3! 6gcd(6, 6)6≠1所以3!在模6下没有逆元我们无法直接通过乘以逆元的方式来计算n! / (m! * (n-m)!) mod 6。因此核心矛盾在于组合数公式中的除法在模运算下需要转化为乘逆元而当分母与模数不互质时逆元不存在直接计算路径被阻断。2.3 扩展卢卡斯的核心思路分解与重组扩展卢卡斯定理的智慧在于“分而治之”。既然整个模数p让我们无法求逆我们就把p分解成若干个质数幂的乘积p p1^k1 * p2^k2 * ... * pt^kt。根据中国剩余定理如果我们能分别计算出C(n, m) mod pi^ki的值就可以唯一地确定C(n, m) mod p的值。于是问题转化为如何计算C(n, m) mod q^k其中q是质数k是正整数这才是扩展卢卡斯定理要解决的核心子问题。即使在这个子问题下分母与模数q^k仍然可能不互质因为分母的阶乘里含有质因子q但情况已经大大简化我们可以专门处理质因子q。3. 攻坚核心子问题计算 C(n, m) mod q^k计算C(n, m) mod q^k是扩展卢卡斯最精妙的部分。我们无法直接做除法那就把分子和分母中的质因子q都“提取”出来分开处理。定义n!的“提取”形式为n! q^e * u其中u是一个与q互质的整数e是n!中质因子q的个数。这个e有一个经典公式可以快速计算e floor(n/q) floor(n/q^2) floor(n/q^3) ...。那么组合数可以写成C(n, m) n! / (m! * (n-m)!) (q^{e_n} * u_n) / (q^{e_m} * u_m * q^{e_{n-m}} * u_{n-m}) q^{e_n - e_m - e_{n-m}} * (u_n * (u_m)^{-1} * (u_{n-m})^{-1})这里(u_m)^{-1}和(u_{n-m})^{-1}是在模q^k意义下的逆元。由于u_m和u_{n-m}已经与q互质因为我们把所有的q都提出来了所以它们在模q^k下必然存在逆元这就巧妙地绕开了逆元不存在的障碍。于是计算C(n, m) mod q^k就分解为三个步骤计算指数部分E e_n - e_m - e_{n-m}。计算互质部分U u_n * inv(u_m) * inv(u_{n-m}) mod q^k。最终结果C(n, m) ≡ q^E * U (mod q^k)。注意如果E k那么q^E模q^k就是0组合数模q^k就是0。否则我们需要计算q^E * U mod q^k。现在问题的关键变成了如何高效地计算u_n即n!中剔除所有质因子q后再对q^k取模的结果3.1 快速计算 u_n观察规律与分组直接计算n!再剔除q的因子效率太低。我们需要一个更聪明的办法。考虑n22, q3, k2我们要计算22! mod 9并剔除所有因子3。我们可以把1*2*3*...*22分成三类数相乘q的倍数3, 6, 9, 12, 15, 18, 21。每个数提取一个因子q后变成1, 2, 3, 4, 5, 6, 7。这恰好是floor(22/3)7的阶乘即7!。并且这个过程可以递归进行因为7!里可能还包含q的倍数。与q互质的数1, 2, 4, 5, 7, 8, 10, 11, 13, 14, 16, 17, 19, 20, 22。这些数可以直接参与模q^k的乘法。“尾巴”部分在n不是q^k的整数倍时会有一段不完整的周期。在这个例子里22除以9的余数部分即22 % 9 4对应的数1, 2, 3, 4需要特殊处理其中3要剔除因子3。更一般地我们可以总结出计算F(n, q, k)即u_n mod q^k的递归公式F(n, q, k) [F(floor(n/q), q, k) * (循环节乘积)^{floor(n / q^k)} * (尾巴部分乘积)] mod q^k其中F(floor(n/q), q, k)递归处理所有q的倍数提取因子后形成的阶乘。循环节乘积考虑一个长度为q^k的周期[1, 2, ..., q^k]。在这个周期里所有与q互质的数的乘积记为prod。因为模q^k下的乘法具有周期性所以每q^k个数与q互质的部分的乘积都是prod mod q^k。我们可以预处理这个prod。尾巴部分乘积最后不足一个完整周期的n % q^k个数需要计算它们中与q互质的部分的乘积并剔除其中q的因子虽然尾巴里通常不会有完整的q因子但严谨起见需要处理。通过这个递归过程我们可以在O(log_{q} n)的时间内计算出u_n mod q^k。实操心得预处理“循环节乘积”prod是优化的关键。对于固定的q和kprod只需要计算一次。计算时要注意prod是1到q^k中所有与q互质的数的乘积对q^k取模。由于q^k可能很大直接连乘可能会溢出需要在乘法过程中不断取模。3.2 一个完整的子问题计算示例假设我们要计算C(22, 8) mod 9即q3, k2。计算指数 Ee_22 floor(22/3)floor(22/9)729e_8 floor(8/3)floor(8/9)202e_14 floor(14/3)floor(14/9)415E 9 - 2 - 5 2。因为E2等于k2所以3^2 mod 9 0这里有个关键点q^E * U mod q^k当E k时结果为0。但我们需要先计算U因为最终是(q^E mod q^k) * U mod q^k。3^2 mod 9 0所以无论U是多少结果都是0。但让我们继续完成计算来验证流程。计算互质部分 U需要计算u_22, u_8, u_14 mod 9。预处理循环节乘积prod对于模数9一个周期是[1,2,3,4,5,6,7,8,9]剔除3, 6, 9q的倍数剩下1,2,4,5,7,8。prod 1*2*4*5*7*8 mod 9。计算1*22,2*48,8*540 mod 94,4*728 mod 91,1*88。所以prod 8 mod 9。计算F(22, 3, 2)n22。floor(22/3)7进入递归F(7,3,2)。n7。floor(7/3)2进入递归F(2,3,2)。n2。floor(2/3)0递归基础情况F(0, q, k) 1。回溯计算F(2,3,2)周期数cyc floor(2/9)0尾巴tail 2 % 9 2。尾巴部分1到2中与3互质的数乘积为1*22。所以F(2,3,2) F(0,3,2) * prod^0 * 2 mod 9 1 * 1 * 2 2。回溯计算F(7,3,2)cyc floor(7/9)0tail 7 % 9 7。尾巴部分1到7中与3互质的数乘积剔除3,6乘积为1*2*4*5*7280 mod 91因为280/931*91。所以F(7,3,2) F(2,3,2) * prod^0 * 1 mod 9 2 * 1 * 1 2。回溯计算F(22,3,2)cyc floor(22/9)2tail 22 % 9 4。尾巴部分1到4中与3互质的数乘积剔除3乘积为1*2*48 mod 9。所以F(22,3,2) F(7,3,2) * prod^2 * 8 mod 9 2 * (8^2 mod 9) * 8 mod 9。8^264 mod 91。所以F(22,3,2) 2 * 1 * 8 mod 9 16 mod 9 7。因此u_22 7。类似地可以计算出u_8 F(8,3,2)u_14 F(14,3,2)。计算过程略假设我们算出u_8 4 mod 9,u_14 2 mod 9这里仅为示例实际需计算。计算逆元inv(u_8) mod 9即inv(4) mod 9因为4*728 mod 91所以逆元为7。inv(u_14) mod 9即inv(2) mod 9因为2*510 mod 91所以逆元为5。U u_22 * inv(u_8) * inv(u_14) mod 9 7 * 7 * 5 mod 9 245 mod 9 2。因为245/927*92组合结果C(22,8) mod 9 ≡ 3^2 * U mod 9 9 * 2 mod 9 18 mod 9 0。可以看到因为指数E2等于k2导致q^E项模q^k为0最终结果为0。这符合组合数C(22,8)是一个很大的数且显然能被9整除的直观理解因为分子中3的因子数远多于分母。4. 算法实现递归、预处理与细节打磨理解了原理实现起来就有了清晰的路线图。整个算法分为几个层次清晰的函数。4.1 快速幂与扩展欧几里得求逆元这是基础工具。虽然我们用到了求逆元但仅限于与q互质的数u对q^k取模的情况这保证了逆元一定存在。我们可以用扩展欧几里得算法求解a*x mod*y gcd(a, mod) 1中的x即为a模mod的逆元。// 快速幂 (a^b % mod) 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); // 解 ax mod*y 1 return (x % mod mod) % mod; }4.2 计算 n! 中质因子 q 的个数 e这个函数是纯数学计算没有模运算。long long get_factor_cnt(long long n, long long q) { long long cnt 0; while (n) { cnt n / q; n / q; } return cnt; }4.3 核心递归函数计算 F(n, q, k) u_n mod q^k这是算法的心脏。我们需要预处理循环节乘积prod。// 计算 n! 中剔除所有因子 q 后的结果对 q^k 取模 long long factorial_mod(long long n, long long q, long long qk) { if (n 0) return 1; long long res 1; // 1. 处理完整周期 // prod 需要预处理1*2*...*qk 中所有与 q 互质的数的乘积 % qk long long cyc n / qk; res res * pow_mod(prod, cyc, qk) % qk; // 2. 处理不完整的尾巴周期 long long tail n % qk; for (long long i 1; i tail; i) { if (i % q ! 0) { // 只乘与 q 互质的数 res res * i % qk; } } // 3. 递归处理 q 的倍数提取因子后形成的阶乘 res res * factorial_mod(n / q, q, qk) % qk; return res; }注意事项prod的计算需要小心。qk可能比较大比如5^53125直接连乘1到qk可能会很慢。一种优化方法是利用模运算的周期性只计算一个周期的乘积并缓存。另外在递归调用factorial_mod(n/q, q, qk)时参数n/q是整数除法确保了递归深度为O(log_q n)。4.4 计算 C(n, m) mod q^k这个函数组合前面的所有部分。long long C_mod_qk(long long n, long long m, long long q, long long k) { if (m 0 || m n) return 0; long long qk 1; for (int i 0; i k; i) qk * q; // 计算 q^k注意可能溢出可用快速幂 long long e_n get_factor_cnt(n, q); long long e_m get_factor_cnt(m, q); long long e_nm get_factor_cnt(n - m, q); long long e e_n - e_m - e_nm; // 如果 q 的指数已经超过 k直接返回 0 if (e k) return 0; long long u_n factorial_mod(n, q, qk); long long u_m factorial_mod(m, q, qk); long long u_nm factorial_mod(n - m, q, qk); long long U u_n * inv(u_m, qk) % qk * inv(u_nm, qk) % qk; long long res pow_mod(q, e, qk) * U % qk; return res; }4.5 利用中国剩余定理合成最终答案最后一步我们对p做质因数分解p p1^k1 * p2^k2 * ...对每个质因子幂pi^ki调用C_mod_qk得到余数a_i和模数m_i pi^ki。然后使用中国剩余定理求解同余方程组x ≡ a_i (mod m_i)。中国剩余定理的求解可以通过迭代方式实现 假设当前已经合并了前i-1个方程得到解x ≡ A (mod M)其中M m1*m2*...*m_{i-1}。 现在要加入第i个方程x ≡ a_i (mod m_i)。 我们需要找t使得A t*M ≡ a_i (mod m_i)。 即t*M ≡ a_i - A (mod m_i)。 这相当于求解M在模m_i下的逆元inv_M然后t (a_i - A) * inv_M mod m_i。 新的解为x A t * M新的模数为M M * m_i。long long CRT(const vectorlong long a, const vectorlong long m) { // a 是余数数组m 是模数数组要求 m 两两互质 long long x 0, M 1; for (int i 0; i a.size(); i) { // 合并方程x ≡ a[i] (mod m[i]) long long t (a[i] - x) % m[i]; if (t 0) t m[i]; long long inv_M inv(M % m[i], m[i]); // 求 M 模 m[i] 的逆元 if (inv_M -1) return -1; // 理论上不会发生因为 m 两两互质 t t * inv_M % m[i]; x x t * M; M M * m[i]; } return x % M; }4.6 主函数逻辑最终扩展卢卡斯的主函数流程如下输入n, m, p。对p进行质因数分解得到质因子列表q[]和指数列表k[]。对于每个(q[i], k[i])计算res_i C_mod_qk(n, m, q[i], k[i])以及模数mod_i q[i]^k[i]。将所有(res_i, mod_i)作为同余方程组调用CRT函数求解。输出结果。5. 边界条件、优化与常见踩坑点实现看起来清晰但魔鬼藏在细节里。下面是我在多次实现中总结出的关键点和易错点。5.1 数值溢出无处不在的“刺客”这是实现扩展卢卡斯时最大的挑战。n和m可以非常大10^18p也可能很大。中间乘法溢出在factorial_mod函数中res res * i % qk即使res和i都小于qk它们的乘积也可能超过long long的范围约9e18。如果qk接近1e9那么res*i就可能溢出。解决方案使用快速乘龟速乘或__int128。在无法使用__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; }将res * i % qk替换为mul_mod(res, i, qk)。注意这会带来一定的常数开销。计算q^k时溢出在C_mod_qk中计算qk pow(q, k)时q^k本身可能就超出了long long范围。但通常题目中p是在long long范围内的所以q^k作为p的因子也不会超限。不过在循环for (int i0; ik; i) qk * q;中如果q很大最后一次乘法可能导致溢出。更安全的做法是用一个循环并在每次乘法后判断是否超过p或LLONG_MAX/q。中国剩余定理合并时溢出x x t * M和M M * m_i这两步乘法极易溢出。M是所有已合并模数的乘积最终会等于p所以M和m_i都在long long范围内但它们的乘积在计算过程中可能溢出。同样需要使用快速乘。x (x mul_mod(t, M, M_new)) % M_new; // M_new M * m_i M M_new;5.2 递归深度与性能factorial_mod是递归函数递归深度约为log_q n。对于n1e18,q2深度约为60可以接受。但递归调用本身有开销。在极端情况下可以考虑写成非递归形式但递归的清晰性更好。主要的性能瓶颈在于factorial_mod函数中尾巴部分的循环for (i1 to tail)如果qk很大比如10^6这个循环会非常慢。但通常q^k作为p的因子p本身不会太大比如1e9以内所以qk也不会太大。如果遇到p很大且含有大质数幂的情况这个算法本身的时间复杂度就会变高这是理论上的限制。5.3 预处理“循环节乘积” prod 的优化在factorial_mod中我们假设prod是已知的。我们需要为每一组(q, k)预处理这个值。计算prod需要遍历1到q^k复杂度O(q^k)这在q^k较大时不可接受。优化技巧注意到prod是1到q^k中所有与q互质的数的乘积模q^k。这可以通过一个周期性的性质来快速计算。实际上prod ≡ -1 mod q^k当q为奇质数或q2且k3时这是一个广义威尔逊定理的推论。但对于q2, k3有prod ≡ 1 mod 2^k。我们可以利用这个性质直接得到prod避免遍历。若q为奇质数则prod q^k - 1。若q2且k1prod1。若q2且k2prod1模4下与2互质的数有1,3乘积为3≡-1 mod 4这里需要验证1*33模4余3即-1。所以对于q2, k2prod3。若q2且k3prod1。 严格来说prod是1到q^k中所有与q互质的数的乘积模q^k。对于奇质数q这个乘积确实同余于-1。对于q2需要单独处理。最稳妥的方法还是在初始化时用循环计算一次并缓存因为q^k通常不会太大。如果q^k真的很大就需要用到这个数论结论来优化。5.4 指数 e k 时的处理在C_mod_qk中如果e k我们直接返回0。这是正确的因为q^e项模q^k为0。但有一点需要注意即使e k我们仍然需要计算U吗理论上不需要了因为乘以0后结果总是0。但为了代码清晰和逻辑完整先判断e并提前返回0是更高效的做法。5.5 对 p 进行质因数分解这是算法的第一步。我们需要得到p ∏ pi^ki且pi是质数。可以使用试除法从2到sqrt(p)进行遍历。由于p可能很大1e9级别这个分解是O(sqrt(p))在p很大时可能成为瓶颈。但在算法竞赛中通常p不会太大或者其质因子很少。如果p非常大如10^18则需要更高效的分解算法如 Pollard-Rho这超出了扩展卢卡斯的范畴通常题目会避免这种情况。6. 总结与扩展思考扩展卢卡斯定理是一个将复杂问题分解、递归求解的经典案例。它完美地结合了数论中的几个核心工具质因数分解、阶乘质因子计数、递归、模逆元、中国剩余定理。掌握它不仅是为了解决一道特定的算法题更是对模运算和组合数计算有了更深的理解。回顾整个流程它的核心思想始终是“解决不可逆问题”通过提取公共质因子将原本不可逆的分母转化为可逆的部分从而在质数幂模数下完成计算最后用中国剩余定理拼回原模数下的答案。在实际编码中最需要警惕的是数值溢出。几乎每一个乘法操作都要考虑是否可能超出数据类型的范围养成使用快速乘的习惯在数论算法中至关重要。其次对于边界情况如m0,mn,ek要处理得当。虽然扩展卢卡斯定理在理论上能解决任意模数p下的组合数取模问题但其时间复杂度并非总是最优。当p的质因子很多或者某个q^k很大时算法的效率会下降。在一些特殊情况下可能有更优的专用算法。但对于通用场景扩展卢卡斯定理无疑是工具箱里一件强大而可靠的武器。最后理解算法的最好方式就是实现它。我建议你在理解上述所有步骤后关闭这篇文章尝试自己从头实现一遍。过程中你一定会遇到这里提到或没提到的各种问题而解决这些问题的过程才是真正将知识内化的关键。当你成功通过一道需要扩展卢卡斯的题目时那种成就感会让你觉得之前所有的思考和调试都是值得的。