【知识点9】组合数

1. 关于n!的一个问题

我们来讨论一个关于 n! 的一个问题——求n!中有多少个质因子p。

这个问题是什么意思呢？举个例子，6!=1x2x3x4x5x6，于是6!中有4个质因子2，因为2、4、6中各有1个2、2个2、1个2；而6!中有两个质因子3，因为3、6中均各有1个3。

对于这个问题，直观的想法是计算从1~n的每个数各有多少个质因子p，然后将结果累加，时间复杂度为 $O(nlogn)$，如下面的代码：

//计算n!中有多少个质因子p
int cal(int n,int p){
	int ans = 0;
	for(int i=2;i<=n;++i){	//遍历2~n 
		int temp = i;
		while(temp % p == 0){	//只要temp还是p的倍数 
			ans++;	//p的个数加1
			temp /= p;	//temp除以p 
		}
	}
	return ans; 
} 

但是这种做法对n很大的情况（例如n是$10^{18}$）是无法承受的，我们需要寻求速度更快的方法。

现在考虑10!中质因子2的个数，如图5-2所示：

显然10!中有因子$2^1$的数的个数为5，有因子$2^2$的数的个数为2，有因子$2^3$的数的个数为1，因此10!中质因子2的个数为 5 + 2 + 1 = 8。

仔细思考便可以发现此过程可以推广为：n!中有$(\frac{n}{p}+\frac{n}{p^2}+\frac{n}{p^3}+...)$个质因子p，其中除法均为向下取整。于是便得到了 $O(logn)$ 的算法，代码如下：

//计算n！中有多少个质因子p
int cal(int n,int p){
	int ans = 0;
	while(n){
		ans += n / p;	//累加n/p^k
		n /= p;			//相当于分母多乘一个p 
	}
	return ans; 
} 

//计算n！中有多少个质因子p
int cal(int n,int p){
	if(n<p){	//n<p时1~n中不可能有质因子p 
		return 0;	
	}else{
		return n / p + cal(n/p,p);		//返回n/p加上(n/p)！中质因子p的个数 
	}
} 

上面就是递推和递归两种实现求n!中有多少个质因子p。

有了这个算法，我们就能很快计算出n!的末尾有多少个零：

末尾0的个数就是指这个数总共有几个10因子，而10又能是2和5的乘积。
$n!$ 中质因子2的个数肯定大于质因子5的个数，所以质因子5的个数就是所要求的结果。

2. 组合数的计算

组合数$C_n^m$是指n个不同元素中选出m个元素的方案数(m ≤ n)，一般也可以写成$C(n,m)$，其定义式为 $C_n^m=\frac{n!}{m!(n-m)!}$，由三个整数的阶乘得到。通过定义可以知道，组合数满足$C_n^m=C_n^{n-m}$，且有$C_n^n=C_n^0=1$成立。本节讨论如下两个问题：

1. 如何计算$C_n^m$。
2. 如何计算$C_n^m\%p$。

2.1 如何计算$C_n^m$

《算法笔记》给出了三种计算$C_n^m$的方法：

方法一：通过定义式直接计算

由$C_n^m=\frac{n!}{m!(n-m)!}$可知，需要先计算n!，然后令其分别除以m!和(n-m)!即可。但显而易见的是，由于阶乘相当庞大，因此通过这种方式计算组合数能接受的数据范围会很小，即便使用long long类型来存储也只能承受n ≤ 20的数据范围。代码如下：

long long C(long long n,long long m){
	long long ans = 1;
	for(long long i = 1;i <= n;++i){
		ans *= i;
	}
	for(long long i = 1;i <= m;++i){
		ans /= i;
	}
	for(long long i = 1;i <= n - m;++i){
		ans /= i;
	}
	return ans;
} 

方法二：通过递推公式计算⭐⭐⭐⭐⭐

这种方式才是最常见的。

long long C(long long n,long long m){
	if(m == 0 || m == n){
		retunr 1;
	}else{
		return C(n-1,m) + C(n-1,m-1);
	}
}

在这种计算方法下完全不涉及阶乘运算，但是会产生另一个问题：重复计算。因此不妨记录下已经计算过的C(n,m)，这样当下次再次碰到时就可以作为结果直接返回了。如下面的递归代码：

long long res[67][67] = {0};
long long C(long long n,long long m){
	if(m == 0 || m == n){
		return 1;
	}
	if(res[n][m] != 0){
		return res[n][m];
	}
	return res[n][m] = C(n-1,m) + C(n-1,m-1);	//赋值给res[n][m]并返回 
}

或是下面这种把整张表都计算出来的递推代码：

const int n = 60;
void calC(){
	for(int i = 0;i <= n;++i){
		res[i][0] = res[i][i] = 1;	//初始化边界 
	}
	
	for(int i=2;i <= n;++i){
		for(int j = 1;j <= i / 2;j++){
			res[i][j] = res[i-1][j] + res[i-1][j-1];	//递推计算C(i,j)
			res[i][i-j] = res[i][j];	//C(i,i-j) = C(i,j) 
		}
	} 
} 

稍加画图可以发现，使用递归计算C(n,m)的时间复杂度和具体的数据有关，但单次计算C(n,m)不会超过$O(n^2)$，而递推计算所有C(n,m)的时间复杂度显然是$O(n^2)$，因此读者应当根据实际需要来选择使用递归还是递推。

方法三：通过定义式的变形来计算

这样，只要能保证每次除法都是整除，就能用这种“边乘边除”的方法避免连续乘法的溢出问题。那么，怎么证明每次除法都是整除呢？

事实上这等价于证明 $\frac{(n-m+1)\times (n-m+2)\times ...\times (n-m+i)}{1\times 2\times ...\times i}(1\le i\le m)$ 是个整数，不过这个结论显然成立，因为该式就是$C_{n-m+i}^i$的定义式展开的结果，而$C_{n-m+i}^i$显然是一个整数。

由此可以写出相应代码，时间复杂度为$O(m)$：

long long C(long long n,long long m){
	long long ans = 1;
	for(long long i = 1;i <= m;i++){
		ans = ans * (n - m + i) / i;	//注意：一定要先乘再除！！！ 
	}
	return ans; 
}

对比一下方法二和方法三在哪些情况下会溢出：

• 方法二在n=67、m=33时开始溢出；
• 方法三在n=62、m=31时开始溢出。

2.2 如何计算$C_n^m\%p$

下面会给出四种方法，它们有各自的适用的数据范围，需要依照具体情况选择使用，但是一般来说方法一已经能够满足需要。

方法一：通过递推公式计算⭐⭐⭐⭐⭐

这种方法基于第一个问题的方法二，也是最容易、最实用的一种。 只需要在原先的代码中适当的地方对$p$取模即可。为了说明问题方便，此处假设两倍的$p$不会超过int型。

在这种做法下，算法可以很好地支持$m\le n\le 1000$的情况，并且对$p$的大小和素性没有额外限制（例如$p\le 10^9$都是可以的）。代码如下：
递归：

int res[1010][1010] = {0};
int C(int n,int m,int p){
	if(m == 0||m == n){
		return 1;
	}
	if(res[n][m] != 0){
		return res[n][m];
	}
	return res[n][m] = (C(n-1,m,p)+C(n-1,m-1,p)) % p;
}

递推：

void calC(){
	for(int i=0;i <= n;i++){
		res[i][0] = res[i][i] = 1;		//初始化边界 
	}
	for(int i = 2;i <= n;i++){
		for(int j=1;j<=i/2;++j){
			res[i][j] = (res[i-1][j] + res[i-1][j-1]) % p;
			res[i][i-j] = res[i][j];
		}
	} 
}

参考文档

算法笔记