卢卡斯定理:大数组合取模的高效算法原理与实现 1. 项目概述为什么我们需要卢卡斯定理在算法竞赛和组合数学的实际问题中我们经常需要计算组合数 C(n, m) 对某个质数 p 取模的结果。比如在计算一个复杂概率模型的方案数或者解决某些计数类动态规划问题时组合数取模是绕不开的一步。最直接的想法是使用公式 C(n, m) n! / (m! * (n-m)!)然后分别计算阶乘再求逆元。这个方法在 n 和 m 比较小的时候比如几千以内是可行的。但是当 n 和 m 的规模达到 10^18 级别而模数 p 是一个相对较小比如 10^5 以内的质数时问题就来了。你不可能去计算一个天文数字的阶乘时间和空间都不允许。这时候卢卡斯定理Lucas‘ Theorem就闪亮登场了。它提供了一种将大规模组合数取模问题分解为若干个规模仅在模数 p 以内的小问题的方法从而使得计算成为可能。简单来说它让“用大炮打蚊子”变成了“用若干个小弹弓精准射击”是处理大数组合模运算的利器。接下来我会带你彻底搞懂它的原理、证明并给出清晰可靠的代码模板和避坑指南。2. 卢卡斯定理的核心原理与证明2.1 定理的表述卢卡斯定理的内容非常简洁对于质数 p 和任意非负整数 n, m有 C(n, m) ≡ C(n mod p, m mod p) * C(n/p, m/p) (mod p)其中C(n/p, m/p) 这一项可以继续递归地用卢卡斯定理进行计算直到 m/p 为 0 为止。一个直观的例子计算 C(13, 5) mod 3。这里 p3, n13, m5。首先计算 n mod p 13 mod 3 1 m mod p 5 mod 3 2。所以有 C(1, 2)。在组合数学中当 m n 时C(n, m) 0。所以 C(1, 2) 0。因此根据卢卡斯定理C(13, 5) mod 3 ≡ 0 * C(13/3, 5/3) 0 * C(4, 1) 0。我们可以验证一下C(13, 5) 1287 1287 mod 3 0。结果正确。这个例子展示了定理的基本形式但要想真正掌握并放心使用我们必须理解它为什么成立。2.2 关键引理组合数与二项式系数证明卢卡斯定理需要一个关键的引理对于质数 p以及任意满足 1 ≤ k ≤ p-1 的整数 k组合数 C(p, k) 能被 p 整除即 C(p, k) ≡ 0 (mod p)。证明C(p, k) p! / (k! * (p-k)!)。分子 p! 包含因子 p而分母 k! 和 (p-k)! 都不包含因子 p因为 k 和 p-k 都小于 p。因此整个分数化简后因子 p 依然保留所以 C(p, k) 是 p 的整数倍。这个引理意味着在模 p 意义下(1 x)^p 的展开式有非常简单的形式 (1 x)^p ≡ 1 x^p (mod p) 因为中间的各项 C(p, k) * x^k 在模 p 下都等于 0。2.3 卢卡斯定理的完整证明生成函数法这是最经典和优美的证明方法。我们考虑生成函数 (1 x)^n 在模 p 意义下的两种展开方式。证明过程第一种展开直接二项式定理 (1 x)^n Σ_{m0}^{n} C(n, m) * x^m。我们关心的是系数 C(n, m) 模 p 的值。第二种展开利用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 其中每个 n_i 和 m_i 都是 0 到 p-1 之间的整数。现在我们观察 (1 x)^n (1 x)^n (1 x)^{n_0} * [(1 x)^p]^{n_1} * [(1 x)^{p^2}]^{n_2} * ... * [(1 x)^{p^k}]^{n_k}根据关键引理(1 x)^p ≡ 1 x^p (mod p)。那么对于更高的幂次有 (1 x)^{p^t} [(1 x)^p]^{p^{t-1}} ≡ (1 x^p)^{p^{t-1}} ≡ 1 x^{p^t} (mod p) 这里用到了模运算的性质和归纳法核心是每次升幂内部的 x 项指数都会乘以 p。因此在模 p 意义下 (1 x)^n ≡ (1 x)^{n_0} * (1 x^p)^{n_1} * (1 x^{p^2})^{n_2} * ... * (1 x^{p^k})^{n_k} (mod p)比较系数 我们现在想得到 x^m 项的系数。注意 m 的 p 进制表示为 m m_0 m_1 * p m_2 * p^2 ... m_k * p^k。 在右边的乘积中因子 (1 x)^{n_0} 贡献 x^{m_0} 的系数即 C(n_0, m_0)。因子 (1 x^p)^{n_1} 贡献 (x^p)^{m_1} x^{m_1*p} 的系数即 C(n_1, m_1)。因子 (1 x^{p^2})^{n_2} 贡献 (x^{p^2})^{m_2} x^{m_2*p^2} 的系数即 C(n_2, m_2)。以此类推。要最终得到 x^m 项我们必须从每个因子中分别取出恰好能构成 m 的对应 p 进制位的项。因为不同因子贡献的 x 的指数是 p 的不同幂次即 p^0, p^1, p^2, ...它们彼此独立互不干扰。所以x^m 项的系数就是所有这些独立组合数系数的乘积 C(n_0, m_0) * C(n_1, m_1) * ... * C(n_k, m_k) (mod p)得出结论 由于第一种展开中 x^m 的系数是 C(n, m)第二种展开得到的是 Π C(n_i, m_i)因此 C(n, m) ≡ Π_{i0}^{k} C(n_i, m_i) (mod p)这个连乘式子恰恰就是卢卡斯定理递归形式的结果C(n, m) ≡ C(n mod p, m mod p) * C(n/p, m/p) (mod p)。递归下去就是将 n 和 m 不断除以 p取每一位的数字进行组合数计算。注意这个证明中有一个隐含条件即如果某一位上 m_i n_i那么 C(n_i, m_i) 0从而导致整个乘积为 0。这对应了组合数 C(n, m) 在模 p 下为 0 的情况。在实际计算中我们需要在递归基或预处理阶乘时处理这种情况。3. 算法实现与代码模板详解理解了原理实现就清晰了。卢卡斯定理的实现通常分为两部分预处理阶乘和逆元用于快速计算小规模组合数以及递归计算卢卡斯定理本身。3.1 预处理阶乘与逆元由于卢卡斯定理最终会将问题分解为计算多个 C(n_i, m_i)其中 n_i, m_i p。我们可以预处理出 0 到 p-1 的所有阶乘 fact[i] 模 p 的值以及对应的逆元 invfact[i]。为什么要预处理逆元因为计算组合数公式 C(a, b) a! / (b! * (a-b)!)在模 p 下除法需要转换为乘以逆元。即 C(a, b) ≡ fact[a] * invfact[b] * invfact[a-b] (mod p)。预处理后每次小组合数计算就是 O(1) 的。预处理步骤时间复杂度 O(p)初始化 fact[0] 1。循环计算 fact[i] fact[i-1] * i % p。计算 invfact[p-1] pow(fact[p-1], p-2, p)。这里利用费马小定理a^{p-1} ≡ 1 (mod p)所以 a 的逆元是 a^{p-2}。逆推计算 invfact[i] invfact[i1] * (i1) % p。因为 invfact[i] ≡ 1/(i!) ≡ 1/((i1)! / (i1)) ≡ invfact[i1] * (i1) (mod p)。3.2 核心递归函数 Lucas(n, m)这是卢卡斯定理的递归实现。算法流程递归基如果 m 0根据定义C(n, 0) 1直接返回 1。分解与递归否则计算 n_i n % p, m_i m % p。如果 m_i n_i根据组合数定义C(n_i, m_i) 0那么整个结果就是 0直接返回 0。否则计算小规模组合数C(n_i, m_i)使用预处理的阶乘和逆元fact[n_i] * invfact[m_i] % p * invfact[n_i - m_i] % p。递归调用计算大规模部分C(n/p, m/p)即Lucas(n/p, m/p)。合并结果将两部分相乘并对 p 取模返回结果。时间复杂度分析递归的深度是 n 的 p 进制位数即 O(log_p n)。每次递归中的小组合数计算是 O(1)。因此总时间复杂度为 O(log_p n p)其中 O(p) 是预处理的时间。当 p 相对较小时这个算法非常高效。3.3 完整C代码模板下面是一个包含详细注释的模板适用于大多数算法竞赛场景。#include iostream using namespace std; typedef long long ll; // 快速幂用于计算逆元a^b mod p ll qpow(ll a, ll b, ll p) { ll res 1; while (b) { if (b 1) res res * a % p; a a * a % p; b 1; } return res; } // 预处理阶乘和阶乘的逆元 ll fact[100005], invfact[100005]; // 数组大小根据模数p的最大值设定 void init(ll p) { fact[0] 1; for (int i 1; i p; i) { fact[i] fact[i-1] * i % p; } // 费马小定理求 (p-1)! 的逆元 invfact[p-1] qpow(fact[p-1], p-2, p); // 逆推得到所有阶乘逆元 for (int i p-2; i 0; --i) { invfact[i] invfact[i1] * (i1) % p; } } // 计算小组合数 C(a, b) % p要求 a, b p ll C_small(ll a, ll b, ll p) { if (b a) return 0; // 组合数定义b不能大于a // C(a, b) a! / (b! * (a-b)!) return fact[a] * invfact[b] % p * invfact[a - b] % p; } // 卢卡斯定理递归函数 ll Lucas(ll n, ll m, ll p) { if (m 0) return 1; // 递归基 // C(n, m) % p C(n%p, m%p) * Lucas(n/p, m/p, p) % p return C_small(n % p, m % p, p) * Lucas(n / p, m / p, p) % p; } int main() { int T; // 询问次数 ll n, m, p; cin T; while (T--) { cin n m p; // 每次输入的模数p可能不同 init(p); // 每次根据新的p进行预处理 cout Lucas(n, m, p) endl; } return 0; }4. 模板使用中的关键细节与避坑指南模板看起来简单但实际使用时细节决定成败。下面是我在多次比赛中总结出的经验和常见问题。4.1 模数p必须是质数这是卢卡斯定理成立的前提。如果p不是质数关键引理(1x)^p ≡ 1x^p (mod p)不一定成立整个证明的基石就塌了。在使用前务必确认p是质数。在竞赛中题目通常会明确给出“素数p”。如果p可能不是质数则需要使用扩展卢卡斯定理ExLucas那是另一个更复杂的算法。4.2 预处理数组的大小模板中fact和invfact数组的大小设置为100005这是假设模数p最大为10^5量级。你必须根据题目中p的最大可能值来调整这个数组大小。通常卢卡斯定理适用于“大n, 小p”的场景p一般在10^5到10^6以内。如果p更大预处理O(p)的时间可能会超时此时卢卡斯定理可能不再是最优选择。4.3 多组询问的初始化注意模板的main函数中对于每组数据(n, m, p)都调用了init(p)。这是因为不同的询问可能对应不同的模数p。如果所有询问的模数p相同那么init(p)只需要调用一次可以大幅提升效率。务必根据题目描述判断。4.4 关于long long的使用n和m可能非常大10^18所以必须使用long long或__int128。在计算n % p和n / p时long long是安全的。在乘法运算fact[a] * invfact[b] % p中虽然fact[a]和invfact[b]都小于p但两者的乘积可能超过int范围因此在乘法前最好先转为long long或直接使用long long类型的数组。我的模板中全部使用了ll(long long)。4.5 递归深度与栈溢出递归深度是O(log_p n)。对于极大的n如10^18和极小的p如2深度约为60这在任何评测系统中都是安全的不会导致栈溢出。但如果你非常担心也可以写成非递归的循环形式ll Lucas_iterative(ll n, ll m, ll p) { ll res 1; while (n m) { res res * C_small(n % p, m % p, p) % p; if (res 0) return 0; // 提前终止优化 n / p; m / p; } return res; }4.6 特判与边界条件m n 的情况在组合数中如果 m n则 C(n, m) 0。这个检查发生在C_small函数中if (b a) return 0;。在递归过程中如果某一位出现m_i n_i会立即返回0并且由于乘法性质最终结果也是0。这是正确的。m 0 的情况在Lucas函数中作为递归基处理返回1。p 非常小的情况当 p 很小比如2、3、5时预处理数组也很小算法会运行得飞快。但要注意此时组合数取模后为0的概率可能会变高这在一些计数问题中可能有特殊含义。5. 典型应用场景与问题剖析卢卡斯定理不是孤立的算法它总是作为工具嵌入到更大的问题中。理解它的应用场景才能更好地识别何时该用它。5.1 场景一大组合数取模的直接计算这是最直白的应用。题目直接要求计算 C(n, m) % p其中 1 ≤ n, m ≤ 10^18 1 ≤ p ≤ 10^5 且 p 为质数。例题特征输入格式通常就是 T组数然后每组给出 n, m, p。解法直接套用上述模板即可。5.2 场景二复杂计数问题中的子问题很多动态规划或排列组合问题其状态转移方程或最终答案表达式中包含组合数且 n, m 可能非常大。例题特征问题最终转化为求 Σ C(a_i, b_i) 或类似形式且 a_i, b_i 范围很大。解法首先推导出问题的组合数学表达式。如果发现表达式中的组合数满足“大n, 小质数p”的条件就可以在计算每个组合数时使用卢卡斯定理。通常需要结合其他算法如数位DP、容斥原理等。示例模型求在 [L, R] 区间内有多少个数的二进制表示中恰好有 k 个 1。这可以转化为数位DP。在DP过程中当我们确定高位后低位可以自由选择需要计算在剩余位数中选若干个位置填1的方案数即组合数。由于 n剩余位数可能很大但模数p通常是结果对某个质数取模不大就可以用卢卡斯定理快速计算这个组合数模 p 的值。5.3 场景三与费马小定理/欧拉定理的对比选择初学者容易混淆何时用费马小定理求逆元何时用卢卡斯定理。费马小定理/扩展欧几里得求逆元适用于计算C(n, m) % p其中n 和 m 可以很大但 p 更大通常需要 n p。因为你需要计算 n! % p如果 n p那么 n! 中包含了因子 p模 p 后为 0导致逆元不存在因为分母有 p 的因子不可逆。所以这种方法要求 n, m p。卢卡斯定理正是为了解决n, m 远大于 p的情况。它将大问题化归到 p 以内的小问题。选择策略读题看数据范围。如果 n, m ≤ 10^6 p ~ 10^97用预处理阶乘逆元费马小定理。如果 n, m ≤ 10^18 p ≤ 10^6用卢卡斯定理。如果 n, m ≤ 10^18 p 不是质数用扩展卢卡斯定理ExLucas。如果 n, m ≤ 5000甚至可以用杨辉三角递推。6. 性能优化与扩展讨论6.1 预处理优化针对固定模数如果模数 p 固定且有多组询问预处理只需做一次。可以将init(p)放在所有询问之前。更进一步如果 p 是常用质数如 1e97, 998244353甚至可以预先写好它们的阶乘和逆元数组当然对于1e97n通常不会超过它直接用法一更常见。6.2 记忆化递归在递归函数Lucas(n, m, p)中参数是(n, m)。对于不同的询问可能会重复计算相同的(n, m)对。如果询问次数极多可以考虑用mappairll, ll, ll存储已经计算过的结果避免重复递归。但通常来说递归深度很浅记忆化带来的提升有限反而增加了 map 的开销需要根据实际情况权衡。6.3 扩展卢卡斯定理ExLucas简介当模数 p 不是质数时标准卢卡斯定理失效。此时需要使用扩展卢卡斯定理。其核心思想是将模数 p 分解质因数p p1^k1 * p2^k2 * ... * pt^kt。然后分别计算 C(n, m) mod pi^ki最后用中国剩余定理CRT合并结果。计算 C(n, m) mod p^k 是 ExLucas 的难点。它需要处理阶乘中 p 因子的剔除因为分母可能包含 p在模 p^k 下不一定有逆元。具体步骤是将 n!, m!, (n-m)! 中的因子 p 全部提取出来单独计算。对于剔除 p 因子后的部分由于它与 p^k 互质可以用扩展欧几里得求逆元。最后将两部分结合。ExLucas 的实现比 Lucas 复杂得多时间复杂度也高。除非题目明确要求否则在竞赛中较少遇到。但了解其存在性和解决思路是必要的。6.4 调试技巧与测试数据自己编写卢卡斯定理代码时如何验证正确性小数据暴力验证写一个暴力计算组合数用高精度或直接算然后取模的程序与你的卢卡斯算法在 n, m 较小比如100时进行对拍。利用已知性质C(n, m) C(n, n-m)。用你的程序验证这个等式。帕斯卡恒等式C(n, m) C(n-1, m-1) C(n-1, m)。选择中等大小的 n, m 进行验证。边界测试m 0, m n。n 很大m 很小或很大。p 2最小的质数。测试m_i n_i导致结果为 0 的情况。7. 从理论到实战一道例题的完整分析让我们通过一道虚构但典型的题目来串联所有知识点。题目给定质数 p (p ≤ 10007)和 T (T ≤ 100) 组询问每组询问给出 n, m (0 ≤ m ≤ n ≤ 10^18)求 C(n, m) mod p。分析数据范围分析n, m 高达 10^18远超 p (≤10007)。这是典型的“大n, 小质数p”场景明确指向卢卡斯定理。算法选择直接使用标准卢卡斯定理。预处理阶乘和逆元数组大小为 p最大10007完全可行。实现细节由于 p 在每组询问中可能不同我们需要对每组询问重新调用init(p)。虽然 T100p10007预处理 O(T*p) 的复杂度约为 10^6可以接受。如果题目保证所有询问 p 相同则只需初始化一次。编写代码直接套用第3部分的模板。测试测试1p10007, n123456789, m98765432。用程序计算。测试2p7, n100, m50。可以手算验证将100和50转化为7进制。100 的 7 进制202 (因为 249 07 2 100)50 的 7 进制101 (因为 149 07 1 50)根据卢卡斯定理C(100,50) mod 7 ≡ C(2,1) * C(0,0) * C(2,1) mod 7。C(2,1)2, C(0,0)1, C(2,1)2。乘积为 4。所以结果应为 4。用程序验证。通过这样完整的分析、实现和验证流程你就能确保卢卡斯定理的代码在实战中万无一失。记住在竞赛中看到巨大的 n, m 和较小的质数 p你的第一反应就应该是卢卡斯定理。