求组合数的四种方法

求组合数

公式法

  • 公式
    Cab=Ca1b1+Ca1b
  • 公式理解
    在b个要抽取的物品中任选一个x物品,第一次抽取情况无非两种,抽到了x和没有抽到x,如果抽到了x那么就只需要在剩下的a-1个中再抽取b-1个即可,如果没有抽到x,就在剩下的a-1个中抽取b个,两次相加即为抽取的总种数。
  • 代码思路及实现
    我们只需要将数据范围n内的所有组合数都预处理出来,就可以实现求任一个组合数。时间复杂度o(n2),本质是一个从1到n的等差数列求和。
    给定n组数,每组包含一对a,b,求每一组Cab
#include<iostream>
using namespace std;
const int N=2010;
long long c[N][N];
const int mod=1e9+7;
int main()
{
    int n;
    cin>>n;
    for(int i=0;i<=2000;i++)
    {
        for(int j=0;j<=i;j++)
        {
            if(j==0)
            {
                c[i][j]=1;
            }
            else
            {
                c[i][j]=(c[i-1][j]+c[i-1][j-1])%mod;
            }
        }
    }
    while(n--)
    {
        int a,b;
        scanf("%d %d",&a,&b);
        printf("%ld\n",c[a][b]);
    }
    
    
    return 0;
}
  • 代码细节
    注意当b=0时的边界判定

定义法

  • 算法思路
    由组合数定义Cab=a(a1)(a2)(ab+1)b!=a!b!(ab)!,这种结果一般会很大所以会要求模上一个数,所以我们算出a!后就不去算除法而是用逆元来算。
  • 代码实现
    0ba105如果 用公式法时间复杂度o(n2)会超时,而定义法时间复杂度为o(nlogn)不会超时。
#include<iostream>
using namespace std;
const int mod=1e9+7;
const int N=1e5+10;
typedef long long LL;
LL fant[N],infant[N];//fant[i]代表i的阶乘,infant[i]代表i的阶乘的逆元
LL quick_mod(int a,int b,int p)//这道题mod是一个质数所以可以用快速幂来求
{
   LL res=1;
   while(b)
   {
       if(b&1)
       {
           res=res*a%p;
       }
       b>>=1;
       a=(LL)a*a%p;
   }
   return res;
}
int main()
{
   int n;
   cin>>n;
   fant[0]=infant[0]=1;
   for(int i=1;i<N;i++)
   {
       fant[i]=(LL)fant[i-1]*i%mod;
       infant[i]=(LL)infant[i-1]*quick_mod(i,mod-2,mod)%mod;//第29行
   }
   while(n--)
   {
       int a,b;
       scanf("%d %d",&a,&b);
       printf("%lld\n",(LL)fant[a]*infant[b]%mod*infant[a-b]%mod);
   }
   
   
   return 0;
}
  • 代码细节
    第29行求逆元的过程本质上就是abcab1c1(modp),即除以两个数同余先乘上一个数的逆元再乘上后一个数的逆元。

卢卡斯定理

  • 定理
    CabCa mod pb mod pCa/pb/p(modp)
  • 数学推导
    要用到一个公式(1+x)p1+xp(modp),这个公式的证明把左边展开即可:(1+x)p=Cp0+Cp1x+CP2x2++Cppxp因为从Cp1xcpp1xp1都是p的倍数模上一个p之后都为0,所以(1+x)p1+xp(modp)
    将a和b用p进制表示为:
    a=akpk+ak1pk1+ak2pk2++a0p0
    b=bkpk+bk1pk1+bk2pk2++b0p0
    对于(1+x)a(1+x)a=(1+x)akpk+ak1pk1+ak2pk2++a0p0=(1+x)akpk(1+x)ak1pk1(1+x)a0p0=(1+xpk)ak(1+xpk1)ak1(1+xp0)a0,由等式左边可以看出xb的系数为Cab,而xb=xbkpk+bk1pk1+bk2pk2++b0p0=(1+x)bkpk(1+x)bk1pk1(1+x)b0p0,而从等式右边可以看出xbkpk的系数为Cakbk以此类推得xb的系数Cab=CakbkCak1bk1Ca0b0(modp),因为a,b是用p进制表示的则有:a0=a mod p,b0=b mod p,a/p实际意义为a在p进制下向右移一位,b/p实际意义为b在p进制下向右移一位。那么对Ca/pb/p再进行一次上面的操作可以得到Ca/pb/p=CakbkCak1bk1Ca1b1(modp)
    综上有:Cab=Ca/pb/pCa mod pa mod p(modp)

-代码实现
卢卡斯算法时间复杂度对于0ba1018,0p105o(klogpnnlog2p)这里的数据组数k一般比较小

#include<iostream>
#include<algorithm>

using namespace std;
typedef long long LL;
int n;
int qmi(int a,int b,int p)
{
    int res=1;
    while(b)
    {
        if(b&1)
        {
            res=(LL)res*a%p;
        }
        b>>=1;
        a=(LL)a*a%p;
    }
    return res;
}
int C(LL a,LL b,int p)
{
    int res=1;
    for(int i=1,j=a;i<=b;i++,j--)
    {
        res=(LL)res*j%p;
        res=(LL)res*qmi(i,p-2,p)%p;
    }
    return res;
}
int lucas(LL a,LL b,int p)
{
    if(b>a)
    {
        return 0;
    }
    if(a<p)
    {
        return C(a,b,p);
    }
    return (LL)C(a%p,b%p,p)*lucas(a/p,b/p,p)%p;//a%p和b%p一定小于p所以直接用C来算组合数即可
}
int main()
{
    cin>>n;
    while(n--)
    {
        LL a,b;
        int p;
        scanf("%lld %lld %d",&a,&b,&p);
        printf("%d\n",lucas(a,b,p));
    }
    return 0;
}
  • 代码细节
    当a,b都小于p时候才用C函数来算,其余都继续递归直到小于p。

大精度计算

  • 针对类型
    可计算a,b较大且结果不模上一个数的组合数
  • 算法思路
    Cab=a!b!(ab)!对其分解质因数可得Cab=a!b!(ab)!=p1α1p2α2pnαn(p1p2pn)对于任一个pi它的次数都等于分子a!所含pi的个数减去分母b!(ab)!所含的pi的个数,求一个阶乘的某个质因子的个数有公式:cnt(a!)=ap+ap2+一直加到分母比a大时停止(因为分母比a大时向下取整为0加了和没加一样),则结果为Cab=i=1npiαi
  • 实现步骤
    先筛出2-a的所有质数,然后算出每个质数对应的次数,再算出i=1npiαi即可
  • 代码实现
    输入a,b求出Cab(0ba5000)
#include<iostream>
#include<algorithm>
#include<vector>

using namespace std;
typedef long long LL;
const int N=5010;
int primes[N];
int sum[N];
int cnt;
bool str[N];
int get(int x,int p)
{
    int res=0;
    while(x)
    {
        res+=x/p;
        x/=p;
    }
    return res;
}
vector<int> mul(vector<int> a,int b)
{
    vector<int> c;
    int t=0;
    for(int i=0;i<a.size();i++)
    {
        t+=a[i]*b;
        c.push_back(t%10);
        t/=10;
    }
    while(t)//这里得写成while,当b大于10的时候t在这里就可能大于10从而需要多次进位
    {
        c.push_back(t%10);
        t/=10;
    }
}
void get_primes(int n)
    return c;
{
    for(int i=2;i<=n;i++)
    {
        if(!str[i])
        {
            primes[cnt++]=i;
        }
        for(int j=0;i<=n/primes[j];j++)
        {
            str[i*primes[j]]=true;
            if(i%primes[j]==0)
            {
                break;
            }
        }
    }
}
int main()
{
    int a,b;
    cin>>a>>b;
    get_primes(a);
    for(int i=0;i<cnt;i++)
    {
        int p=primes[i];
        sum[p]=get(a,p)-get(b,p)-get(a-b,p);
    }
    vector<int> res;
    res.push_back(1);
    for(int i=0;i<cnt;i++)
    {
        for(int j=0;j<sum[primes[i]];j++)
        {
            res=mul(res,primes[i]);
        }
    }
    for(int i=res.size()-1;i>=0;i--)
    {
        printf("%d",res[i]);
    }
    return 0;
}

-代码细节
注意求大精度数相乘的模板


  1. 当p为质数时,有公式:Cpb0(modp)证明:Cpb=p!b!(pb)!=p(p1)!b(b1)!(pb)!=pbCp1b1由于p是质数,则b与p互质,b在模p条件下存在逆元,则有:Cpb=pinfact(b)Cp1b1(modp)显然是一个P的倍数则模上p后为0. ↩︎

  2. ap就表示a!中的a项乘积中是p的倍数的个数即a项乘积中能提出的p的个数,同理ap2就表示a!中的a项乘积中是p2的倍数的个数即a项乘积中能提出的p2的个数,以此类推就可以得到a!中含p的个数,这里你可能会有疑问,不应该加上二倍p2的个数,三倍p3的个数...才是所含p的个数吗,其实加上一倍就可以了,因为比如p2的倍数的个数在算p的倍数的时候就已经被算过了,毕竟p2的k倍数一定是p的k*p倍数。 ↩︎

posted @   Taco_gu  阅读(138)  评论(0编辑  收藏  举报
相关博文:
阅读排行:
· 在鹅厂做java开发是什么体验
· 百万级群聊的设计实践
· WPF到Web的无缝过渡:英雄联盟客户端的OpenSilver迁移实战
· 永远不要相信用户的输入:从 SQL 注入攻防看输入验证的重要性
· 浏览器原生「磁吸」效果!Anchor Positioning 锚点定位神器解析
点击右上角即可分享
微信分享提示