组合数求解

不同数据范围下的组合数求解办法。

递推

题目链接:AcWing885 求组合数I

给定 $n$ 组询问,每组询问给定两个整数 $a,b$,请你输出 $C_a^b \ \mathrm{mod} \ 10^9+7$ 的值。

输入格式

第一行包含整数 $n$。

接下来 $n$ 行,每行包含一组 $a$ 和 $b$。

输出格式

共 $n$ 行,每行输出一个询问的解。

数据范围

$1≤n≤10000$, $1≤b≤a≤2000$

输入样例

1
2
3
4
3
3 1
5 3
2 2

输出样例

1
2
3
3
10
1

解法

递推公式:$C_a^b = C_{a-1}^{b} + C_{a-1}^{b-1}$

预处理时间复杂度:$2000 \times 2000 = 4 \times 10^6$;查询:$O(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
28
#include <iostream>

using namespace std;

const int N = 2010, mod = 1e9 + 7;

int n;
int c[N][N];

void init() {
    for (int i = 0; i < N; i++) {
        for (int j = 0; j <= i; j++) {
            if (!j) c[i][j] = 1;
            else c[i][j] = (c[i - 1][j] + c[i - 1][j - 1]) % mod;
        }
    }
}

int main() {
    init();
    cin >> n;
    while (n--) {
        int a, b;
        scanf("%d%d", &a, &b);
        printf("%d\n", c[a][b]);
    }
    return 0;
}

逆元

题目链接:AcWing885 求组和数II

给定 $n$ 组询问,每组询问给定两个整数 $a,b$,请你输出 $C_a^b \ \mathrm{mod} \ 10^9+7$ 的值。

输入格式

第一行包含整数 $n$。

接下来 $n$ 行,每行包含一组 $a$ 和 $b$。

输出格式

共 $n$ 行,每行输出一个询问的解。

数据范围

$1≤n≤10000$, $1≤b≤a≤10^5$

输入样例

1
2
3
4
3
3 1
5 3
2 2

输出样例

1
2
3
3
10
1

解法

公式:$C_a^b = \frac{a!}{b!(a-b)!}$

预处理出 $10^5$ 内阶乘及其逆元模 $p$ 的值。由于 $p$ 为质数,考虑费马小定理,那么 $a$ 的模 $p$ 逆元为 $a^{p-2}\ % \ p$,快速幂求解即可。

数组 $\mathrm{fact}[i]$ 表示 $i$ 的阶乘模 $p$ 的值,$\mathrm{infact}[i]$ 表示 $i$ 的阶乘的逆元模 $p$ 的值

预处理时间复杂度:$10^5 \times \log_2{10^9} = 3 \times 10^6$;查询:$O(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
28
29
30
31
32
33
34
35
36
37
38
39
40
#include <bits/stdc++.h>

using namespace std;

typedef long long LL;

const int N = 100010, mod = 1e9 + 7;

int n;
int fact[N], infact[N];

int ksm(int a, int k) {
    int res = 1;
    while (k) {
        if (k & 1) res = (LL)res * a % mod;
        a = (LL)a * a % mod;
        k >>= 1;
    }
    return res;
}

void init() {
    fact[0] = infact[0] = 1;
    for (int i = 1; i < N; i++) {
        fact[i] = (LL)fact[i - 1] * i % mod;
        infact[i] = (LL)infact[i - 1] * ksm(i, mod - 2) % mod;
        // or: infact[i] = ksm(fact[i], mod - 2);
    }
}

int main() {
    init();
    cin >> n;
    while (n--) {
        int a, b;
        cin >> a >> b;
        printf("%d\n", (LL)fact[a] * infact[b] % mod * infact[a - b] % mod);
    }
    return 0;
}

Lucas定理

题目链接:AcWing885 求组和数III

给定 $n$ 组询问,每组询问给定三个整数 $a,b,p$,其中 $p$ 是质数,请你输出 $C_a^b \ \mathrm{mod} \ p$ 的值。

输入格式

第一行包含整数 $n$。

接下来 $n$ 行,每行包含一组 $a,b,p$。

输出格式

共 $n$ 行,每行输出一个询问的解。

数据范围

$1≤n≤20$, $1≤b≤a≤10^{18}$, $1≤p≤10^5$

输入样例

1
2
3
4
3
5 3 7
3 1 5
6 4 13

输出样例

1
2
3
3
3
2

解法

数据范围:$1≤b≤a≤10^{18}, 1≤p≤10^5$

Lucas定理:对于质数 $p$,则 $1 \le b \le a$ 有:

$$ C_a^b = C_{a \ \mathrm{mod} \ p}^{b \ \mathrm{mod} \ p} \cdot C_{\lfloor a / p \rfloor}^{\lfloor b / p \rfloor} \ (\mathrm{mod} \ p) $$

时间复杂度:$O(\log_pN \cdot p \cdot \log p \cdot n)$,$N=1 \text{~}10^{18}, p=1 \text{~} 10^5, n=20$

 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
#include <bits/stdc++.h>

using namespace std;

typedef long long LL;

int p;

int ksm(int a, int k) {
    int res = 1;
    while (k) {
        if (k & 1) res = (LL)res * a % p;
        a = (LL)a * a % p;
        k >>= 1;
    }
    return res;
}

int C(int a, int b) {
    int res = 1;
    for (int i = 1, j = a; i <= b; i++, j--) {
        res = (LL)res * j % p;
        res = (LL)res * ksm(i, p - 2) % p;
    }
    return res;
}

int lucas(LL a, LL b) {
    if (a < p && b < p) return C(a, b);
    return (LL)C(a % p, b % p) * lucas(a / p, b / p) % p;
}

int main() {
    int n; cin >> n;
    while (n--) {
        LL a, b;
        cin >> a >> b >> p;
        cout << lucas(a, b) << endl;
    }
    return 0;
}

分解质因数

题目链接:AcWing885 求组和数IV

输入 $a,b$,求 $C_a^b$ 的值。

注意结果可能很大,需要使用高精度计算。

输入格式

共一行,包含两个整数 $a$ 和 $b$。

输出格式

共一行,输出 $C_a^b$ 的值。

数据范围

$1≤b≤a≤5000$

输入样例

1
5 3

输出样例

1
10

解法

公式展开:$C_a^b = \frac{a!}{b!(a-b)!}$

主要思路是求出答案中每个质因子 $p$ 的个数 $s$,高精度乘起来。$s$ 的求法可参照 AcWing197 阶乘分解

 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
47
48
49
50
51
52
53
#include <iostream>
#include <algorithm>
#include <vector>

#define vint vector<int>

using namespace std;

const int N = 5010;

int a, b;
int primes[N], cnt;
bool st[N];

void euler(int n) {
    for (int i = 2; i <= n; i++) {
        if (!st[i]) primes[cnt++] = i;
        for (int j = 0; primes[j] <= n / i; j++) {
            st[primes[j] * i] = true;
            if (i % primes[j] == 0) break;
        }
    }
}

vint mul(vint &a, int b) {
    vint c;
    int t = 0;
    for (int i = 0; i < a.size() || t; i++) {
        if (i < a.size()) t += a[i] * b;
        c.push_back(t % 10);
        t /= 10;
    }
    return c;
}

int main() {
    euler(N - 1);
    cin >> a >> b;
    vint res(1, 1);
    for (int i = 0; i < cnt; i++) {
        int p = primes[i];
        int s = 0;
    	// s加上a!中p的个数,减去b!中p的个数,减去(a-b)!中p的个数
        for (int k = a; k; k /= p) s += k / p;
        for (int k = b; k; k /= p) s -= k / p;
        for (int k = a - b; k; k /= p) s -= k / p;
        // 高精度乘法将s个p乘到答案中
        for (int j = 0; j < s; j++) res = mul(res, p);
    }
    for (int i = res.size() - 1; i >= 0; i--) printf("%d", res[i]);
    puts("");
    return 0;
}

参考文献

  1. AcWing算法基础课
  2. 常用代码模板4——数学知识

若有疏误,敬请指出。

updatedupdated2020-05-112020-05-11
加载评论
点击刷新