每日算法 — 使用java实现卢卡斯定理:组合数取模与数论分治

算法概述

卢卡斯定理(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、椭圆曲线等算法中的模运算组件
  • 组合计数问题:在模意义下计算方案数,如路径计数、选取方案等

卢卡斯定理将”不可能完成的大数计算”转化为”可控的小规模迭代”,是数论分治思想的典型范例。