Lucas定理精讲：大组合数取模的算法原理与实战

1. 从一个竞赛题说起系数、组合数与模运算的碰撞前几天在牛客的算法竞赛题目里又看到了一个老朋友——“系数”问题。这类问题通常长这样给你一个二项式展开式比如(x y)^n问你展开后x^k * y^(n-k)这一项的系数是多少。当然如果只是求系数那就是简单的组合数C(n, k)。但竞赛题从来不会这么简单它往往会加上一个模数p并且这个p可能不是质数或者n和k的规模巨大比如n可以大到10^18。这时直接计算组合数再取模的路就走不通了因为中间结果会溢出阶乘和逆元的计算也会因为模数非质数而失效。这就是Lucas定理的经典应用场景。它像一把精巧的钥匙专门用来解决“大组合数对大质数取模”的计算难题。今天我们就以这道题为引子彻底拆解Lucas定理的原理、实现细节以及那些在竞赛和实际编码中容易踩的坑。2. 问题本质为什么大组合数取模这么棘手在深入Lucas定理之前我们必须先理解它要解决的核心矛盾。计算组合数C(n, m) n! / (m! * (n-m)!)在计算机中通常有两种思路一是利用阶乘公式直接计算二是利用递推公式如杨辉三角。但当n和m很大比如超过10^7时这两种方法都会遇到天花板。2.1 直接计算的瓶颈溢出与精度即使不考虑取模用double或long double存储阶乘精度也会迅速丢失。用整数类型则必然面临溢出问题。例如100!的结果是一个158位的天文数字远超任何基本数据类型的表示范围。2.2 取模运算带来的新问题除法与逆元当我们引入取模运算希望计算C(n, m) % p时情况变得更复杂。取模运算对加法和乘法是友好的(ab)%p (a%p b%p)%p(a*b)%p (a%p * b%p)%p。但是它对除法不友好。我们不能直接计算(n! / (m! * (n-m)!)) % p因为除法在模运算下没有直接的定义。解决方案是引入模逆元。如果p是一个质数并且a不是p的倍数那么根据费马小定理a^(p-1) ≡ 1 (mod p)因此a的逆元inv(a) a^(p-2) % p。这样我们就可以把除法转化为乘法(a / b) % p (a * inv(b)) % p。因此对于质数模数p计算C(n, m) % p的标准流程是预处理出1!到n!对p取模的结果fact[i]。预处理出1!到n!的逆元inv_fact[i]通常通过递推inv_fact[i] inv_fact[i1] * (i1) % p来高效计算。组合数C(n, m) fact[n] * inv_fact[m] % p * inv_fact[n-m] % p。这个算法的时间复杂度是O(n)预处理O(1)查询非常高效。但它有一个致命前提n必须小于p。2.3 当 n 大于等于 p 时卢卡斯定理的舞台为什么n p是前提考虑一个极端例子计算C(p, 1) % p。根据组合数公式它应该等于p! / (1! * (p-1)!) p。在模p下结果应该是0。但是如果我们用上述的逆元法来计算fact[p] p! % p 0因为p!包含因子p。inv_fact[1]和inv_fact[p-1]都存在。最终计算0 * inv_fact[1] * inv_fact[p-1] % p会得到0。看起来结果是对的但这里存在一个严重的逻辑漏洞当fact[p] 0时inv_fact[p]是不存在的因为0没有乘法逆元。实际上我们预处理阶乘和逆元时通常只处理到p-1因为p! % p 0从p开始的阶乘模p都是0整个逆元链在这里就断了。所以当n p时fact[n]很可能为0导致基于逆元的算法完全失效。这就是Lucas定理要解决的问题。它不直接计算巨大的C(n, m)而是巧妙地将n和m用p进制表示将大问题分解为若干个规模小于p的子问题而这些子问题正是我们擅长的、可以用逆元法高效解决的C(n_i, m_i) % p。3. Lucas定理的数学原理与证明思路Lucas定理的描述非常简洁优美 对于非负整数n,m和质数p将n和m表示为p进制数n n_k * p^k ... n_1 * p n_0m m_k * p^k ... 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。3.1 一个直观的例子假设p7,n100,m30。 首先将n和m转化为7进制100 ÷ 7 14 ... 2-n0 214 ÷ 7 2 ... 0-n1 02 ÷ 7 0 ... 2-n2 2所以100的7进制是(2,0,2)_7即2*7^2 0*7^1 2*7^0。同理30的7进制是(4,2)_7即30 4*7 2补齐高位为(0,4,2)_7。 根据Lucas定理C(100, 30) % 7 ≡ C(2,0) * C(0,4) * C(2,2) % 7因为C(0,4)中40所以该项为0。 因此C(100, 30) % 7 ≡ 1 * 0 * 1 ≡ 0 (mod 7)。 这个结果告诉我们C(100,30)能被7整除。3.2 定理的证明思路生成函数法理解证明能帮助我们更深刻地把握定理的适用边界。一个经典的证明利用的是二项式系数与生成函数的关系。核心观察是在模p下二项式定理(1x)^n有特殊的性质。根据费马小定理对于质数p有(1x)^p ≡ 1 x^p (mod p)。这是因为展开式(1x)^p Σ C(p, i) x^i中除了i0和ip的项系数为1中间项C(p, i)都包含因子p在模p下均为0。现在将n写成p进制n n0 n1*p n2*p^2 ...。那么(1x)^n (1x)^{n0} * [(1x)^p]^{n1} * [(1x)^{p^2}]^{n2} * ...利用上面的性质在模p下(1x)^p ≡ 1x^p(1x)^{p^2} ((1x)^p)^p ≡ (1x^p)^p ≡ 1x^{p^2}以此类推。 所以(1x)^n ≡ (1x)^{n0} * (1x^p)^{n1} * (1x^{p^2})^{n2} * ... (mod p)我们现在考虑等式两边x^m项的系数。左边(1x)^n中x^m的系数就是C(n, m)。 右边是若干个形如(1x^{p^i})^{n_i}的乘积。要得到x^m我们需要从每个因式中分别取出x^{m_i * p^i}项并且满足m m0 m1*p m2*p^2 ...这恰好是m的p进制表示。而从(1x^{p^i})^{n_i}中取出x^{m_i * p^i}项的系数是C(n_i, m_i)。 因此比较两边x^m的系数就得到了C(n, m) ≡ Π C(n_i, m_i) (mod p)这个证明清晰地揭示了Lucas定理的本质它将一个关于(1x)^n的模p计算分解为多个相互独立的、规模更小的关于(1x^{p^i})^{n_i}的计算。注意证明中关键的一步(1x)^p ≡ 1x^p (mod p)依赖于p是质数。如果p不是质数这个同余式不成立因此Lucas定理只适用于模数 p 为质数的情况。对于合数模数需要用到其扩展形式如中国剩余定理CRT结合质因数分解这又是另一个复杂的话题了。4. Lucas定理的算法实现与代码细节理解了原理实现起来就清晰了。算法流程分为两步子问题计算实现一个函数C(n, m, p)用于计算当n, m p时的组合数C(n, m) % p。这里就用到了我们第二节提到的预处理阶乘和逆元的方法。递归分解实现Lucas(n, m, p)函数递归地将大问题C(n, m) % p分解为C(n%p, m%p, p) * Lucas(n/p, m/p, p) % p直到n或m为0。4.1 基础工具快速幂与逆元在实现C(n, m, p)之前我们需要两个基础工具快速幂用于计算a^b % p复杂度O(log b)。逆元计算根据费马小定理a在模质数p下的逆元为a^(p-2) % p可以用快速幂计算。// 快速幂模板 (a^b % p) long long qpow(long long a, long long b, long long p) { long long res 1; while (b) { if (b 1) res res * a % p; a a * a % p; b 1; } return res; } // 求逆元 (费马小定理要求 p 是质数且 a 不是 p 的倍数) long long inv(long long a, long long p) { return qpow(a, p - 2, p); }4.2 计算小范围组合数 C(n, m, p)这里假设n, m p。我们采用预处理阶乘数组fact和阶乘逆元数组inv_fact的方式实现O(1)查询。// 预处理阶乘和阶乘逆元范围是 [0, p-1] vectorlong long fact, inv_fact; void init_fact(int n, long long p) { fact.resize(n 1); inv_fact.resize(n 1); fact[0] 1; for (int i 1; i n; i) { fact[i] fact[i - 1] * i % p; } // 计算 n! 的逆元然后倒推得到所有阶乘的逆元 inv_fact[n] inv(fact[n], p); for (int i n - 1; i 0; --i) { inv_fact[i] inv_fact[i 1] * (i 1) % p; } } // 计算 C(n, m) % p, 调用前需确保 n, m p 且已调用 init_fact(p-1, p) long long comb(long long n, long long m, long long p) { if (m n) return 0; // 利用公式 C(n, m) n! / (m! * (n-m)!) return fact[n] * inv_fact[m] % p * inv_fact[n - m] % p; }4.3 实现Lucas定理递归函数这是核心函数它处理任意大的n和m。// Lucas定理主函数 long long lucas(long long n, long long m, long long p) { if (m 0) return 1; // C(n, 0) 1 // 递归分解C(n, m) % p C(n%p, m%p) * Lucas(n/p, m/p) % p return comb(n % p, m % p, p) * lucas(n / p, m / p, p) % p; }4.4 完整代码示例与调用将以上部分组合起来并处理边界情况。#include iostream #include vector using namespace std; vectorlong long fact, inv_fact; long long qpow(long long a, long long b, long long p) { /* 同上 */ } long long inv(long long a, long long p) { /* 同上 */ } void init_fact(int n, long long p) { /* 同上 */ } long long comb(long long n, long long m, long long p) { /* 同上 */ } long long lucas(long long n, long long m, long long p) { if (m 0) return 1; // 注意这里 comb 的参数 n%p, m%p 一定小于 p是安全的 return comb(n % p, m % p, p) * lucas(n / p, m / p, p) % p; } int main() { int T; // 测试用例数 cin T; while (T--) { long long n, m, p; cin n m p; // 初始化阶乘表范围只需到 p-1 init_fact(p - 1, p); // 计算 C(n, m) % p cout lucas(n, m, p) endl; // 注意如果有多组测试且 p 不同需要每次重新初始化 fact 和 inv_fact } return 0; }5. 实战中的边界条件与易错点分析代码写出来能跑通样例只是第一步在竞赛和工程中各种边界情况和细节才是真正的挑战。5.1 模数 p 不是质数怎么办这是最常被忽略的一点。Lucas定理要求 p 必须是质数。如果题目给的p是合数比如p10007它确实是质数但p1000就不是直接套用上述代码会得到错误结果。对于合数模数标准做法是对p进行质因数分解p p1^a1 * p2^a2 * ... * pk^ak。对每个质因子幂pi^ai计算C(n, m) mod pi^ai。这里需要用到扩展Lucas定理(ExLucas)其思想是将阶乘n!中所有pi的因子提出来单独计算剩余部分再用逆元处理。实现起来比普通Lucas复杂得多。最后利用中国剩余定理(CRT) 将k个同余方程的解合并得到C(n, m) mod p。在竞赛中如果p不是质数但n, m不大比如10^5级别有时可以直接用质因数分解组合数的方法或者用动态规划预处理避免使用Lucas定理。5.2 当 m_i n_i 时的处理根据Lucas定理在分解后的某一层如果m_i n_i则C(n_i, m_i) 0从而导致整个乘积为0。这在我们的comb函数中通过if (m n) return 0;已经处理了。这意味着只要m的p进制表示中有任何一位的数字大于n对应位的数字那么C(n, m)就能被p整除。这是一个非常有用的性质可以快速判断组合数的整除性。5.3 递归深度与栈溢出Lucas函数的递归深度等于n的p进制位数即O(log_p n)。对于n高达10^18p最小为2最大深度约为60这在绝大多数编程环境栈空间默认几MB下是安全的。但如果p很大接近n递归深度为1更没问题。所以通常不必担心栈溢出。当然写成非递归的迭代形式也是可以的但递归写法更直观。5.4 预处理的范围与多组测试这是一个极易出错的细节。我们的init_fact函数预处理了0到p-1的阶乘和逆元。这是因为在comb函数中参数n和m是n%p和m%p它们严格小于p。因此预处理到p-1就足够了预处理到n是巨大浪费n可能10^18而p通常10^5左右。在多组测试用例且每组模数p相同时可以只初始化一次阶乘表大幅提升效率。但如果每组测试的p不同必须在每组测试开始前根据新的p重新初始化fact和inv_fact数组。忘记重新初始化会导致使用错误的模数计算阶乘和逆元结果必然错误。一种安全的做法是将初始化过程直接放在lucas函数或每组计算的开头。5.5 时间复杂度分析预处理阶乘和逆元O(p)。这是算法的主要开销因为p通常需要能存入数组所以p一般在10^5到10^6量级。快速幂求逆元O(log p)但通过递推预处理阶乘逆元我们将单次求逆元优化到了O(1)。Lucas递归O(log_p n)次递归调用每次调用执行一次O(1)的comb查询。 因此总时间复杂度为O(p log_p n)对于单次查询O(p)的预处理是瓶颈。如果有Q次查询p相同则均摊复杂度为O(p Q * log_p n)非常高效。6. 回到原题牛客集训营B题“系数”解析虽然题目正文没有给出但根据标题“2021牛客寒假算法基础集训营6 B系数”和关键词“lucas定理”我们可以推断出题目的典型面貌。这类“系数”题通常是这样题目描述给定一个二项式(ax by)^n求展开后x^k * y^(n-k)项的系数对某个质数p取模的结果。其中n很大10^18k在0到n之间a,b是常数p是一个质数如1e97。解题思路根据二项式定理(axby)^n展开后x^k * y^{n-k}项是C(n, k) * (ax)^k * (by)^{n-k}。因此该项的系数为C(n, k) * a^k * b^{n-k}。我们需要计算[C(n, k) * a^k * b^{n-k}] % p。由于n巨大C(n, k)需要用 Lucas定理计算。a^k % p和b^{n-k} % p可以用快速幂计算。最终答案为lucas(n, k, p) * qpow(a, k, p) % p * qpow(b, n-k, p) % p。需要注意的陷阱常数 a 或 b 可能是负数在取模运算中需要先将a和b对p取模调整到[0, p-1)范围内避免负数取模带来问题。(a % p p) % p是一个安全的做法。模数 p 可能不是 1e97虽然1e97很常见但一定要用题目给的模数不能想当然。特判 k n 的情况虽然根据组合数定义此时系数为0但程序中也应处理直接返回0。代码实现框架#include bits/stdc.h using namespace std; typedef long long ll; // ... 省略 qpow, inv, init_fact, comb, lucas 的实现 ... int main() { int T; cin T; while (T--) { ll n, k, a, b, p; cin n k a b p; // 处理负数常数 a (a % p p) % p; b (b % p p) % p; // 初始化阶乘表 init_fact(p-1, p); // 计算答案 ll comb_val lucas(n, k, p); ll pow_a qpow(a, k, p); ll pow_b qpow(b, n - k, p); ll ans comb_val * pow_a % p * pow_b % p; cout ans endl; } return 0; }7. 不止于竞赛Lucas定理的应用场景与扩展Lucas定理的价值远不止于解决算法竞赛题。它在密码学、组合数学和某些特定领域的工程计算中也有应用。7.1 判断组合数的奇偶性一个经典结论是组合数C(n, m)为奇数当且仅当在二进制下m的每一位都不大于n的对应位。这其实就是p2时的Lucas定理特例。因为C(0,0)1,C(0,1)0,C(1,0)1,C(1,1)1。所以C(n, m) % 2等于所有二进制位对应组合数的乘积。只有每一位C(n_i, m_i)都不为0即m_i n_i最终结果才为1奇数。7.2 数论问题的分解有些数论问题涉及到大组合数模小质数。Lucas定理提供了一种将“大数”问题转化为“小数”问题的思路。例如证明某个组合数公式在模p下成立可能只需要验证所有n_i, m_i p的情况成立即可这大大简化了证明。7.3 扩展Lucas (ExLucas) 简介当模数p不是质数时我们需要扩展Lucas。其核心思想是将p质因数分解为∏ pi^ai。对于每个pi^ai计算C(n, m) mod pi^ai。计算时需要将n!,m!,(n-m)!中所有pi的因子提取出来单独计算这部分对模数的贡献可能为0剩余部分由于与pi互质可以求逆元。这个过程需要用到阶乘模质数幂的计算技巧例如n! mod p^k可以递归计算。最后用中国剩余定理合并所有结果。ExLucas的实现复杂度远高于普通Lucas在竞赛中通常作为“模板题”出现。理解其原理有助于深化对模运算和组合数计算的认识。7.4 一个重要的优化当 p 很小时如果模数p非常小比如p2, 3, 5那么n和m的p进制表示会非常长递归深度log_p n很大。但与此同时预处理阶乘表fact的规模p又非常小。此时虽然递归调用次数多但每次comb计算是O(1)的总复杂度O(p log_p n)仍然可以接受。极端情况下p2递归深度约60预处理只需计算fact[0]和fact[1]速度极快。8. 调试技巧与常见问题排查即使理解了所有原理实现时仍可能遇到各种问题。以下是一些实用的调试技巧。8.1 验证小数据用暴力计算动态规划或直接公式计算验证n, m 15且p较小如p7, 11时你的Lucas算法结果是否正确。这是检验算法逻辑最基本的方法。8.2 中间结果打印在递归函数中打印n, m, n%p, m%p, comb(n%p, m%p)等中间值与手算的p进制分解过程对比看分解和计算是否正确。8.3 检查模数 p 是否为质数写一个简单的质数判断函数试除法到sqrt(p)确保输入的p是质数。如果p不是质数你的结果肯定不对。这是使用Lucas定理的先决条件。8.4 多组数据初始化这是最隐蔽的bug之一。如果你的程序对第一组数据正确后续数据随机错误很可能是fact和inv_fact全局数组没有针对新的模数p重新初始化。确保每次p变化时都重新resize和计算这两个数组。8.5 数据范围与溢出注意long long的溢出。在qpow,comb函数中两个long long相乘可能溢出即使最后会取模。例如a * a % p如果a接近10^9a*a就会超过64位整数范围。安全的做法是使用慢速乘或直接使用__int128如果编译器支持。在竞赛中常见的处理方法是// 使用 __int128 临时存储 long long mul(long long a, long long b, long long p) { return (__int128)a * b % p; } // 或者在快速幂和组合数计算中使用 long long res (__int128)fact[n] * inv_fact[m] % p * inv_fact[n-m] % p;如果不支持__int128则需要实现一个O(log b)的慢速乘函数来替代直接乘法。8.6 特判 m0 或 mn在递归终点lucas函数中我们判断if (m 0) return 1;。同样在comb函数中C(n,0)C(n,n)1。这些特判能避免不必要的计算也符合数学定义。掌握Lucas定理就像是获得了一把处理大组合数模质数问题的瑞士军刀。它背后的数论思想进制分解、同余性质非常深刻。从理解为什么需要它到掌握其证明和实现再到注意各种边界条件和应用扩展这个过程本身就是一个很好的算法思维训练。下次再遇到“系数”和“大组合数取模”你应该能自信地拿起Lucas定理这把利器了。在实际编码时我习惯把qpow,inv,init_fact,comb,lucas这几个函数作为一个模板保存好并特别注意模数变化时的重新初始化问题这能避免很多不必要的麻烦。