组合数的四种求解方式，递推，快速幂，Lucas，高精度&线性筛

本文介绍求组合数的四种方法

---

一、递推法——杨辉三角

1.1 问题描述

有q（q <= 10000）组询问，每组询问两个整数n、m（1 <= m <= n <= 2000），求C(n, m) mod (1e9 + 7)的值。

1.2 递推法

$$C_n^m=C_{n-1}^m+C_{n-1}^{m-1}$$

即第i个数选或不选

选：那么从剩下的n - 1个数中再挑出m - 1个数

不选：那么从剩下的n - 1个数中挑出m 个数

二者累加，即上述公式

时间复杂度O(n^2)

1.3 代码实现

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

int C[N][N];

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

二、快速幂

2.1 问题描述

有q（q <= 10000）组询问，每组询问两个整数n、m（1 <= m <= n <= 1e5），求C(n, m) mod (1e9 + 7)的值。

注意该问题的数据范围扩大到1e5了，再次使用递推法必然会MLE和TLE

2.2 快速幂

费马小定理

若p为质数，且a，p互质，则a ^ (p-1) ≡ 1(mod p)

而1e9 + 7是竞赛中非常常见的一个质数，那么给定 a，我们能够通过快速幂计算出其逆元为 a ^ (1e9 + 7 - 2) mod (1e9 + 7)

那么我们可以O(n log mod)预处理出阶乘fact[]和阶乘的逆元inv[]
$$\begin{array}{l}fact[i]=i\ast fact[i-1](mod1e9+7)\\ inv[i]=\frac{1}{i!}(mod1e9+7)\\ =(\frac{1}{(i-1)!}(mod1e9+7))\times (\frac{1}{i}(mod1e9+7))(mod1e9+7)\\ =inv[i-1]\times i^{1e9+5}(mod1e9+7)\end{array}$$
然后对于任意组合数查询，可以O(1)计算出C(n，m) = fact[n] * inv[m] * inv[n - m] % mod

预处理时间复杂度O(n logmod)

查询时间复杂度O(1)

2.3 代码实现

typedef long long LL;
const int N = 1e5 + 10, mod = 1e9 + 7;

LL fact[N], inv[N];

LL qp(LL a, LL b){
	LL res = 1;
	while(b){
		if(b & 1) res = res * a % mod;
		a = a * a % mod, b >>= 1;
	}
	return res;
}

void init(){
	fact[0] = inv[0] = 1;
	for(int i = 1; i < N; i ++){
		fact[i] = i * fact[i - 1] % mod;
		inv[i] = inv[i - 1] * qp(i, mod - 2) % mod;
	}
}

LL getC(LL m, LL n){
	return fact[n] * inv[n - m] % mod * inv[m] % mod;
}

三、卢卡斯定理

3.1 大组合数取模问题

给定正整数n，m，p的值，求出C(n，m) (mod p)的值。

其中 1 <= m <= n <= 1e18，1 <= p <= 1e5，保证 p 为质数

3.2 卢卡斯(Lucas)定理

3.2.1 定理内容

$$C_n^m\equiv C_{n/p}^{m/p}{C}_{nmodp}^{mmodp}(modp),其中p为质数$$

n mod p 和 m mod p都是小于 p 的数, p <= 1e5，这个我们可以直接快速幂求解

C(n / p, m / p)可以继续Lucas定理递归求解

递归出口：m = 0，返回1

3.2.2引理1

$$C_p^x\equiv 0(modp),0<x<p$$

证明：
$$\begin{array}{lc}{C}_p^x& =\frac{p!}{x!(p-x)!}\\ & =\frac{p}{x}\frac{(p-1)!}{(x-1)!(p-x)!}\\ C_p^x& \equiv p\cdot inv(x)C_{p-1}^{x-1}\equiv 0(modx)\end{array}$$

3.2.3 引理2

$$\begin{array}{l}(1+x{)}^p\equiv 1+x^p(modp)\end{array}$$

证明：
$$\begin{array}{l}(1+x{)}^p={\sum }_{i=0}^p{C}_p^i{x}^i\\ 由引理1，可知除了1，x^p,其它项modp为0\\ 故，(1+x{)}^p\equiv 1+x^p(modp)\end{array}$$

3.2.4 证明

$$
\begin{array}{l}
令n=ap+b,\ m=cp+d \
(1+x)^{n} \equiv \sum_{i=0}^{n}C_{n}{i}x^{i} \ (mod\ p)\cdots ① \
(1+x)^{n} \equiv (1+x)^{ap+b} \
\equiv ((1+x)^{p}){a}(1+x)^{b} \
\equiv (1+x^{p}){a}(1+x)^{b} \
\equiv \sum_{i=0}^{a}C_{a}{i}x^{ip}\cdot \sum_{j=0}^{b}C{b}{j}x^{j}(mod\ p)\cdots ②\
①中x^{{m}系数为C_{n}}{m} \
②中x^{m}=x{cp}\cdot x^{{d}的系数为C_{a}}{c}C_{b}^{d}\
故C_{n}^{m}\equiv C_{a}^{c}C_{b}{d}(mod\ p)\
故C_{n}^{m}\equiv C_{n/p}^{m/p}C_{n % p}^{m % p}(mod\ p)。证毕\
时间复杂度:O(p\log p + \log_{p}n)

\end{array}
$$

3.3 代码实现

typedef long long LL;
const int N = 1e5 + 10;
int n, m, p;
LL fact[N], inv[N];

LL qp(LL a, LL b){
	LL res = 1;
	while(b){
		if(b & 1) res = res * a % p;
		a = a * a % p, b >>= 1;
	}
	return res;
}

void init(){
	fact[0] = inv[0] = 1;
	for(int i = 1; i <= p; i ++){
		fact[i] = i * fact[i - 1] % p;
		inv[i] = inv[i - 1] * qp(i, p - 2) % p;
	}
}

LL getC(LL m, LL n){
    if(n < m) return 0;
	return fact[n] * inv[n - m] % p * inv[m] % p;
}

int lucas(int m, int n){
	if(!m) return 1;
	return lucas(m / p, n / p) * getC(m % p, n % p) % p;
}

3.4 OJ模板

P3807 【模板】卢卡斯定理/Lucas 定理 - 洛谷 | 计算机科学教育新生态 (luogu.com.cn)

注意getC对n < m特判

#include <iostream>
#include <cstring>
#include <algorithm>
using namespace std;
typedef long long LL;
const int N = 1e5 + 10;
int n, m, p;
LL fact[N], inv[N];

LL qp(LL a, LL b){
	LL res = 1;
	while(b){
		if(b & 1) res = res * a % p;
		a = a * a % p, b >>= 1;
	}
	return res;
}

void init(){
	fact[0] = inv[0] = 1;
	for(int i = 1; i <= p; i ++){
		fact[i] = i * fact[i - 1] % p;
		inv[i] = inv[i - 1] * qp(i, p - 2) % p;
	}
}

LL getC(LL m, LL n){
    if(n < m) return 0;
	return fact[n] * inv[n - m] % p * inv[m] % p;
}

int lucas(int m, int n){
	if(!m) return 1;
	return lucas(m / p, n / p) * getC(m % p, n % p) % p;
}

int main(){
	//freopen("in.txt", "r", stdin);
	int _ = 1;
	cin >> _;
	
	while(_--){
		cin >> n >> m >> p;
		init();
		cout << lucas(n, n + m);
		if(_)
			cout  << '\n';
	}
	return 0;
}

四、高精度、线性筛

4.1 大组合数不取模问题

给定两个整数n，m（1 <= m <= n <= 10000），求C（n, m）的值

组合数增长速度很快，不取模的情况下，需要用高精度存取。

C(100, 50)为30位数，C(1000, 500)是300位，C(10000，5000)是3009位数

4.2 算法流程

1. 筛质数，筛出2~n内的质数
2. 枚举质因子
   1. 计算C(n, m)内的质数个数s
   2. 利用高精度计算答案乘s次p

关于如何快速计算阶乘质因子p的次数可以做这道题：197. 阶乘分解 - AcWing题库

就是贡献法：对于质数p，p的倍数有多少个？p^2的倍数有多少个……我们直接累加Σ n / p^i，就是pc中的c，第一次，全体p的倍数贡献了1，第二次全体p^2倍数也累加了1，第三次……这样保证了不重不漏

时间复杂度的话由于涉及到调和级数的极限，但是1e4量级远达不到，我们以n = 10000，m = 5000的情况计算最坏时间复杂度

质数个数1229，Σs = 851，len = 3009，851 * 3009 大概也就2e6，效率还是很高的

4.3 代码实现

#include <iostream>
#include <cstring>
#include <algorithm>
using namespace std;

const int N = 10010;
int primes[N], cnt;
bool isprime[N];

void init(int n){
	for(int i = 2; i <= n; i++){
		if(!isprime[i]) primes[cnt++] = i;
		for(int j = 0; i * primes[j] <= n; j++){
			isprime[i * primes[j]] = 1;
			if(!(i % primes[j])) break;
		}
	}
}

int getnum(int n, int p){
	int res = 0;
	while(n) res += n / p, n /= p;
	return res;
}

int getps(int n, int m, int p){
	return getnum(n, p) - getnum(m, p) - getnum(n - m, p);
}

void mul(int *a, int b, int& len){
	int t = 0;
	for(int i = 0; i < len; i++){
		t += a[i] * b;
		a[i] = t % 10;
		t /= 10;
	}
	while(t)
		a[len] = t % 10, t /= 10, len++;
}

int main(){
	int n, m;
	cin >> n >> m;
	init(n);
	int A[N]{0}, len = 1;
	A[0] = 1;
	for(int i = 0; i < cnt; i++){
		int p = primes[i];
		int s = getps(n, m, p);
		while(s--)
			mul(A, p, len);
	}
	for(int i = len - 1; i >= 0; i--)
			cout << A[i];
	return 0;
}