卢卡斯定理

What · 是什么

卢卡斯定理（Lucas's theorem）：对素数 p 与非负整数 n,m， C(n,m) ≡ C(⌊n/p⌋,⌊m/p⌋) · C(n mod p, m mod p) (mod p)。 其中规定 C(a,b)=0 当 b>a。等价地，把 n,m 写成 p 进制展开，则 C(n,m) ≡ ∏ C(n_i, m_i) (mod p)。由此可在 O(p + T log_p n) 内处理任意大的 n（可远超 64 位），只要模数是小的素数。前提：p 为素数且通常 p ≤ 1e6~1e7（需要预处理 p 以内阶乘逆元）。

Why · 为什么

How · 实现

伪代码

# 预处理 0..p-1 阶乘与逆元（O(p)）
fact[0..p-1], inv_fact[0..p-1]
# 单次小组合数（p 以内）
def comb(a, b):
  if b < 0 or b > a: return 0
  return fact[a] * inv_fact[b] % p * inv_fact[a-b] % p
# Lucas 递归/迭代
def lucas(n, m, p):
  res = 1
  while n or m:
    res = res * comb(n % p, m % p) % p
    n //= p; m //= p
  return res

Python 实现

def lucas_comb(n, m, p):
    """Lucas 定理 + CRT 组合的完整实现。
    p 必须为素数。
    该函数同时验证了与 CRT 结合的经典写法
    （多素因子模数场景请见 extended-lucas 节点）。
    """
    if m < 0 or m > n:
        return 0
    if p == 1:
        return 0

    fact = [1] * p
    for i in range(1, p):
        fact[i] = fact[i - 1] * i % p
    inv_fact = [1] * p
    inv_fact[p - 1] = pow(fact[p - 1], p - 2, p)
    for i in range(p - 2, -1, -1):
        inv_fact[i] = inv_fact[i + 1] * (i + 1) % p

    def comb(a, b):
        if b < 0 or b > a:
            return 0
        return fact[a] * inv_fact[b] % p * inv_fact[a - b] % p

    res = 1
    while n or m:
        res = res * comb(n % p, m % p) % p
        n //= p
        m //= p
    return res


def crt_merge_lucas(n, m, mod):
    """Lucas + CRT：mod 可为合数但其每个素因子都是小素数。
    对每个素因子 p_i 用 Lucas 求余数，再 CRT 合并。
    注意：p_i 若以高次幂出现需用 exLucas（见 extended-lucas 节点）。
    """
    import math
    from math import gcd

    mm = mod
    primes = []
    d = 2
    while d * d <= mm:
        if mm % d == 0:
            primes.append(d)
            while mm % d == 0:
                mm //= d
        d += 1
    if mm > 1:
        primes.append(mm)

    rs = []
    mods = []
    for p in primes:
        if p > 1_000_000:
            raise ValueError("p 过大无法预处理阶乘")
        rs.append(lucas_comb(n, m, p))
        mods.append(p)

    total = 1
    for m in mods:
        total *= m
    x = 0
    for r, m in zip(rs, mods):
        mi = total // m
        x = (x + r * mi * pow(mi, -1, m)) % total
    return x


if __name__ == "__main__":
    import math

    # 与 math.comb 交叉验证（小 p 全范围）
    for p in [2, 3, 5, 7, 11]:
        for n in range(0, 80):
            for m in range(0, n + 1):
                assert lucas_comb(n, m, p) == math.comb(n, m) % p, (n, m, p)

    # 经典样例
    assert lucas_comb(10, 3, 7) == 120 % 7 == 1
    assert lucas_comb(100, 50, 97) == math.comb(100, 50) % 97
    # 大 n（远超 64 位）也能算：n=1e5 的精确组合数仍可由 math.comb 求出并取模交叉验证
    n = 100000
    assert lucas_comb(n, n // 3, 999983) == math.comb(n, n // 3) % 999983
    assert lucas_comb(n, n // 2, 998244353) == math.comb(n, n // 2) % 998244353

    # Lucas+CRT 组合：mod 为 squarefree 合数
    assert crt_merge_lucas(10, 3, 6) == 120 % 6
    assert crt_merge_lucas(100, 40, 35) == math.comb(100, 40) % 35

    print("All tests passed.")

易错点

• p 是素数才可用 Lucas；合数直接套会错（要用 exLucas）
• 预处理阶乘数组要覆盖 0..p-1（大小恰为 p），越界越界访问
• n<m 时 C(n,m)=0，不能直接算（comb 里 b>a 要返回 0）
• 递归写法尾递归深度 O(log_p n)，n 很大时用迭代避免递归栈问题（Python 无尾递归优化）
• CRT 合并时模数互素（素因子）才可用公式，模数为素因子乘积要保证 squarefree
• 阶乘逆元用费马小定理需 p 素数；p=2 时逆元计算仍正确但注意边界

例题

信源

• OI Wiki✓ https://oi-wiki.org/math/number-theory/lucas/
• CP-Algorithms✓ https://cp-algorithms.com/combinatorics/binomial-coefficients.html
• 书籍✓ CPH ch21 Number theory（组合数模素数 p）