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

最新推荐文章于 2022-01-02 16:01:33 发布

liuzibujian

最新推荐文章于 2022-01-02 16:01:33 发布

阅读量1.9w

点赞数 35

本文链接：https://blog.csdn.net/liuzibujian/article/details/81346595

版权

前言

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

递推（杨辉三角）

先给出递推公式： $C_n^m=C_{n-1}^m+C_{n-1}^{m-1}$
证明：从n个数中选m个数，若第n个数不选，则相当于从n-1个数中选m个数，即 $C_{n-1}^m$ ；若第n个数选，则相当于从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];
}

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

乘法逆元

先强调一遍，这种方法只适用于对答案模一个大质数的情况。
一般组合数会非常非常大，可能long long都会爆，所以题目会让我们模一个大质数，然后输出。既然要取模，那就要用到数论中的内容了。
观察组合数的公式： $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记为 $a^{-1}$ 。当然，a和p要互质，不然这个乘法逆元不存在。但乘法逆元有什么用呢？看如下的式子（以下等式均在mod p的前提下）：
$\frac{a}{b}$ = $\frac{a}{b}*b*b^{-1}$ =a* $b^{-1}$ （先是乘以1，再是约掉b）
这样我们就把除法取模转化成了乘以它的乘法逆元再取模。
那么乘法逆元如何计算呢？

递推

求i在mod p下的乘法逆元（gcd(i,p)=1, $i<p$ ）：
令p=k*i+r，则k=p/i，r=p%i
显然，k*i+r=0(mod p)
两边同乘 $i^{-1}*r^{-1}$ ，得到k* $r^{-1}+i^{-1}$ =0(mod p)
移项，再进行处理，得 $i^{-1}=-k*r^{-1}=(p-k)*r^{-1}$ (mod p)
把k和r用p和i代掉，得 $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];
}

费马小定理

费马小定理： $a^{p-1}$ =1(mod p)（p是质数，gcd(a,p)=1)
求a在mod p下的乘法逆元可以求 $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为质数下的特殊情况，既然费马小定理可以，那么欧拉定理也可以。但使用乘法逆元求组合数的前提是模一个大质数，能用简单的费马小定理，为什么要用麻烦的欧拉定理呢？所以我在此只是提一提，并不想多讲。

刚刚介绍了三种求乘法逆元的方法，现在回到主题，如何求组合数。
还记得刚才的公式吗？
$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，且是个质数。

卢卡斯定理

先把这种方法的适用条件写在前面，适用于对答案模一个质数的情况。
和上面的乘法逆元求组合数的条件对比一下，只相差了一个大字。也就是说这个质数不用很大，非常小也行。
介绍一下卢卡斯定理，公式： $C_n^m=C_{n/p}^{m/p}*C_{n\%p}^{m\%p}$ （mod p），要求p为质数。
公式很好记吧，至于证明呢，看懂了也下次也不会证，会证了也没什么用（主要是我不会），我就不证明了。
观察公式，有没有发现，在n和m都小于p的时候，公式一点用都没有（ $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);
}

总结

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