算法概述
卢卡斯定理(Lucas Theorem)是数论中用于计算大组合数取模的经典定理。当需要计算 $C(n, m) \mod p$(其中 $p$ 为质数)且 $n, m$ 非常大时(例如 $n, m \leq 10^{18}$),直接计算阶乘会导致数值溢出,而卢卡斯定理通过将问题分解为 $p$ 进制下的子问题,巧妙地将大规模组合数取模转化为若干小规模组合数的乘积。
该定理在组合数学、密码学和算法竞赛中具有广泛应用,是解决大数组合数问题的核心工具之一。
数学原理
卢卡斯定理
设 $p$ 为质数,将 $n$ 和 $m$ 表示为 $p$ 进制数:
$$
n = n_k p^k + n_{k-1} p^{k-1} + \cdots + n_1 p + n_0
$$
$$
m = m_k p^k + m_{k-1} p^{k-1} + \cdots + m_1 p + m_0
$$
其中 $0 \leq n_i, m_i < p$。则:
$$
C(n, m) \equiv \prod_{i=0}^{k} C(n_i, m_i) \pmod{p}
$$
关键推论
若存在某一位 $m_i > n_i$,则 $C(n_i, m_i) = 0$,因此 $C(n, m) \equiv 0 \pmod{p}$。
小组合数计算
对于 $C(n_i, m_i)$ 其中 $n_i, m_i < p$,可通过预处理阶乘和逆元在 $O(1)$ 时间内求解:
$$
C(n, m) = \frac{n!}{m!(n-m)!} \equiv n! \cdot (m!)^{-1} \cdot ((n-m)!)^{-1} \pmod{p}
$$
其中逆元利用费马小定理计算:$a^{-1} \equiv a^{p-2} \pmod{p}$。
Java完整实现
import java.util.Arrays;
/**
* 卢卡斯定理(Lucas Theorem)实现大组合数取模
* 时间复杂度:单次查询 O(log_p(n) * p),预处理 O(p)
* 适用条件:模数p为质数
*/
public class LucasTheorem {
private final int p; // 质数模数
private final long[] fact; // 阶乘数组 fact[i] = i! % p
private final long[] invFact; // 阶乘逆元数组 invFact[i] = (i!)^{-1} % p
/**
* 构造函数:预处理阶乘和逆元
* @param p 质数模数,必须满足 p <= 2*10^6 以控制预处理空间
*/
public LucasTheorem(int p) {
if (!isPrime(p)) {
throw new IllegalArgumentException("Modulus must be prime, got: " + p);
}
this.p = p;
this.fact = new long[p];
this.invFact = new long[p];
// 预处理阶乘
fact[0] = 1;
for (int i = 1; i < p; i++) {
fact[i] = fact[i - 1] * i % p;
}
// 费马小定理求阶乘逆元:invFact[p-1] = (p-1)!^{p-2} mod p
invFact[p - 1] = modPow(fact[p - 1], p - 2, p);
// 递推求其他逆元:invFact[i-1] = invFact[i] * i % p
for (int i = p - 1; i >= 1; i--) {
invFact[i - 1] = invFact[i] * i % p;
}
}
/**
* 快速幂:计算 base^exp mod mod
*/
private long modPow(long base, long exp, long mod) {
long result = 1;
base %= mod;
while (exp > 0) {
if ((exp & 1) == 1) {
result = result * base % mod;
}
base = base * base % mod;
exp >>= 1;
}
return result;
}
/**
* 判断是否为质数
*/
private boolean isPrime(int n) {
if (n < 2) return false;
if (n == 2) return true;
if (n % 2 == 0) return false;
for (int i = 3; i * i <= n; i += 2) {
if (n % i == 0) return false;
}
return true;
}
/**
* 计算小组合数 C(n, m) mod p,其中 0 <= n, m < p
*/
private long smallC(int n, int m) {
if (m < 0 || m > n) return 0;
// C(n,m) = n! * inv(m!) * inv((n-m)!)
return fact[n] * invFact[m] % p * invFact[n - m] % p;
}
/**
* 卢卡斯定理:计算 C(n, m) mod p
* @param n 上标(可大到10^18)
* @param m 下标(可大到10^18)
* @return C(n, m) mod p
*/
public long lucas(long n, long m) {
if (m < 0 || m > n) return 0;
long result = 1;
// 按p进制分解,逐位计算
while (n > 0 || m > 0) {
int ni = (int)(n % p);
int mi = (int)(m % p);
if (mi > ni) return 0; // 该位无法选取,组合数为0
result = result * smallC(ni, mi) % p;
n /= p;
m /= p;
}
return result;
}
/**
* 扩展版本:处理多个质数模数的CRT合并(中国剩余定理)
* 用于计算 C(n,m) mod M,其中M为若干不同质数的乘积
*/
public static class ExtendedLucas {
private final int[] primes;
private final LucasTheorem[] lucasInstances;
public ExtendedLucas(int[] primes) {
this.primes = primes.clone();
this.lucasInstances = new LucasTheorem[primes.length];
for (int i = 0; i < primes.length; i++) {
lucasInstances[i] = new LucasTheorem(primes[i]);
}
}
/**
* 中国剩余定理合并多个同余方程
* x ≡ a_i (mod p_i)
*/
public long crtMerge(long[] remainders) {
long result = 0;
long M = 1;
for (int prime : primes) {
M *= prime;
}
for (int i = 0; i < primes.length; i++) {
long Mi = M / primes[i];
// 求 Mi 在模 primes[i] 下的逆元
long invMi = modInverse(Mi % primes[i], primes[i]);
result = (result + remainders[i] * Mi % M * invMi) % M;
}
return (result % M + M) % M;
}
private long modInverse(long a, long mod) {
return modPow(a, mod - 2, mod);
}
private long modPow(long base, long exp, long mod) {
long result = 1;
base %= mod;
while (exp > 0) {
if ((exp & 1) == 1) {
result = result * base % mod;
}
base = base * base % mod;
exp >>= 1;
}
return result;
}
/**
* 计算 C(n,m) mod (p1*p2*...*pk)
*/
public long solve(long n, long m) {
long[] rems = new long[primes.length];
for (int i = 0; i < primes.length; i++) {
rems[i] = lucasInstances[i].lucas(n, m);
}
return crtMerge(rems);
}
}
public static void main(String[] args) {
// 示例1:基础卢卡斯定理
System.out.println("=== Lucas Theorem Demo ===");
int p = 1009; // 质数模数
LucasTheorem lucas = new LucasTheorem(p);
// 测试用例:大组合数取模
long[][] testCases = {
{10, 3}, // 小数值验证
{100, 50}, // 中等数值
{1000000000L, 500000000L}, // 超大数值
{1000000000000L, 500000000000L}, // 极大数值
};
for (long[] tc : testCases) {
long n = tc[0], m = tc[1];
long result = lucas.lucas(n, m);
System.out.printf("C(%d, %d) mod %d = %d%n", n, m, p, result);
}
// 示例2:验证正确性(小数值对比直接计算)
System.out.println("\n=== Verification for small values ===");
LucasTheorem smallLucas = new LucasTheorem(97);
for (int n = 0; n <= 20; n++) {
for (int m = 0; m <= n; m++) {
long direct = directC(n, m) % 97;
long lucasResult = smallLucas.lucas(n, m);
if (direct != lucasResult) {
System.out.printf("MISMATCH: C(%d,%d) direct=%d, lucas=%d%n", n, m, direct, lucasResult);
}
}
}
System.out.println("Small value verification passed!");
// 示例3:扩展卢卡斯 + 中国剩余定理
System.out.println("\n=== Extended Lucas with CRT ===");
int[] primeMods = {3, 5, 7, 11}; // 互质模数,合并后M=1155
ExtendedLucas extLucas = new ExtendedLucas(primeMods);
long n = 100, m = 50;
long result = extLucas.solve(n, m);
System.out.printf("C(%d, %d) mod %d = %d%n", n, m, 1155, result);
// 验证:分别检查模每个质数的结果
System.out.println("Verification per modulus:");
for (int prime : primeMods) {
long subResult = extLucas.lucasInstances[0].lucas(n, m); // 使用对应实例
// 重新计算以展示
LucasTheorem lt = new LucasTheorem(prime);
System.out.printf(" mod %d: C(%d,%d) = %d%n", prime, n, m, lt.lucas(n, m));
}
// 性能测试
System.out.println("\n=== Performance Test ===");
LucasTheorem perfLucas = new LucasTheorem(1000003);
long start = System.currentTimeMillis();
for (int i = 0; i < 100000; i++) {
perfLucas.lucas(1000000000000L + i, 500000000000L + i);
}
long elapsed = System.currentTimeMillis() - start;
System.out.println("100000 queries: " + elapsed + "ms");
System.out.println("Average: " + (elapsed / 100000.0) + "ms per query");
}
// 辅助方法:小数值直接计算组合数用于验证
private static long directC(int n, int m) {
if (m < 0 || m > n) return 0;
long result = 1;
for (int i = 1; i <= m; i++) {
result = result * (n - m + i) / i;
}
return result;
}
}
关键要点解析
1. 阶乘逆元的递推求法
除了使用费马小定理对每个阶乘单独求逆元外,更优的方法是:
– 先求 $invFact[p-1] = (p-1)!^{p-2} \mod p$(费马小定理)
– 再利用 $invFact[i-1] = invFact[i] \cdot i \mod p$ 递推
这样将逆元预处理优化到 $O(p)$。
2. p进制分解的终止条件
在 lucas() 方法中,循环条件为 while (n > 0 || m > 0)。若只写 n > 0,当 $n < m$ 且 $n$ 先变为0时可能遗漏判断;确保 $m$ 的每一位都被检查是关键。
3. 模数非质数的处理
卢卡斯定理要求模数为质数。若模数 $M$ 可分解为不同质因子的乘积 $M = p_1 \cdot p_2 \cdots p_k$,可分别对每个质因子使用卢卡斯定理,再通过中国剩余定理(CRT)合并结果。
复杂度分析
| 步骤 | 时间复杂度 | 空间复杂度 |
|---|---|---|
| 预处理(阶乘+逆元) | $O(p)$ | $O(p)$ |
| 单次卢卡斯查询 | $O(\log_p n)$ | $O(1)$ |
| 扩展版本(k个质因子) | $O(k \cdot \log_p n)$ | $O(k \cdot p)$ |
应用场景
- 组合数大数取模:$n, m$ 达到 $10^{18}$ 级别时的模运算
- 二项式系数模质数:概率计算、多项式展开中的系数提取
- 密码学基础:RSA、椭圆曲线等算法中的模运算组件
- 组合计数问题:在模意义下计算方案数,如路径计数、选取方案等
卢卡斯定理将”不可能完成的大数计算”转化为”可控的小规模迭代”,是数论分治思想的典型范例。