不同数据范围下的组合数求解办法。
递推
题目链接：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$
输入样例
输出样例
解法
递推公式：$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$
输入样例
输出样例
解法
公式：$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≤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$
输入样例
输出样例
解法
公式展开：$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;
}
|
参考文献
- AcWing算法基础课
- 常用代码模板4——数学知识
若有疏误，敬请指出。
