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

1 基础知识

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

1. 数a和数b的乘积模上p，等于数a模上p和数b模上p的乘积。即，
   $$(a\cdot b)modp=(amodp)\cdot (bmodp)$$
2. 数a除以数b的结果模上p，并不等于数a模上p除以数b模上p。即，
   $$(a/b)modp\ne (amodp)/(bmodp)$$

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

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

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

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}modp=？$$
其中$n$和$k$的数据范围在$10^{18}$之内，而$p$是质数，且它的范围在$10^6$以内。

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

Lucas定理：当模数p是质数时，有以下等式成立，
$$C_n^{k}modp=C_{nmodp}^{kmodp}\cdot C_{n/p}^{k/p}modp$$
其中$kmodp$和$nmodp$是$p$以内的数，可直接计算组合数$C_{nmodp}^{kmodp}$；而对于$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\le i\le 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;
}

2 模板

暂无。。。

3 工程化

暂无。。。