数组合数学

组合数

字数统计: 4k字   |   阅读时长: 20分
, , ,

摘要: 求组合数的代码模板和例题

【对算法，数学，计算机感兴趣的同学，欢迎关注我哈，阅读更多原创文章】
我的网站：潮汐朝夕的生活实验室
我的公众号：算法题刷刷
我的知乎：潮汐朝夕
我的github：FennelDumplings
我的leetcode：FennelDumplings

在文章 组合数学1-排列组合 中我们详细梳理了组合数学的基本模型，以及排列组合的相关公式。本文我们研究编程求组合数的几个算法和代码模板。主要内容如下：

  • 组合数的公式: $\dbinom{n}{m} = \frac{n!}{m!(n-m)!}$
  • 性质
    1. $\dbinom{n}{m} = \dbinom{n}{n-m}$
    2. $\dbinom{n}{m} = \dbinom{n-1}{m} + \dbinom{n-1}{m-1}$
    3. $\dbinom{n}{0} + \dbinom{n}{1} + … + \dbinom{n}{n} = 2^{n}$
  • 组合数$\dbinom{n}{m}$求法及模板
    • (a). 杨辉三角：利用性质2，递推法。求 $0 \leq y \leq x \leq n$ 的所有 $\dbinom{x}{y}$, 时间复杂度 $O(N^{2})$
    • (b). 利用 $ln(ab) = ln(a) + ln(b)$ 做中转
  • 大组合数模质数 $\dbinom{n}{m} \mod p$求法及模板
    • (c). 1~n 均存在模 p 乘法逆元，可以先求 $n! \mod p$，再求 $m!(n-m)! \mod p$ 的逆元，再乘起来，时间复杂度 $O(N)$，逆元的写法：同余
    • (d). 将每个$1 \leq k \leq n$ 的 $k! \mod p$ 及其逆元都预处理出来，然后 $O(1)$ 响应查询。
  • 大组合数 $\dbinom{n}{m}$ 的高精度求法及模板
    • (e). 对 $\dbinom{n}{m}$ 进行高精度计算，高精度的写法：高精度运算

$1 求组合数

求 $\dbinom{n}{m}$ 的几个方法。

1.1 递推法(杨辉三角)

二维数组 c 计算完之后可以多次查询。

时间复杂度 $O(NM)$

1
2
3
4
5
6
7
8
9
10
11
12
13
ll C(int n, int m)
{
    // 调用方保证 n >= m
    if(m > n - m)
        return C(n, n - m);
    vector<vector<ll>> c(n + 1, vector<ll>(m + 1, (ll)0));
    for(int i = 0; i <= n; ++i)
        c[i][0] = 1;
    for(int i = 1; i <= n; ++i)
        for(int j = 1; j <= min(i, m); ++j)
            c[i][j] = (c[i - 1][j] + c[i - 1][j - 1]);
    return c[n][m];
}

取模的版本：

1
2
3
4
5
6
7
8
9
10
11
12
13
ll C(int n, int m, int MOD)
{
    // 调用方保证 n >= m
    if(m > n - m)
        return C(n, n - m, MOD);
    vector<vector<ll>> c(n + 1, vector<ll>(m + 1, (ll)0));
    for(int i = 0; i <= n; ++i)
        c[i][0] = 1;
    for(int i = 1; i <= n; ++i)
        for(int j = 1; j <= min(i, m); ++j)
            c[i][j] = (c[i - 1][j] + c[i - 1][j - 1]) % MOD;
    return c[n][m];
}

1.2 用对数做中转

基于杨辉三角的递推法有两个问题，一个是很容易溢出，另一个是时间复杂度比较弱。

取对数做中转可以解决溢出的问题。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
using ll = long long;

long double lnC(int n, int m)
{
    if (m > n)
    {
        return 0;
    }
    if (m * 2 < n)
    {
        m = n - m;
    }
    long double s1 = 0;
    for (int i = m + 1; i <= n; ++i)
    {
        s1 += log((long double)i);
    }
    long double s2 = 0;
    int ub = n - m;
    for (int i = 2; i <= ub; ++i)
    {
        s2 += log((long double)i);
    }
    return s1 - s2;
}

ll C(int n, int m)
{
    if (m > n)
    {
        return 0;
    }
    return round(exp(lnC(n, m)));
}

$2 求大组合数模质数

求 $\dbinom{n}{m} \mod p$ 的几个方法。

2.1 逆元法

(1) 用 $m!(n-m)! \mod p$ 的逆元

1~n 均存在模 p 乘法逆元，可以先求 $n! \mod p$，再求 $m!(n-m)! \mod p$ 的逆元，再乘起来，时间复杂度 $O(N)$。

求逆元直接使用代码模板，接口如下。算法细节参考 同余。

1
2
// 求 b 模 p 的乘法逆元 (b 与 p 互质)
ll inv(ll b, ll p)

代码模板

没有给 p 的时候可以尝试取 p 为 INT 范围最大的质数 2147483629。下面代码中的 inv() 是求逆元的函数，见上面的代码模板。

时间复杂度 $O(N)$。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
// 求 n! mod p
ll fac(ll n, ll p)
{
    if(n == 0)
        return 1;
    ll ans = 1;
    for(int i = 1; i <= n; ++i)
        ans = (ans * i) % p;
    return ans;
}

ll C(ll n, ll m, ll p)
{
    if(m == 0) return 1;
    // 先求 $n! \mod p$
    ll ans1 = fac(n, p);
    // 再求 $m!(n-m)! \mod p$ 的逆元
    ll ans2 = (fac(m, p) * fac(n - m, p)) % p;
    ans2 = inv(ans2, p);
    // 再乘起来
    return (ans1 * ans2) % p;
}

(2) 用 $m \mod p$ 的逆元

直接求 $P(n, m) = n(n - 1)\cdots(n - (m - 1))$，再乘以 $m!$ 的乘法逆元。

而不计算阶乘。

代码模板

1
2
3
4
5
6
7
8
9
10
11
int C(ll n, int m, int p)
{
    // n < 0 和 n < m 的情况在调用方处理
    n %= p;
    int ans = 1;
    for(int i = 0; i < m; ++i)
        ans = (ll)ans * (n - i) % p;
    for(int i = 1; i <= m; ++i)
        ans = (ll)ans * inv(i, p) % p;
    return ans;
}

如果 m 的范围不大，也可以预处理 m % p 的逆元，比预处理阶乘的逆元优秀。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
// 预处理逆元

vector<int> invs;

void init(int N, int p)
{
    invs.assign(N, -1);
    for(int i = 1; i < N; ++i)
        invs[i] = inv(i, p);
}

int C(ll n, int m, int p)
{
    // n < 0 和 n < m 的情况在调用方处理
    n %= p;
    int ans = 1;
    for(int i = 0; i < m; ++i)
        ans = (ll)ans * (n - i) % p;
    for(int i = 1; i <= m; ++i)
        ans = (ll)ans * invs[i] % p;
    return ans;
}

2.2 预处理阶乘和逆元

逆元法的过程中，将 0 <= k <= n 的阶乘 k! % p 及其逆元存在 jc 和 jc_inv 中，可以在 $O(N\log N)$ 预处理后，$O(1)$ 响应 y <= x <= n 的查询 C(n, m, p) = (jc[n] * jx_inv[m] * jc_inv[n - m]) % p。

代码模板

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
using ll = long long;

vector<ll> jc, jc_inv;

void preprocess(ll n, ll p)
{
    // [1, n] 的阶乘及其逆元
    // jc[i] = i! % p
    jc = vector<ll>(n + 1, 1);
    for(ll i = 1; i <= n; ++i)
        jc[i] = (jc[i - 1] * i) % p;
    jc_inv = vector<ll>(n + 1, 1);
    // 0 不存在逆元
    for(ll i = 1; i <= n; ++i)
        jc_inv[i] = inv(jc[i], p);
}

ll C(ll n, ll m, ll p)
{
    // 先求 $n! \mod p$
    ll ans1 = jc[n];
    // 再求 $m!(n-m)! \mod p$ 的逆元
    ll ans2 = (jc_inv[m] * jc_inv[n - m]) % p;
    // 再乘起来
    return (ans1 * ans2) % p;
}

其中 jc_inv[i] = inv(jc[i], p); 还可以换成下面这种递推的形式：

1
jc_inv[i] = (ll)jc_inv[i - 1] * inv(i, p) % p;

$3 求高精度大组合数

方案1：直接用支持除法的高精度

struct BigInteger 为代码模板，算法细节以及实现见 高精度运算。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
BigInteger fac(BigInteger n)
{
    if(n == 0)
        return 1;
    BigInteger ans = 1;
    for(int i = 1; i <= n; ++i)
        ans = (ans * i);
    return ans;
}

BigInteger C(BigInteger n, BigInteger m)
{
    if(m == 0) return 1;
    BigInteger ans1 = fac(n);
    BigInteger ans2 = (fac(m) * fac(n - m));
    return ans1 / ans2;
}

方案2：避免除法

用 阶乘分解 的算法，把分子，分母分解质因数，在数组中保存各项质因子的指数，然后把分子分母各个质因子的指数对应相减，再把剩余质因子乘起来。

$4 题目: 62. 不同路径

一个机器人位于一个 m x n 网格的左上角 （起始点在下图中标记为 “Start” ）。

机器人每次只能向下或者向右移动一步。机器人试图达到网格的右下角（在下图中标记为 “Finish” ）。

问总共有多少条不同的路径？

提示：

1
2
1 <= m, n <= 100
题目数据保证答案小于等于 2 * 1e9

示例 1：
输入：m = 3, n = 7
输出：28

示例 2：
输入：m = 3, n = 2
输出：3
解释：
从左上角开始，总共有 3 条路径可以到达右下角。

  1. 向右 -> 向下 -> 向下
  2. 向下 -> 向下 -> 向右
  3. 向下 -> 向右 -> 向下

示例 3：
输入：m = 7, n = 3
输出：28

示例 4：
输入：m = 3, n = 3
输出：6

写出组合数

组合计数结果如下。

下面的算法 1 是与组合数无关的动态规划算法。算法 2 ~ 6 分别为前面介绍的求组合数的方法，直接用代码模板。

算法1. 动态规划

1
2
3
4
5
6
7
8
9
10
11
12
13
14
class Solution {
public:
    int uniquePaths(int m, int n) {
        vector<vector<int> > dp(m, vector<int>(n, 0));
        dp[0][0] = 1;
        for(int i = 0; i < m; ++i)
            for(int j = 0; j < n; ++j)
            {
                if(i > 0) dp[i][j] += dp[i - 1][j];
                if(j > 0) dp[i][j] += dp[i][j - 1];
            }
        return dp[m - 1][n - 1];
    }
};

算法2. 组合数-杨辉三角法

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
ll C(int n, int m)
{
    // 调用方保证 n >= m
    if(m > n - m)
        return C(n, n - m);
    vector<vector<ll>> c(n + 1, vector<ll>(m + 1, (ll)0));
    for(int i = 0; i <= n; ++i)
        c[i][0] = 1;
    for(int i = 1; i <= n; ++i)
        for(int j = 1; j <= min(i, m); ++j)
            c[i][j] = (c[i - 1][j] + c[i - 1][j - 1]);
    return c[n][m];
}

class Solution {
public:
    int uniquePaths(int m, int n) {
        return C(n + m - 2, n - 1);
    }
};

算法3. 组合数-杨辉三角法,对数作中转

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
long double lnC(int n, int m)
{
    if (m > n)
    {
        return 0;
    }
    if (m * 2 < n)
    {
        m = n - m;
    }
    long double s1 = 0;
    for (int i = m + 1; i <= n; ++i)
    {
        s1 += log((long double)i);
    }
    long double s2 = 0;
    int ub = n - m;
    for (int i = 2; i <= ub; ++i)
    {
        s2 += log((long double)i);
    }
    return s1 - s2;
}

ll C(int n, int m)
{
    if (m > n)
    {
        return 0;
    }
    return round(exp(lnC(n, m)));
}

class Solution {
public:
    int uniquePaths(int m, int n) {
        return C(n + m - 2, n - 1);
    }
};

算法4. 逆元法

代码中 inv(b, p) 为求 b 模 p 的逆元。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
using ll = long long;

// 求 n! mod p
ll fac(ll n, ll p)
{
    if(n == 0)
        return 1;
    ll ans = 1;
    for(int i = 1; i <= n; ++i)
        ans = (ans * i) % p;
    return ans;
}

ll C(ll n, ll m, ll p)
{
    if(m == 0) return 1;
    // 先求 $n! \mod p$
    ll ans1 = fac(n, p);
    // 再求 $m!(n-m)! \mod p$ 的逆元
    ll ans2 = (fac(m, p) * fac(n - m, p)) % p;
    ans2 = inv(ans2, p);
    // 再乘起来
    return (ans1 * ans2) % p;
}

// INT 范围最大的质数 2147483629
class Solution {
public:
    int uniquePaths(int m, int n) {
        const ll MOD = 2147483629;
        return C(n + m - 2, n - 1, MOD);
    }
};

算法5. 逆元+预处理

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
// INT 范围最大的质数 2147483629
class Solution {
public:
    int uniquePaths(int m, int n) {
        const ll MOD = 2147483629;
        preprocess(n + m - 2, MOD);
        return C(n + m - 2, n - 1, MOD);
    }

private:
    vector<ll> jc;
    vector<ll> jc_inv;

    void preprocess(ll n, ll p)
    {
        // [1, n] 的阶乘及其逆元
        // jc[i] = i! % p
        jc = vector<ll>(n + 1, 1);
        for(ll i = 1; i <= n; ++i)
            jc[i] = (jc[i - 1] * i) % p;
        jc_inv = vector<ll>(n + 1, 1);
        // 0 不存在逆元
        for(ll i = 1; i <= n; ++i)
            jc_inv[i] = inv(jc[i], p);
    }

    ll C(ll n, ll m, ll p)
    {
        if(m == 0) return 1;
        // 先求 $n! \mod p$
        ll ans1 = jc[n];
        // 再求 $m!(n-m)! \mod p$ 的逆元
        ll ans2 = (jc_inv[m] * jc_inv[n - m]) % p;
        // 再乘起来
        return (ans1 * ans2) % p;
    }
};

算法6. 高精度

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
BigInteger fac(BigInteger n)
{
    if(n == 0)
        return 1;
    BigInteger ans = 1;
    for(BigInteger i = 1; i <= n; i += 1)
        ans = (ans * i);
    return ans;
}

BigInteger C(BigInteger n, BigInteger m)
{
    if(m == 0) return 1;
    BigInteger ans1 = fac(n);
    BigInteger ans2 = (fac(m) * fac(n - m));
    return ans1 / ans2;
}


class Solution {
public:
    int uniquePaths(int m, int n) {
        BigInteger ans_big = C(n + m - 2, n - 1);
        stringstream ss;
        ss << ans_big;
        int ans_int;
        ss >> ans_int;
        return ans_int;
    }
};

$5 卢卡斯定理

若 p 是质数，对于任意整数 $1 \leq m \leq n$，有：

即把 m, n 表示成 p 进制数，对 p 进制下的每一位分别计算组合数，最后再乘起来。证明需要用到母函数理论。

1
Lucas(n, m, p) = C(n % p, m % p, p) * Lucas(n / p, m / p, p)

卢卡斯定理求 $\dbinom{n}{m} \mod p$ 代码模板

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
// 卢卡斯定理模板
// Lucas 定理实现 C(n,m) % p 的代码：
using ll = long long;
ll exp_mod(ll a, ll b, ll p)
{
    //快速幂取模
    ll res = 1;
    while(b != 0)
    {
        if(b & 1) res = (res * a) % p;
        a = (a * a) % p;
        b >>= 1;
    }
    return res;
}

ll Comb(ll a, ll b, ll p)
{
    //求组合数C(a,b) % p
    if(a < b)   return 0;
    if(a == b)  return 1;
    if(b > a - b)   b = a - b;

    ll ans = 1, ca = 1, cb = 1;
    for(ll i = 0; i < b; ++i)
    {
        ca = (ca * (a - i)) % p;
        cb = (cb * (b - i)) % p;
    }
    ans = (ca * exp_mod(cb, p - 2, p)) % p;
    return ans;
}

ll Lucas(ll n, ll m, ll p)
{
    //Lucas定理求C(n,m)%p
    ll ans = 1;

    while(n && m && ans)
    {
        ans = (ans * Comb(n % p, m % p, p)) % p;
        n /= p;
        m /= p;
    }
    return ans;
}

Lucas定理 + 预处理 (代码模板)

组合数用【逆元 + 预处理】来做，这部分需要预处理阶乘及其逆元，其中所需逆元用递推法来预处理。

预处理过程: inv[i] = -(p/i)*inv[i%p]，边界 inv[1] = 1，具体可以参考 同余。

此方法适用于 p 为不大的素数且需要多次调用的情况。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
vector<int> fac, invfac;
vector<int> inv;

void preprocess(ll p)
{
    // 预处理 [1, p) 的阶乘及其逆元
    fac.assign(p, 1);
    invfac.assign(p, 1);
    for(int i = 2; i < p; ++i)
    {
        fac[i] = ((ll)fac[i - 1] * i) % p;
        invfac[i] = (ll)invfac[i - 1] * inv[i] % p;
    }
}

void get_inv(ll p)
{
    // 预处理 [1, p) 的逆元
    inv.assign(p, 1);
    for(int i = 2; i < p; ++i)
        inv[i] = p - ((p / i) * inv[p % i]) % p;
}

int Comb(int n, int m, int p)
{
    if(n < 0 || m < 0 || n < m)
        return 0;
    return (ll)fac[n] * invfac[m] % p * invfac[n - m] % p;
}

ll Lucas(ll n, ll m, ll p)
{
    // Lucas定理求C(n,m)%p
    ll ans = 1;

    while(n && m && ans)
    {
        ans = (ans * Comb(n % p, m % p, p)) % p;
        n /= p;
        m /= p;
    }
    return ans;
}

Lucas定理例题: 1259. 不相交的握手

偶数 个人站成一个圆，总人数为 num_people 。每个人与除自己外的一个人握手，所以总共会有 num_people / 2 次握手。

将握手的人之间连线，请你返回连线不会相交的握手方案数。

由于结果可能会很大，请你返回答案 模 1e9+7 后的结果。

提示：

1
2
2 <= num_people <= 1000
num_people % 2 == 0

示例 1：
输入：num_people = 2
输出：1

示例 2：
输入：num_people = 4
输出：2
解释：总共有两种方案，第一种方案是 [(1,2),(3,4)] ，第二种方案是 [(2,3),(4,1)] 。

示例 3：
输入：num_people = 6
输出：5

示例 4：
输入：num_people = 8
输出：14

写出组合数

组合计数结果如下，其中 p = 1e9 + 7。

下面的算法 1 和 2 求大组合数模质数的方法，算法 3 是卢卡斯定理，直接用代码模板。

算法1. 杨辉三角

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
using ll = long long;

ll C(int n, int m, int MOD)
{
    // 调用方保证 n >= m
    if(m > n - m)
        return C(n, n - m, MOD);
    vector<vector<ll>> c(n + 1, vector<ll>(m + 1, (ll)0));
    for(int i = 0; i <= n; ++i)
        c[i][0] = 1;
    for(int i = 1; i <= n; ++i)
        for(int j = 1; j <= min(i, m); ++j)
            c[i][j] = (c[i - 1][j] + c[i - 1][j - 1]) % MOD;
    return c[n][m];
}

class Solution {
public:
    int numberOfWays(int num_people) {
        using ll = long long;
        int MOD = 1e9 + 7;
        if(num_people == 0) return 0;
        int n = num_people / 2;
        ll ans = (C(2 * n, n, MOD) - C(2 * n, n - 1, MOD) + MOD) % MOD;
        return (int)ans;
    }
};

算法2. 逆元法

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
// 求 n! mod p
ll fac(ll n, ll p)
{
    if(n == 0)
        return 1;
    ll ans = 1;
    for(int i = 1; i <= n; ++i)
        ans = (ans * i) % p;
    return ans;
}

ll C(ll n, ll m, ll p)
{
    if(m == 0) return 1;
    // 先求 $n! \mod p$
    ll ans1 = fac(n, p);
    // 再求 $m!(n-m)! \mod p$ 的逆元
    ll ans2 = (fac(m, p) * fac(n - m, p)) % p;
    ans2 = inv(ans2, p);
    // 再乘起来
    return (ans1 * ans2) % p;
}

class Solution {
public:
    int numberOfWays(int num_people) {
        using ll = long long;
        int MOD = 1e9 + 7;
        if(num_people == 0) return 0;
        int n = num_people / 2;
        ll ans = (C(2 * n, n, MOD) - C(2 * n, n - 1, MOD) + MOD) % MOD;
        return (int)ans;
    }
};

算法3. 卢卡斯定理

1
2
3
4
5
6
7
8
9
10
11
class Solution {
public:
    int numberOfWays(int num_people) {
        using ll = long long;
        int MOD = 1e9 + 7;
        if(num_people == 0) return 0;
        int n = num_people / 2;
        ll ans = (Lucas((ll)2 * n, n, (ll)MOD) - Lucas((ll)2 * n, (ll)(n - 1), (ll)MOD) + MOD) % MOD;
        return (int)ans;
    }
};
Share