acwing算法基础之数学知识--求组合数进阶版

1 基础知识

请明确如下关于取余的基本定理：

数a和数b的乘积模上p，等于数a模上p和数b模上p的乘积。即，
$\cdot b ) \ mod \ p = (a \ mod \ p) \cdot (b \ mod \ p)$
数a除以数b的结果模上p，并不等于数a模上p除以数b模上p。即，
$(a/b)\ mod \ p \neq (a \ mod \ p) / (b \ mod \ p)$

（一）
题目要求：求组合数模上p的结果，即
$C_n^k \ mod\ p =?$
其中 $p = 1 e 9 + 7$ ，是一个质数。

重新考虑组合数 $C_n^k$ 的计算公式，
$C_n^k \ mod \ p=\frac{n!}{k!\cdot (n-k)!} \ mod \ p$
记数 $k!$ 模p的乘法逆元为x，数 $(n - k)!$ 模p的乘法逆元为y，则上式可写成，
$\frac{n!}{k!\cdot (n-k)!} \ mod \ p=n! \cdot x\cdot y \ mod \ p=(n! \ mod \ p)\cdot (x \ mod \ p) \cdot (y \ mod \ p)$
那么，考虑数 $k!$ 模p的乘法逆元x，由于p是质数，故x可由下式计算，
$x \ mod\ p = (k!)^{p-2} \ mod \ p$
观察可以推导出其递推公式，
$(k!)^{p-2} \ mod \ p = ((k-1)!)^{p-2}\cdot k^{p-2} \ mod \ p$
而对于 $k^{p-2}\ mod \ p$ ，可以快速幂在 $O (l o g N)$ 时间复杂度下求解。

故综合上述，可以预处理出阶乘和阶乘的逆元，那么答案可以表示如下，
$\cdot infact[k] \cdot infact[n-k] \ mod \ p$
将上述过程，用代码表述如下，

const int N = 1e5 + 10, mod = 1e9 + 7;
int fact[N], infact[N];int qmi(int a, int k, int p) {long long res = 1;while (k) {if (k & 1) {res = res * a % p;} k >>= 1;a = (long long)a * a % p;}return res;
}void init() {fact[0] = infact[0] = 1;for (int i = 1; i < N; ++i) {fact[i] = (long long)fact[i-1] * i % mod;infact[i] = (long long)infact[i-1] * qmi(i, mod - 2, mod) % mod;}return;
}

（二）
题目要求：
$C_n^{k} \ mod \ p = ？$
其中 $n$ 和 $k$ 的数据范围在 $10^{18}$ 之内，而 $p$ 是质数，且它的范围在 $10^6$ 以内。

对于上述特别大的组合数求解，一般引入Lucas定理。

Lucas定理：当模数p是质数时，有以下等式成立，
$C_n^k \ mod \ p =C_{n \ mod \ p} ^ {k \ mod \ p } \cdot C_{n/p}^{k/p} \ mod \ p$
其中 $k\ mod\ p$ 和 $n\ mod\ p$ 是 $p$ 以内的数，可直接计算组合数 $C_{n\ mod\ p}^{k \ mod \ p}$ ；而对于 $C_{n/p}^{k/p}$ ，则递归使用Lucas定理计算。

故，代码如下所示，

int qmi(int a, int k, int p) {long long res = 1;while (k) {if (k & 1) res = res * a % p;k >>= 1;a = (long long)a * a % p;}return res;
}int C(int a, int b, int p) {if (b > a) return 0; //无效值，返回0long long res = 1;for (int i = 1, j = a; i <= b; ++i, --j) {res = res * j % p;res = res * qmi(i, p - 2, p) % p;}return res;
}int Lucas(long long a, long long b, int p) {if (a < p && b < p) return C(a, b, p); //终止条件return (long long)C(a % p, b % p, p) * Lucas(a / p, b / p, p) % p;
}

（三）
题目要求：
$C_n^k=?$
此处不模上数p了，且 $n$ 和 $k$ 的数据范围在 $10^4$ 以内。

上式可以写成，
$C_n^k=\frac{n!}{k!\cdot (n-k)!}$
考虑任意一个数 $a$ 的阶乘 $a!$ 的分解质因子，
$a!=p_1^{\alpha_1}\cdot p_2^{\alpha_2}\cdots p_k^{\alpha_k}$
对于 $\alpha_i$ ，其中 $0\leq i \leq k$ ，可以通过如下式快速计算，
$\alpha_i=\lfloor \frac{a}{p_i} \rfloor + \lfloor \frac{a}{p_i^2} \rfloor + \lfloor \frac{a}{p_i^3} \rfloor + \cdots$
那么，可以快速求解出 $C_n^k$ 的分解质因子，然后利用高精度乘法将它们相乘即可。

代码如下，

#include <iostream>
#include <vector>using namespace std;const int N = 5010;
int primes[N], cnt;
bool st[N];
int sum[N];void get_primes(int n) {//求n以内的质数for (int i = 2; i <= n; ++i) {if (!st[i]) {primes[cnt++] = i;}for (int j = 0; primes[j] <= n / i; ++j) {st[i * primes[j]] = true;if (i % primes[j] == 0) break;}}return;
}int get(int a, int p) {//求a!中质因子p的幂int res = 0;while (a) {res += a / p;a /= p;}return res;
}vector<int> mul(vector<int> a, int b) {vector<int> c;int t = 0;for (int i = 0; i < a.size() || t; ++i) {if (i < a.size()) {t = t + a[i] * b;}c.emplace_back(t % 10);t /= 10;        }while (c.size() > 1 && c.back() == 0) {c.pop_back();}return c;
}int main() {int a, b;cin >> a >> b;get_primes(a);for (int i = 0; i < cnt; ++i) {int p = primes[i];sum[i] = get(a, p) - get(b, p) - get(a - b, p);}vector<int> res = {1};for (int i = 0; i < cnt; ++i) {int p = primes[i];for (int j = 0; j < sum[i]; ++j) {res = mul(res, p);}}for (int i = res.size() - 1; i >= 0; --i) {cout << res[i];}cout << endl;return 0;
}