组合数的几种求法

一、利用杨辉三角求。
组合数有一个公式C(n,m)=C(n-1,m)+C(n-1,m-1);,仔细看,这不就是杨辉三角了吗?
但是这种只适合n,m都比较小的情况,不然会爆内存。

int a[2010][2010];
for(int i=0;i<2005;i++)a[i][0]=1;
for(int i=1;i<2005;i++){
    for(int j=1;j<=i;j++){
        a[i][j]=(a[i-1][j-1]+a[i-1][j])%1007;
    }
}

二、费马小定理
C(n,m)=n!/(m!*(n-m)!),求C(n,m)%MOD,除法运算的过程中无法直接求,最后取模的话容易超long long,这里可以用到费马小定理。
费马小定理:b^(p-1)%p==1 (p为素数)
费马小定理及其证明。
a/b%p = a / b * inv[a] * inv[b]%p = a * inv[b] % p
这样,除法就转化为乘法了,乘法满足同余定理,就不怕超long long了。
那怎么求inv[b]呢?
由费马小定理可知,b^(p-1)==1(mod p)
所以b*b^(p-2)==1(mod p)
所以inv[b]=b^(p-2)%p
我们可以快速幂求出来。
但是还是有点慢。
我们设 i = p/b,j = p%b。
所以p==0(mod p) ==> i*b+j == 0(mod p)
==> -i * b == j (mod p)
==> -i * inv[j] == inv[b] (左右同时乘以inv[b]*inv[j])
==>inv[b] == (p - p/b)*inv[p%b]%p
递推公式就求出来了。
贴代码。

ll inv(ll x){
    return x==1?1:(mod-mod/x)*inv(mod%x)%mod;
}
ll C(ll n,ll m){
    if(m<0||n<m)return 0;
    if(m>n-m)m=n-m;
    ll zi=1,mu=1;
    for(ll i=0;i<m;i++){
        zi=zi*(n-i)%mod;
        mu=mu*(i+1)%mod;
    }
    return zi*inv(mu)%mod;
}

三、Lucas定理
Lucas定理:A、B是非负整数,p是质数。A B写成p进制:A=a[n]a[n-1]…a[0],B=b[n]b[n-1]…b[0]。
则组合数C(A,B)与C(a[n],b[n])C(a[n-1],b[n-1])…*C(a[0],b[0]) mod p同余
即:Lucas(n,m,p)=C(n%p,m%p)*Lucas(n/p,m/p,p)

const int mod=1009;
ll inv(ll x){
    return x==1?1:(mod-mod/x)*inv(mod%x)%mod;
}
ll C(ll n,ll m){
    if(m<0||n<m)return 0;
    if(m>n-m)m=n-m;
    ll zi=1,mu=1;
    for(ll i=0;i<m;i++){
        zi=zi*(n-i)%mod;
        mu=mu*(i+1)%mod;
    }
    return zi*inv(mu)%mod;
}
ll lucas(ll n,ll m){
     ll ans=1;
     while(n&&m&&ans){
        ans=ans*C(n%mod,m%mod)%mod;
        n/=mod,m/=mod;
     }
     return ans;
}

四、分解质因数求组合数
BZOJ 1005 但是涉及到树的prufer编码。
树的prufer编码可以百度,树有唯一的编码,prufer编码也唯一对应一棵树。
那么这一道题就可以转化成prufer编码的排列组合数。
对于-1的点随便放置,就是剩余编码位置的k次方(k是-1点的种类)。
这里就用到了素数表的优化。毕竟大整数太慢了。

#include<bits/stdc++.h>
#define ll long long
#define lowbit(x) x&(-x)
#define eps 1e-8
#define lson l,m,td<<1
#define rson m+1,r,td<<1|1
#define INF 0x3f3f3f3f
#define line puts("");
#define LOCALR freopen("C:\\Users\\Administrator\\Desktop\\in.txt","r",stdin);
#define LOCALW freopen("C:\\Users\\Administrator\\Desktop\\out.txt","w",stdout);
using namespace std;
/*******************************************************************/
/*******************************************************************/
template<class T>
inline int read(T &x){
    char ch;
    bool flag=false;
    if((ch=getchar())==EOF)return -1;
    for(;!isdigit(ch);ch=getchar())if(ch=='-')flag=true;
    for(x=0;isdigit(ch);x=x*10+ch-'0',ch=getchar());
    x=flag?-x:x;
    return 1;
}
template<class T>
inline void write(T x,bool isnewline=false){
    static const int maxlen=100;
    static char s[maxlen];
    if(x<0){putchar('-');x=-x;}
    if(!x){putchar('0');return;}
    int len=0;
    for(;x;x/=10)s[len++]=x%10+'0';
    if(isnewline)puts(s);
    else for(int i=len-1;i>=0;i--)putchar(s[i]);
}
/*******************************************************************/
/*******************************************************************/
const int size=168;
const int prime[size+1]={0,2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,
79,83,89,97,101,103,107,109,113,127,131,137,139,149,151,157,163,167,173,179,181,191,
193,197,199,211,223,227,229,233,239,241,251,257,263,269,271,277,281,283,293,307,311,
313,317,331,337,347,349,353,359,367,373,379,383,389,397,401,409,419,421,431,433,439,
443,449,457,461,463,467,479,487,491,499,503,509,521,523,541,547,557,563,569,571,577,
587,593,599,601,607,613,617,619,631,641,643,647,653,659,661,673,677,683,691,701,709,
719,727,733,739,743,751,757,761,769,773,787,797,809,811,821,823,827,829,839,853,857
,859,863,877,881,883,887,907,911,919,929,937,941,947,953,967,971,977,983,991,997};
int d[10001],ans[size+1];
void C(int n,int m){
    for(int i=n-m+1;i<=n;i++){
        int p=i,j=1;
        while(p>1){
            while(p%prime[j]==0){
                ans[j]++;
                p/=prime[j];
            }
            j++;
        }
    }
    for(int i=1;i<=m;i++){
        int p=i,j=1;
        while(p>1){
            while(p%prime[j]==0){
                ans[j]--;
                p/=prime[j];
            }
            j++;
        }
    }
}
int main(){
    int n,temp,cnt=0;
    scanf("%d",&n);
    int len=n-2;
    for(int i=1;i<=n;i++){
        scanf("%d",&temp);
        if(temp==0||len-(temp-1)<0){printf("0");return 0;}
        if(temp==-1){cnt++;continue;}
        temp--;
        C(len,temp);
        len-=temp;
    }
    if(len){
        int j=1;
        while(cnt>1){
            while(cnt%prime[j]==0){
                ans[j]+=len;
                cnt/=prime[j];
            }
            j++;
        }
    }
    d[1]=1;
    int k=1;
    for(int i=1;i<=size;i++){
        while(ans[i]){
            ans[i]--;
            for(int j=1;j<=k;j++)d[j]*=prime[i];
            for(int j=1;j<=k;j++)if(d[j]>9){d[j+1]+=d[j]/10;d[j]%=10;}
            while(d[k+1]){k++;d[k+1]+=d[k]/10;d[k]%=10;}
        }
    }
    for(int i=k;i>0;i--)printf("%d",d[i]);
    return 0;
}

猜你喜欢

转载自blog.csdn.net/qq_41635132/article/details/81713031