CC++ & Algorithm

模逆元在组合数计算中的应用

极难3
语言版本：通用
概述：用生活例子理解模逆元，学会用C++和Python高效计算大组合数取模。

模逆元大揭秘：轻松计算大组合数取模

你遇到过这样的问题吗？从全班 100 个同学中选 5 个去参加比赛，有多少种选法？答案是 C(100,5) = 75,287,520。这个数字还不算太大。但如果要从 100 万件商品中选 50 万件，结果会大得惊人，直接计算会溢出，而且我们往往只需要结果除以一个大质数（比如 10^9+7）的余数。这时，模逆元 就是一把神奇的钥匙，它能把复杂的除法转换成简单的乘法，让我们又快又准地算出组合数取模的结果。

什么是模逆元？—— 一个“分蛋糕”的故事

想象你有一个大蛋糕，要把蛋糕分给 5 个小朋友，每人分到的块数必须相等。如果蛋糕有 20 块，每人得到 20 ÷ 5 = 4 块。现在蛋糕变成了一个“模世界”，我们只关心每块蛋糕的编号除以 7 的余数。假设你知道蛋糕有 20 块（在模 7 下是 6），但不知道具体块数，你想知道每人分到多少块（模 7 下）。20 ÷ 5 = 4，而 4 模 7 还是 4。如果蛋糕是 13 块（模 7 下也是 6），13 ÷ 5 不是整数，在模世界里除法不能直接进行。这时就需要一个特殊的规则：除以 5 等于乘以 5 的逆元。

在模 p 的算术中，除法不直接进行，因为模 p 下的整数不一定可除（除非除数与模数互质）。模逆元就是用来把除法转换成乘法的一个工具。比如在模 7 下，你需要找到一个数 x，使得 5 × x ≡ 1 (mod 7)，这个 x 就是 5 的逆元。因为 5 × 3 = 15 ≡ 1 (mod 7)，所以逆元是 3。那么 13 ÷ 5 相当于 13 × 3 = 39 ≡ 4 (mod 7)，正好等于 (13 mod 7 = 6) × 3 = 18 ≡ 4。你看，逆元使得除法变成了乘法。

生活中还有哪些类似情况？

  • 假如你有一包零食，想平均分给 4 个朋友，但零食的总数模 11 后为 7（即不知道具体总数，只知道除以 11 余 7），问每人分到多少（模 11 下）？需要先找到 4 的逆元（4 × 3 = 12 ≡ 1 mod 11，所以逆元是 3），然后 7 × 3 = 21 ≡ 10 (mod 11)，每人分到 10（模 11）。
  • 再比如考试排名：如果班级有 30 人，你考了第 6 名，你想知道名次除以 7 的余数是多少？6 ÷ 7 余 6，没问题。但如果你想知道排名在模 7 下是第几名，而排名是 13，13 ÷ 7 不能整除，就需要逆元来帮忙啦。

组合数与模逆元的关系

组合数 C(n,k) 表示从 n 个物品中选 k 个的方法数，计算公式是：

C(n,k)=n!k! (n−k)!C(n,k) = \frac{n!}{k!\,(n-k)!}

当 n 和 k 很大（比如 n = 10^6）并且结果需要模一个素数 p（比如 p = 10^9+7）时，直接计算阶乘再除法会遇到大整数溢出和除法不整除的问题。因为模 p 下，除法必须用逆元代替。所以我们会把公式变成：

C(n,k)mod  p=n!×(k!)−1×((n−k)!)−1mod  pC(n,k) \mod p = n! \times (k!)^{-1} \times ((n-k)!)^{-1} \mod p

这里 (x!)−1(x!)^{-1} 表示 x! 的模逆元。因为 p 是素数（通常取 1e9+7），根据费马小定理，任何与 p 互质的数 a 的逆元就是 a^{p-2} mod p。这样我们就可以用快速幂来计算逆元了。

为什么 p 要选素数？
因为费马小定理只适用于素数的模数。如果 p 不是素数，逆元不一定存在，就算存在也得用扩展欧几里得算法求。而在竞赛中常用 1e9+7 这个素数，它足够大，且不是 n 的倍数。

数学原理：费马小定理

费马小定理说：如果 p 是素数，a 是整数且不是 p 的倍数，那么

ap−1≡1(modp)a^{p-1} \equiv 1 \pmod{p}

两边同时除以 a（乘以 a 的逆元）得到：

ap−2≡a−1(modp)a^{p-2} \equiv a^{-1} \pmod{p}

所以 a 的逆元就是 a^{p-2} mod p。这要求 p 是素数且 a 不是 p 的倍数，正好符合组合数中分母的阶乘（因为 p > n，所以所有阶乘都不是 p 的倍数）。

巧妙吗？ 用一次快速幂（O(log p)）就能求出逆元，而 p 通常约 10^9，log p 大约 30 次乘法，非常快。

计算步骤（实现思路）

  1. 预处理阶乘：计算 fact[i] = i! mod p，对于 i 从 0 到 n。
  2. 预处理逆元阶乘：有两种方法：
    • 递推法（高效）：先算出 invFact[n] = (n!)^{-1} mod p，然后递推 invFact[i-1] = invFact[i] * i mod p。为什么？因为 (i-1)! × i = i!，两边取倒数得到 (i-1)! 的逆 = i! 的逆 × i。
    • 逐个求逆元法：对每个阶乘单独用快速幂求逆，但 O(n log p) 较慢。
  3. 查询：C(n,k) = fact[n] × invFact[k] % p × invFact[n-k] % p。

这样预处理后，每次查询只需要常数时间 O(1)，非常适合多次询问。

新手容易犯的错误

  1. 忘记预处理：直接调用 combination() 而没调用 precompute()，导致 invFact 全是默认值。
  2. 数组越界：预处理时 MAXN 设置太小，而 n > MAXN。建议 MAXN 设为可能的最大 n + 5。
  3. 模数不是素数时直接使用费马小定理：费马小定理只对素数模有效。如果模数不是素数（如 10^6+7 不是素数），会得出错误结果。此时必须用扩展欧几里得求逆元。
  4. 对不互质的数求逆元：如果 a 和 p 不互质（比如 a 是 p 的倍数），逆元不存在。但在组合数中分母阶乘 mod p 不为 0（因为 p > n），所以安全。
  5. 乘法溢出：在 C++ 中，两个 int 相乘可能溢出，应使用 long long 并取模。
  6. 快速幂中的指数类型：exp 可能是 long long，但 p-2 在 p=1e9+7 范围内，int 足够，但为了安全用 long long。

例子：如果 MAXN 设成了 100，但实际 n=1000，那么 precompute 不会计算到 1000，combination 时 fact[1000] 是默认 0，导致结果为 0。所以预处理范围必须覆盖所有查询。

完整可运行代码示例（C++）

下面用 C++ 实现预处理阶乘和逆元，然后计算多次组合数询问。代码中详细注释每一步。

#include <iostream>
#include <vector>
using namespace std;

typedef long long ll;
const int MOD = 1e9 + 7;            // 一个常用素数，大于可能的n
const int MAXN = 1000000;           // 最大的n

// 快速幂：计算 base^exp mod MOD
ll fastPow(ll base, ll exp) {
    ll result = 1;
    while (exp > 0) {
        if (exp & 1)                // 如果exp最低位为1
            result = result * base % MOD;
        base = base * base % MOD;   // base平方
        exp >>= 1;                  // exp右移一位
    }
    return result;
}

// 预处理阶乘和逆元阶乘
vector<ll> fact(MAXN+1);            // 阶乘数组，fact[i] = i! mod MOD
vector<ll> invFact(MAXN+1);         // 逆元阶乘数组，invFact[i] = (i!)^(-1) mod MOD

void precompute() {
    // 阶乘
    fact[0] = 1;
    for (int i = 1; i <= MAXN; i++) {
        fact[i] = fact[i-1] * i % MOD;
    }
    // 逆元阶乘：先计算最大值的逆元
    invFact[MAXN] = fastPow(fact[MAXN], MOD-2); // 费马小定理求逆元
    // 然后递推得到所有逆元阶乘
    for (int i = MAXN; i >= 1; i--) {
        invFact[i-1] = invFact[i] * i % MOD;    // 因为 (i-1)! * i = i!
                                                 // 所以 (i-1)!的逆 = i!的逆 * i
    }
}

// 计算组合数 C(n, k) mod MOD
ll combination(int n, int k) {
    if (k < 0 || k > n) return 0;              // 不合法情况
    return fact[n] * invFact[k] % MOD * invFact[n-k] % MOD;
}

int main() {
    precompute();                              // 预处理！别忘了
    int n = 10, k = 3;
    cout << "C(" << n << "," << k << ") mod " << MOD << " = " << combination(n,k) << endl;
    // 输出应为 120 mod 1e9+7 = 120
    n = 1000; k = 500;
    cout << "C(" << n << "," << k << ") mod " << MOD << " = " << combination(n,k) << endl;

    // 生活例子：从30个同学中选5个，有多少种选法（取模）
    n = 30; k = 5;
    cout << "从30个同学里选5个的方法数 mod " << MOD << " = " << combination(n,k) << endl;
    // 实际 C(30,5)=142506，取模后仍是142506
    return 0;
}

代码解释：

  • fastPow 实现了快速幂，用于计算 a^(p-2) 模 p。
  • precompute 首先计算所有阶乘，然后利用最大阶乘的逆元递推计算出所有逆元阶乘。注意递推公式：invFact[i-1] = invFact[i] * i % MOD，这是因为 (i-1)! × i = i!，两边取倒数得到 (i-1)!的逆 = i!的逆 × i。
  • combination 直接套用公式，注意乘法可能溢出，但用了 long long 和取模就安全了。
  • 主函数测试了两个例子，并加了一个生活场景。

运行结果：

C(10,3) mod 1000000007 = 120
C(1000,500) mod 1000000007 = 330810883  (这个数很大，但计算机能算)
从30个同学里选5个的方法数 mod 1000000007 = 142506

完整可运行代码示例（Python）

Python 自带大整数，但取模运算同样需要逆元。Python 版本使用与 C++ 完全相同的算法。

MOD = 10**9 + 7
MAXN = 1000000

def fast_pow(base, exp):
    """快速幂，返回 base^exp mod MOD"""
    result = 1
    while exp > 0:
        if exp & 1:
            result = result * base % MOD
        base = base * base % MOD
        exp >>= 1
    return result

# 预处理阶乘和逆元阶乘
fact = [1] * (MAXN + 1)          # fact[i] = i! mod MOD
inv_fact = [1] * (MAXN + 1)      # inv_fact[i] = (i!)^(-1) mod MOD

def precompute():
    # 阶乘
    for i in range(1, MAXN+1):
        fact[i] = fact[i-1] * i % MOD
    # 逆元阶乘
    inv_fact[MAXN] = fast_pow(fact[MAXN], MOD-2)
    for i in range(MAXN, 0, -1):
        inv_fact[i-1] = inv_fact[i] * i % MOD

def combination(n, k):
    """返回 C(n,k) mod MOD"""
    if k < 0 or k > n:
        return 0
    return fact[n] * inv_fact[k] % MOD * inv_fact[n-k] % MOD

# 调用预处理
precompute()

# 测试
print(f"C(10,3) mod {MOD} = {combination(10,3)}")
print(f"C(1000,500) mod {MOD} = {combination(1000,500)}")
print(f"从30个同学里选5个的方法数 mod {MOD} = {combination(30,5)}")

Python 代码解释：

  • 与 C++ 逻辑完全一致，只是语法差异。注意 Python 中 & 是按位与，>>= 是右移。
  • 列表 fact 和 inv_fact 预先分配空间，然后按递推计算。
  • combination 直接返回结果，使用整数乘法取模。

运行结果与 C++ 相同。

进阶：使用扩展欧几里得算法求逆元（备用）

如果模数不是素数，或者不方便使用费马小定理（比如模数很大且不是素数），可以用扩展欧几里得算法求逆元。扩展欧几里得算法可以求解方程 ax + by = gcd(a,b)，当 gcd(a,b)=1 时，x 就是 a 模 b 的逆元。这里我们只给出 C++ 实现备用。

// 扩展欧几里得算法，返回 gcd(a,b)，并设置 x,y 使得 ax+by = gcd
ll exgcd(ll a, ll b, ll &x, ll &y) {
    if (b == 0) { x = 1; y = 0; return a; }
    ll g = exgcd(b, a % b, y, x);
    y -= a / b * x;
    return g;
}

// 求 a 模 MOD 的逆元（假设 MOD 与 a 互质）
ll inv(ll a) {
    ll x, y;
    exgcd(a, MOD, x, y);
    return (x % MOD + MOD) % MOD; // 确保结果为正
}

但通常比赛和实际应用中，费马小定理配合快速幂已经足够，因为常用模数 1e9+7 是素数。

总结

模逆元就像模世界里的“除法按钮”，有了它，我们就能安全地计算组合数取模。核心步骤：预处理阶乘和逆元阶乘（利用费马小定理+快速幂），然后 O(1) 回答每个组合询问。这种方法可以处理 n 高达 10^6 甚至 10^7 的问题，只需要线性时间和常数内存。

希望你现在对模逆元在组合数中的应用有了清晰的理解，赶快动手试试吧！

相关知识点指引

  • [快速幂]：求逆元的基础，学会用 O(log p) 计算 a^b mod p。
  • [扩展欧几里得算法]：当模数不是素数时求逆元的通用方法。
  • [卢卡斯定理]：用于处理 n、k 很大且 p 很小（p 是素数）的组合数取模。
  • [组合数递推]：如果只需要一个组合数，可以用杨辉三角递推，但 n 不能太大。
  • [模运算基础]：理解同余、取模的性质，是学习模逆元的前提。

把这些知识点串起来，你就能轻松应对各种组合数取模的题目啦！

例题精讲

1单选题

在模质数p下，计算组合数C(n,k) mod p时，需要用到模逆元来计算以下哪个表达式？

An! * inv(k!) mod p
Bn! * inv((n-k)!) mod p
Cn! * inv(k! * (n-k)!) mod p
Dn! * inv(k!) * inv((n-k)!) mod p
2判断题

在模数p为质数时，利用费马小定理，可以通过快速幂计算a^(p-2) mod p得到a的逆元，从而用于组合数计算。

3填空题
已知预处理了阶乘数组fact[0..n]和阶乘逆元数组inv_fact[0..n]（inv_fact[i] = (i!)^{-1} mod p），请补全计算组合数C(n,k) mod p的函数。
def comb(n, k, p):
    return fact[n] * inv_fact[k] % p * ___ % p
4单选题

以下关于模逆元在组合数计算中的应用，说法正确的是？

A当模数p为质数时，所有整数都有模逆元
B计算C(n,k) mod p时，需要用到k!和(n-k)!的模逆元
C模逆元可以用费马小定理计算当且仅当p是合数
D使用模逆元可以避免所有除法运算
5填空题
下面是用扩展欧几里得算法求模逆元的代码片段，请补全。
def modinv(a, p):
    g, x, y = extended_gcd(a, p)
    if g != 1:
        return None
    else:
        return ___ % p