大组合数计算方法-CSDN博客

本文链接：https://blog.csdn.net/qq_39304630/article/details/77417017

组合数的公式：

$C_n^r=\frac{n*(n-1)*...*(n-r+1)}{r!}=\frac{n!}{r!(n-r)!}$

观察上述式子，读者应该发现了，我们主要还是采用第一个算式，因为乘法操作相对较少，速度比较快，也相对来说不容易溢出。

思路1：

暴力做乘法，并且使用组合数本身的性质进行优化，使用long long，保证在乘法过程中可以连乘更多的数。另外我们也可以优先做除法，再做乘法。当然这种优化方式能带来的支持运算的范围增大是相当小的。

A. $C_n^0=C_n^n=1$

B. $C_n^r=C_n^(n-r)$

#include<bits/stdc++.h>
using namespace std;
int main(){
	int n, r;
	long long mon = 1, son = 1; 
	cin >> n >> r;
	if( r == 0 || r == n) cout << 1;
	if( r > n - r) r = n - r;
	for( int i = 1; i <= r ; ++ i, -- n){
		son *= i;
		mon *= n;
	}
	cout << mon/son ;
	return 0;
}

思路2：

其实我们就是要想办法让数在连乘过程中不要因为达到上界发生溢出而导致错误，所以我们可以考虑在做乘法的过程中同时除。这样会导致当前没法整除，而在计算机中，整数中的整除是下取整的，所以用整数直接做除法是不正确的。解决办法有两个：

A. 我们将除数做缓存，每次可以整除的时候，拿出来除，缩小被除数，并从缓存表中删去这个除数

B. 直接用浮点数存，虽然浮点数会有精度损失，但是损失的量是十分小的，而组合数的结果一定是整数，因此我们将浮点结果上取整就好啦~~ ( 但是double类型的存储范围要比long long广，不过需要同时考虑到浮点数本身带来的精度损失，因此我们只当作一种思维拓展，而不把他们视为正规解法 )。

C.我们如果对上面的公式两边取对数，就可以得到

$\bg_white ln(C_n^r)=ln(n!)-ln(r!)-ln((n-r)!)=\sum_{i=1}^{n}ln(i)-\sum_{i=1}^{r}ln(i)-\sum_{i=1}^{n-r}ln(i)$

使用对数的性质可以帮助我们压缩存储范围，同时改乘除法为加减法，但问题是计算机本身的存储浮点数带来的精度损失问题。

思路3：

套大整数的模板。但是大整数采用数组实现，所以速度是比较慢的，另外，有时会碰到取模的要求，在大整数类型中做除法和取模操作是十分不方便的。

思路4：

采用杨辉三角递推式。对于组合数，有这么一个性质：

$\bg_white C_n^r=C_{n-1}^r+C_{n-1}^{r-1}$

我们以取物问题为例，理解上述公式： $C_n^r$ 表示从 $n$ 个物体中拿出 $r$ 个，那么如果我们针对第n个物体考虑，可以将这个原问题划分为两个小的子问题：

1. 我们将第n个物体当作我们拿出的第r个物体，此时我们需要在前n-1个物体中拿出r-1个，即 $C_{n-1}^{r}$ ；

2. 我们保证不拿第n个物体，此时我们需要在前n-1个物体中拿出r个，即 $C_{n-1}^{r-1}$ ；

#include<bits/stdc++.h>
using namespace std;

int C[1010][1010];

int main(){
	int n, r, i, j;
	cin >> n >> r;
	C[1][0] = C[1][1] = 1;
	for( i = 2; i <= n ; ++ i){
		C[i][0] = C[i][i] = 1;
		for( j = 1; j < n ; ++ j ){
			C[i][j] = C[i-1][j] + C[i-1][j-1];
		}
	}
	cout << C[n][r];
	return 0;
}

相信读者也发现了，递推式本身过于消耗空间，因此递推法能解决的范围也很有限。当然细心的读者应该发现了，这边的递推数组是一个下三角矩阵，因此可以采用坐标转换的方式化为一维存储，可以省去一半的空间。

注意，这个方法需要在打表计算的同时就进行取模操作 $(a+b)%p=(a%p+b%p)%p$ ，因为在不取模进行累加时，数据的增长是异常迅速的，因此很快就会数据溢出报错。

思路5：

用数论的知识，所有的数都可以使用素因数进行分解：即对任意一个数n，存在一个素数序列 $P_i$ ，使得

$n=\prod_{i=1}^{k}P_i^{Q_i}$

我们使用质因数分解的办法，将除数和被除数分解成小的素数因子后，再统计计算结果。因而我们需要一个素数表，我们可以使用欧拉筛快速获得这张表。而也正是因为这样，我们的算法也有了局限性，由于内存的限制，我们没办法获得很大的素数，也就导致了当需要被分解的数中存在大素数因子的时候，我们就束手无测。（或者当素数表使用完之后，继续向后枚举素数，这样会降低算法性能），不过一般情况下，至少处理百万级的因数时，是完全ok的。

#include<bits/stdc++.h>
using namespace std;
const int MaxN = 1000010;

bool isPrime[MaxN];
int Primes[MaxN], cnt, Factor[MaxN];

void EulrSift(){ // 欧拉筛 
	int i, j;
	for( i = 2; i <= MaxN; ++ i ){
		if( !isPrime[i] ) {
			Primes[cnt++] = i;
			for( j = 0;  j < cnt && Primes[j] * i < MaxN ; ++ j){
				isPrime[ Primes[j] * i ] = true; 
				if( i % Primes[j]) break;
			}
		} 
	}
}

void PrimeResolve(int num, int delta){ // 质因数分解 
	for( int i = 0 ; num != 1 && i < cnt ; ++ i ){
		while( num % Primes[i] == 0 && num != 1){
			num /= Primes[i];
			Factor[i] += delta;
		}
	}
}

int quickMul(int a,int b){ // 快速幂 
	int r = 1;
	while( b ){
		if( b & 1 ){
			r *= a;
		}
		b >>= 1;
		a *= a;
	}
	return r;
}

int main(){
	int n, r, m, i, ans = 1;
	cin >> n >> r;
	m = n;
	EulrSift();
	for( i = 1; i <= r ; ++ i, -- n){
		PrimeResolve(n, 1);
		PrimeResolve(i, -1);
	}
	for( i = 0; i < cnt &&  Primes[i] <= m  ; ++ i){
		if ( Factor[i] != 0 ){
			ans *= quickMul(Primes[i], Factor[i]);
		}
	}
	cout << ans;
	return 0;
}

思路6：

从这里开始，我们会接触到求大组合数+取模的问题。观察组合数的运算公式，我们发现运算过程中有除法，但是在取模条件下，我们可以把除法使用逆元的方法转化为乘法来做。费马定理的要求过于苛刻，我们可以采用扩展欧几里得来求。欧拉定理比费马小定理的通用性更强，但是有时候模数过大，可能导致欧拉定理无法使用，另外，对于欧拉定理而言，一般我们不使用他的关键原因在于，本身欧拉定理告诉我们，如果 $gcd(A,P) = 1$ ，那么一定会有：

$A^{\varphi(P)} =1(mod \ P)$

但是即便我们使用欧拉筛取求解这个 $\varphi(P)$ ，也需要 $\dpi{100} \bg_white \dpi{100} \bg_white O(sqrt(P))$ 的复杂度，并不值得。但是我们还是复习一遍欧拉函数吧(233)。

我们针对问题： $C_n^r % p$ ，给出下面两种解决办法：

扩展欧几里得版：

#include<bits/stdc++.h>
using namespace std;

// A*x0 + B*y0 = B*x + (A-A/B*B) *y
// A*x0 + B*y0 = A*y + B * (x-A/B*y) 

int extgcd(int A,int B,int &x,int &y){
	if( B == 0 ){
		x = 1;
		y = 0;
		return A; 
	}
	int gcd = extgcd(B, A%B, x, y);
	int x0 = y;
	int y0 = x - A / B * y;
	x = x0;
	y = y0;
	return gcd;
}

// A * x = 1 ( mod P ) = > A * x + P * y = 1
int getInv(int A,int p){
	int x, y, gcd = extgcd(A,p,x,y);
	if( gcd != 1)return -1; // 扩展欧几里得无法求出逆元 
	p /= gcd;
	return ( x % p + p ) % p; // 取最小正整数解 
}

int main(){
	int i, n, r, p, ans = 1, inv;
	cin >> n >> r >> p;
	for( i = 1; i <= r; ++ i, -- n){
		ans = ans * n % p;
		if( ans % r == 0 ) ans /= r; //如果可以除掉 就不要求逆元 
		else if( r >= 2 ){ // 1不需要求逆元 
			inv = getInv(r,p);
			if( inv == -1 ){
				cout << "扩展欧几里得无法求出逆元,该程序无法求解"<<endl;
				return  0;
			}
			ans = ans * inv % p;
		}
	}
	cout << ans;
	return 0;
}

欧拉函数版：

#include<bits/stdc++.h>
using namespace std;

int cnt, phi[100010], Primes[100010];

int Phi(int n){ // 质因数分解 结合 欧拉函数性质  
	int res = n;
	for( int i = 2; i * i <= n; ++ i){ 
		if( n % i == 0 ){
			res = res / i * ( i - 1 );
		}
		while( n % i == 0 ){
			n /= i;
		}
	}
	if( n > 1 ){ // 除干净了因子后 n = 1 或者是最后一个质因数的幂次 
		res = res / n * ( n - 1 );
	}
	return res;
} 

int quickMul(int A, int B, int P){
	int r = 1;
	while( B ){
		if( B & 1 ){
			r = r * A % P;
		}
		B >>= 1;
		A = A * A % P;
	}
	return r;
}

int gcd(int A,int B){
	while(B^=A^=B^=A%=B);
	return A;
}

int main(){
	int n, r, p, i, ans = 1, phi, inv; // son用来存储那些没有逆元的因子 
	cin >> n >> r >> p;
	phi = Phi(p);
	for( i = 1 ; i <= r ; ++ i, -- n ){
		ans = ans * n ;
		if( gcd(i,p) != 1 ){
			cout << " 欧拉定理无法求逆元，此程序无法求解。 " << endl; 
			return 0;
		}else{
			inv = quickMul(i, phi, p);
			ans = ans * inv % p;
		}
		ans %= p;
	}
	cout << ans ;
}

有些读者可能会想到，本身求逆元是需要时间的，能否在能整除时做除，不能做整除的时候求逆元改为乘法，同时每一步保持取模，答案是否定的，我们不能在乘法、除法、逆元分别取模后再合并。

以 $\frac{12}{3}(mod\ 9)=4$ 为例，如果我们做完分子的取模后，再除掉分母后再取模，那么答案为：

$\frac{12(mod\ 9)}{3}(mod\ 9)=\frac{3}{3}(mod\ 9)=1$

原因在于，取模本身不满足和除法得交换律：

1.因此，我们只能选择，在可以整除的时候，整除，但是一旦除法参与了，那么取模运算必须放到最后进行。

2.全部改用逆元做乘法。

我们使用这些算法的目的就是为了避免乘法溢出，因此方案1肯定是不切实际的啦~

思路7：卢卡斯定理

卢卡斯定理提供了一种在模数为素数的情况下，由于n和r远大于p，所以在除法过程中求逆元的时候很容易遇到与p不满足互质的情况，例如2*p，3*p，从而由于不满足互质导致无法求解逆元，但是卢卡斯定理提供了一种这种情形下还能求解大组合数的方案。

定理内容：对于非负整数 $\bg_white n,r,p$ (其中p为素数)，那么有如下的同余方程：

$\dpi{120} \bg_white C_n^r=\prod_{i=0}^{k}C_{n_i}^{r_i}(mod\ p)$

其中， $n_i,r_i$ 需要满足如下等式：

$n_k*p^k+n_{k-1}*p^{k-1}+......+n_1*p+n_0$

$r_k*p^k+r_{k-1}*p^{k-1}+......+r_1*p+r_0$

不难发现，这就是分别将 $n,r$ 做 $p$ 进制下的分解罢了。

在实际的算法实现上，我们会这么使用这个定理，与上述是等价的：

$Lucas_n^r=C_{n%p}^{r%p}*Lucas_{n/p}^{r/p}$

注意上述公式中的除法，是整除，其实上述公式完全等价于卢卡斯定理原本的公式。

需要注意的是，在这个定理中，如果发生了 $n<r$ 的情况，那么规定 $C_n^r=0$ 。

公式证明：求解 $C_p^x(mod\ p)$ ，设 $x<p,x!=0$

$C_p^x=\frac{p!}{x!(p-x)!}=\frac{p*(p-1)!}{x*(x-1)!*(p-x)!}=\frac{p}{x}*C_{p-1}^{x-1}$

由于我们已经设 $x<p$ 且 $p$ 为素数，因此我们知道 $x$ 肯定有关于 $p$ 的逆元，所以上述式子可以转换为：

$\bg_white \dpi{120} \bg_white C_p^x(mod\ p)=\frac{p}{x}*C_{p-1}^{x-1}(mod\ p)=p*inv(x\ |\ p)*C_{n-1}^{r-1}(mod\ p)$

$inv(x|p)$ 表示x在模p下的逆元。

此时等式右端是 $p$ 的倍数，因此 $C_p^x=0(mod\ p)$ 。

根据二项式定理：

$\bg_white (1+x)^n=\sum_{i=0}^{n}C_n^i*x^i$

我们知道：

$\bg_white (1+x)^p=\sum_{i=0}^{p}C_p^i*x^i$

根据我们上面所证，这个上述式子的展开式中，只有 $i=0$ 和 $i=p$ 的两项模 $p$ 不为0。故我们得到：

$(1+x)^p=1+x^p(mod\ p)$

我们先预设：

$\bg_white \left\{ \begin{aligned} n/p = n_k \\ n\ mod\ p = res_k \\ \end{aligned} \right.$

我们对 $(1+x)^n$ 做p进制分解分解：

$(1+x)^n\\ =(1+x)^{n_k*p+res_k}\\ =(1+x)^{n_k*p}*(1+x)^{res_k}\\ =[(1+x)^p]^{n_k}*(1+x)^{res_k}...(1)\\ =(1+x^p)^{n_k}*(1+x)^{res_k}...(2)\\ =\sum_{i=0}^{n_k}C_{n_k}^i*(x^{p})^{i}*\sum_{j=0}^{res_k}C_{res_k}^{j}*x^j...(3)\\ =\sum_{i=0}^{n_k}\sum_{j=0}^{res_k}C_{n_k}^i*C_{res_k}^{j}*x^{i*p+j}(mod\ p)...(4)$

(1)->(2)：在%P的条件下，我们用了上述推出的结论- $(1+x)^p=1+x^p(mod\ p)$

(2)->(3)：使用两次二项式展开定理

(3)->(4)：多项式的乘法就是分别乘法后求和

注意(*)式中我们用了上面推导的结论： $(1+x)^p=1+x^p(mod\ p)$

所以我们得到了：

$(1+x)^n=\sum_{r=0}^{n}C_n^r*x^r=\sum_{i=0}^{n_k}C_{n_k}^ix^{i*p}*\sum_{j=0}^{res_k}C_{res_k}^jx^j\\=\sum_{i=0}^{n_k}\sum_{j=0}^{res_k}C_{n_k}^i*C_{res_k}^{j}*x^{i*p+j}(mod\ p),(x<p,res_k<p)$

由于上述式是恒等式，因此比较两端 $x^k$ 的系数可知：

$C_n^r*x^r=C_{n_k}^i*C_{res_k}^{j}*x^{i*p+j}(mod\ p),(x<p,res_k<p)$

即：

$C_n^r=C_{n_k}^i*C_{res_k}^{j}(mod\ p),(x<p,res_k<p)$

$n=n_k*p+res_k$

$r=i*p+j$

则

$C_n^r=C_{n_k}^i*C_{res_k}^j=C_{n/p}^{r/p}*C_{n%p}^{r%p}$

证毕。

由于 $r%p<p,n%p<p$ ，因此对于 $C_{n%p}^{r%p}$ 我们就可以比较容易地计算（因为p是质数，所以对与x<p,一定有gcd(x,p)=1，互质条件下可以求逆元）

在卢卡斯定理的加持下，我们可以保证在如下式子中：

$\dpi{120} \bg_white C_n^r=\prod_{i=0}^{k}C_{n_i}^{r_i}(mod\ p)$

一定会有 $n_i,r_i<p$ (按照p进制数理解，每一位上的权值不可能大于它所表示的数的进制值)。

而又由于 $p$ 为素数，因此我们又可以用逆元啦，完美！

#include<bits/stdc++.h>
using namespace std;

int quickMul(int A,int B,int P){
	int res = 1;
	while( B ){
		if( B & 1 ){
			res = res * A % P;
		}
		B >>= 1;
		A = A * A % P;
	}
	return res;
}

int inv(int n,int p){ // 费马小定理求解逆元  n^(p-1)=1(mod p) 
	return quickMul(n, p-2, p);
}

int C(int n,int r,int p){
	int res = 1;
	for( int i = 1; i <= r; ++ i, -- n){
		res =  res * n % p * inv(i, p) % p;
	}
	return res;
}

int Lucas(int n,int r,int p){
	return r == 0 ? 1 : Lucas(n/p, r/p, p) * C (n%p, r%p, p) %  p; 
}

int main(){
	int n, r, p;
	cin >> n >> r >> p;
	if( n < r || n == 0 ) {
		cout << -1 ; // 视具体题目要求而定，毕竟数学上规定 0!=1 
		return 0;
	}
	if( p == 1 ){
		cout << 0; 
		return 0;
	}
	cout << Lucas(n, r, p);
	return 0;
}

来看一道模板题吧~

洛谷 P3807 【模板】卢卡斯定理

提交16.49k

通过7.28k

时间限制1.00s

内存限制125.00MB

题目背景

这是一道模板题。

题目描述

给定整数 n, m, pn,m,p 的值，求出 C_{n + m}^n \bmod pCn+mnmodp 的值。

输入数据保证 pp 为质数。

注: CC 表示组合数。

输入格式

本题有多组数据。

第一行一个整数 T，表示数据组数。

对于每组数据:

一行，三个整数 n, m, p。

输出格式

对于每组数据，输出一行，一个整数，表示所求的值。

输入输出样例

输入 #1复制

2
1 2 5
2 1 5

输出 #1复制

3
3

#include<bits/stdc++.h>
using namespace std;

typedef long long ll;

ll quickMul(ll A,ll B,ll P){
	ll res = 1;
	while( B ){
		if( B & 1 ){
			res = res * A % P;
		}
		B >>= 1;
		A = A * A % P;
	}
	return res;
}

ll inv(ll n,ll p){ // 费马小定理求解逆元  n^(p-1)=1(mod p) 
	return quickMul(n, p-2, p);
}

ll C(ll n,ll r,ll p){
	ll res = 1;
	for( int i = 1; i <= r; ++ i, -- n){
		res =  res * n % p * inv(i, p) % p;
	}
	return res;
}

ll Lucas(ll n,ll r,ll p){
	return r == 0 ? 1 : Lucas(n/p, r/p, p) * C (n%p, r%p, p) %  p; 
}

int main(){
	ll T, n, m, p;
	cin >> T;
	while( T -- ){
		cin >> n >> m >> p;
		n += m;
		if( p == 1 ) cout << 0 << endl;
		else cout << Lucas(n, m, p) << endl;
	}
	return 0;
}

奇技淫巧：

在使用卢卡斯定理的时候，我们还可以预处理阶乘，这样可以加快速度，在预处理阶乘时，我们采用公式二：

$\dpi{100} \bg_white C_n^r=\frac{n!}{r!(n-r)!}$

事实上我们可以得到：

$\bg_white \dpi{100} \bg_white C_n^r(mod\ p)=\frac{n!}{r!(n-r)!}(mod\ p)\\ =n!(mod\ p)*inv(r!\ |\ p)(mod\ p)*inv(n-r!\ |\ p)(mod\ p)\\ =n!(mod\ p)*inv(r!(mod\ p)\ |\ p)(mod\ p)*inv(n-r!(mod\ p)\ |\ p)(mod\ p)$

即我们需要证明： $\frac{x!}{y!}(mod\ p)=x!*inv(y!(mod\ p)\ |\ p)(mod\ p)$

设 $z=y!(mod\ p)$ ，即 $y!=z+k*p$ ，那么有 $z<p$ ，且 $p$ 为素数，故 $z$ 必有逆元，则有：

$z*inv(z\ |\ p)=1(mod\ p)$

那么：

$y!*inv(z)=(z+k*p)*inv(z\ |\ p)=z*inv(z)+k*p*inv(z\ |\ p)\\ =z*inv(z)(mod\ p)=1$

证毕。

#include<bits/stdc++.h>
using namespace std;
const int MaxN = 1e5+10;

typedef long long ll;
ll Factorial[MaxN << 1];

ll quickMul(ll A,ll B,ll P){
	ll res = 1;
	while( B ){
		if( B & 1 ){
			res = res * A % P;
		}
		B >>= 1;
		A = A * A % P;
	}
	return res;
}

ll inv(ll n,ll p){ // 费马小定理求解逆元  n^(p-1)=1(mod p) 
	return quickMul(n, p-2, p);
}

ll C(ll n,ll r,ll p){ // 这里需要注意 n < r 时 是非法状态 返回 0 
	return n < r ? 0 : Factorial[n] * inv(Factorial[r],p) % p * inv(Factorial[n-r],p) % p; 
}

ll Lucas(ll n,ll r,ll p){
	return r == 0 ? 1 : Lucas(n/p, r/p, p) * C (n%p, r%p, p) %  p; 
}

int main(){
	ll T, n, m, p, i;
	cin >> T;
	Factorial[1] = 1;
	while( T -- ){
		cin >> n >> m >> p;
		n += m;
		memset(Factorial, 0, sizeof(Factorial));
		Factorial[1] = 1;
		if( p == 1 ){
			cout << 0 << endl;
		} 
		else {
			m = m > (n >> 1) ? n - m : m;
			for ( i = 2; i <= n; ++ i){
				Factorial[i] = Factorial[i-1] * i % p;
			}
			cout << Lucas(n, m, p) << endl;
		}
	}
	return 0;
}