计算组合数的几种方法总结

前言

组合数就是 Cmn C n m C_n^m，是排列组合中非常重要的一部分，最近做到了几道与计算组合数有关的题，想在此总结一下如何编程计算组合数。
某位大佬在我之前也写过类似的文章，这是他的博客：https://blog.csdn.net/u010582475/article/details/47707739

递推（杨辉三角）

先给出递推公式： Cmn=Cmn−1+Cm−1n−1 C n m = C n − 1 m + C n − 1 m − 1 C_n^m=C_{n-1}^m+C_{n-1}^{m-1}
证明：从n个数中选m个数，若第n个数不选，则相当于从n-1个数中选m个数，即 Cmn−1 C n − 1 m C_{n-1}^m；若第n个数选，则相当于从n-1个数中选m-1个数，即 Cm−1n−1 C n − 1 m − 1 C_{n-1}^{m-1}。那么总方案数即为两者之和。
当然，也可以根据组合数公式，展开来证明其相等，这里就不写出来了，有兴趣的话也可以自己去推。
如果写成数组形式，就是c[n][m]=c[n-1][m]+c[n-1][m-1]，这同时是杨辉三角的递推公式，所以组合数满足杨辉三角。边界是c[0][0]=1。若只要求一个组合数的话，也可以滚动一下这个数组，节约空间，但同时m这一维要倒着枚举，就像01背包一样。程序如下：

#include<iostream>
using namespace std;
int n,m;
long long c[10005];
int main()
{cin>>n>>m;m=min(m,n-m);c[0]=1;for (int i=1;i<=n;i++){for (int j=m;j>=1;j--){c[j]=c[j]+c[j-1];}}cout<<c[m];
}

因为 Cmn=Cn−mn C n m = C n n − m C_n^m=C_n^{n-m}，所以m可以取m和n-m中小的那一个，以节省时间。但复杂度还是过高，约为O( n2 n 2 n^2)。

乘法逆元

先强调一遍，这种方法只适用于对答案模一个大质数的情况。
一般组合数会非常非常大，可能long long都会爆，所以题目会让我们模一个大质数，然后输出。既然要取模，那就要用到数论中的内容了。
观察组合数的公式： Cmn=n!m!(n−m)!=n∗(n−1)∗...∗(n−m+1)m! C n m = n ! m ! ( n − m ) ! = n ∗ ( n − 1 ) ∗ . . . ∗ ( n − m + 1 ) m ! C_n^m=\frac{n!}{m!(n-m)!}=\frac{n*(n-1)*...*(n-m+1)}{m!}
对于取模来说，乘法可以每做一步都取模，但除法却不行。为了解决除法的取模问题，人们发明了乘法逆元这种奇妙的东西。若a*b=1(mod p)，则称b为a在模p下的乘法逆元，一般认为， b<p b < p b

，b记为 a−1 a − 1 a^{-1}。当然，a和p要互质，不然这个乘法逆元不存在。但乘法逆元有什么用呢？看如下的式子（以下等式均在mod p的前提下）：
ab a b \frac{a}{b}= ab∗b∗b−1 a b ∗ b ∗ b − 1 \frac{a}{b}*b*b^{-1}=a* b−1 b − 1 b^{-1}（先是乘以1，再是约掉b）
这样我们就把除法取模转化成了乘以它的乘法逆元再取模。
那么乘法逆元如何计算呢？

递推

求i在mod p下的乘法逆元（gcd(i,p)=1, i<p i < p i

）：
令p=k*i+r，则k=p/i，r=p%i
显然，k*i+r=0(mod p)
两边同乘 i−1∗r−1 i − 1 ∗ r − 1 i^{-1}*r^{-1}，得到k* r−1+i−1 r − 1 + i − 1 r^{-1}+i^{-1}=0(mod p)
移项，再进行处理，得 i−1=−k∗r−1=(p−k)∗r−1 i − 1 = − k ∗ r − 1 = ( p − k ) ∗ r − 1 i^{-1}=-k*r^{-1}=(p-k)*r^{-1}(mod p)
把k和r用p和i代掉，得 i−1=(p−p/i)∗(p%i)−1 i − 1 = ( p − p / i ) ∗ ( p % i ) − 1 i^{-1}=(p-p/i)*(p\% i)^{-1}
若用inv[i]表示i的乘法逆元，则上式可以写成inv[i]=(p-p/i)*inv[p%i]
复杂度为线性。伪代码如下：

long long inv[10000005];
inv[1]=1;
long long ni(int x,int p)
{if (inv[x]!=0) return inv[x];inv[x]=(p-p/x)*ni(p%x,p)%p;return inv[x];
}

费马小定理

费马小定理： ap−1 a p − 1 a^{p-1}=1(mod p)（p是质数，gcd(a,p)=1)
求a在mod p下的乘法逆元可以求 ap−2 a p − 2 a^{p-2}，用快速幂就好了。

long long pow(long long x,long long y,long long p)//快速幂，求x^y%p
{if (y==0) return 1;long long z=pow(x,y/2)%p;if (y%2==0) return z*z%p;return z*z%p*x%p;
}
long long ni(long long x,long long p)//求x在mod p下的乘法逆元
{return pow(x,p-2,p);
}

欧拉定理

费马小定理是欧拉定理在p为质数下的特殊情况，既然费马小定理可以，那么欧拉定理也可以。但使用乘法逆元求组合数的前提是模一个大质数，能用简单的费马小定理，为什么要用麻烦的欧拉定理呢？所以我在此只是提一提，并不想多讲。

刚刚介绍了三种求乘法逆元的方法，现在回到主题，如何求组合数。
还记得刚才的公式吗？
Cmn=n!m!(n−m)!=n∗(n−1)∗...∗(n−m+1)m! C n m = n ! m ! ( n − m ) ! = n ∗ ( n − 1 ) ∗ . . . ∗ ( n − m + 1 ) m ! C_n^m=\frac{n!}{m!(n-m)!}=\frac{n*(n-1)*...*(n-m+1)}{m!}
对于分子，我们可以边乘边mod p，对于分母，只需要求出分母的乘法逆元，再乘以分子就行了。
以下是程序（此程序用的是费马小定理求乘法逆元）：

#include<iostream>
using namespace std;
int n,m,p;
long long pow(long long x,long long y)
{if (y==0) return 1;long long z=pow(x,y/2)%p;if (y%2==0) return z*z%p;return z*z%p*x%p;
}
long long ni(long long x,long long p)
{return pow(x,p-2);
}
long long c(int n,int m,int p)
{long long x=1,y=1;for (int i=n;i>=n-m+1;i--) x=x*i%p;for (int i=1;i<=m;i++) y=y*i%p;return x*ni(y,p)%p;
}
int main()
{cin>>n>>m>>p;cout<<c(n,m,p);
}

再说一遍，p一定要是大质数，起码得大于m。为什么呢？若要求a在模p下的乘法逆元，必须要保证a与p互质。在此题中，要求m!的乘法逆元，那么p必须与m!互质，那么p就要大于m，且是个质数。

卢卡斯定理

先把这种方法的适用条件写在前面，适用于对答案模一个质数的情况。
和上面的乘法逆元求组合数的条件对比一下，只相差了一个大字。也就是说这个质数不用很大，非常小也行。
介绍一下卢卡斯定理，公式： Cmn=Cm/pn/p∗Cm%pn%p C n m = C n / p m / p ∗ C n % p m % p C_n^m=C_{n/p}^{m/p}*C_{n\%p}^{m\%p}（mod p），要求p为质数。
公式很好记吧，至于证明呢，看懂了也下次也不会证，会证了也没什么用（主要是我不会），我就不证明了。
观察公式，有没有发现，在n和m都小于p的时候，公式一点用都没有（ C00=1 C 0 0 = 1 C_0^0=1）。所以这个公式是在n>=p或m>=p的情况下使用的。这样可以减小n和m，使之小于p，再用乘法逆元去求组合数。
其实卢卡斯定理就是一个小优化，大部分和乘法逆元求组合数一样，不过要注意在n=0和m=0时的边界处理。

#include<iostream>
using namespace std;
int n,m,p;
long long pow(long long x,long long y)
{if (y==0) return 1;long long z=pow(x,y/2)%p;if (y%2==0) return z*z%p;return z*z%p*x%p;
}
long long ni(long long x)
{return pow(x,p-2);
}
long long c(int n,int m)
{if (m==0) return 1%p;if (n==0) return 0;if (n>=p||m>=p) return c(n/p,m/p)*c(n%p,m%p)%p;//卢卡斯定理核心语句long long x=1,y=1;for (int i=n;i>=n-m+1;i--) x=x*i%p;for (int i=1;i<=m;i++) y=y*i%p;return x*ni(y)%p;
}
int main()
{cin>>n>>m>>p;cout<<c(n,m);
}

质因数分解

乘法逆元只能处理模数为大质数的情况，卢卡斯定理只能处理模数为质数的情况，那有没有一种方法能处理模数不是质数的情况呢？显然是有的。而且不取模也是可以的。
我们可以把组合数中要乘或除的每一个数分解质因数，再把分母的质因数减掉，最后把剩下的质因数乘起来，边乘边模p就行了。
那如何分解质因数呢？可以用欧拉筛把每个数的最小质因数求出来，把i的最小质因数的编号保存在min_prime[i]里。
具体看代码吧。

#include<iostream>
using namespace std;
int n,m,p,b[10000005],prime[1000005],t,min_prime[10000005];
void euler_shai(int n)//用欧拉筛求出1~n中每个数的最小质因数的编号是多少，保存在min_prime中
{for (int i=2;i<=n;i++){if (b[i]==0) {prime[++t]=i;min_prime[i]=t;}for (int j=1;j<=t&&i*prime[j]<=n;j++){b[prime[j]*i]=1;min_prime[prime[j]*i]=j;if (i%prime[j]==0) break;}}
}
long long c(int n,int m,int p)//计算C(n,m)%p的值
{euler_shai(n);int a[t+5];//t代表1~n中质数的个数 ，a[i]代表编号为i的质数在答案中出现的次数 for (int i=1;i<=t;i++) a[i]=0;//注意清0，一开始是随机数 for (int i=n;i>=n-m+1;i--)//处理分子 {int x=i;while (x!=1){a[min_prime[x]]++;//注意min_prime中保存的是这个数的最小质因数的编号（1~t） x/=prime[min_prime[x]];}}for (int i=1;i<=m;i++)//处理分母 {int x=i;while (x!=1){a[min_prime[x]]--;x/=prime[min_prime[x]];}}long long ans=1;for (int i=1;i<=t;i++)//枚举质数的编号，看它出现了几次{while (a[i]>0){ans=ans*prime[i]%p;a[i]--;}}return ans;
}
int main()
{cin>>n>>m>>p; m=min(m,n-m);//小优化 cout<<c(n,m,p);
}

总结

写了这么多种算法，每种算法都有其优点与局限性。递推写起来快，思维简单，但时间复杂度高。乘法逆元用得比较普遍，因为一般都是模一个大质数，复杂度也几乎是线性的。卢卡斯定理只会在特定的题目里做到，但其实编程复杂度并不高，就是在乘法逆元的基础上加几句话。质因数分解的适用性最广，编程复杂度也最高，这就是完美的代价吧。