组合数取模

求C(n,k）%p

p是素数 

1、1=<k<=10^6,1<=n<=10^18,求出mod p的分子a和分母b，a*b^(p-2)即为所求； FZU2020

2、1<=n,k<=10^18,p<=10^7,直接用lucas定理；hdu3037        hdu3944

p不是素数

3、1=<n,k<=10^6,分解因式再用快速幂  hdu3398

4、1<=n,k<=10^18，P=p1^c1 * p2^c2 * p3^c3 * … *pt ^ ct，pi为素数，1≤pi^ci≤10^7。

先对每一个pi^ci求出去除p因子的结果以及pi的个数，再用拓展欧几里德求逆元，最后再乘以pi，对每一个结果用中国剩余定理合并。详见ACblog   BZOJ2142   2012 Petr Mitrichev Contest 10-B

ACcode：

/**
FZU2020
*/
#include<cstdio>
#include<cstring>

typedef long long LL;
const int MS=10010;

LL res,b;
int n,m,p,t,r,T;
int c[MS],cnt;
int isprime[MS],pri[MS],top=0;

void prepare()
{
    for (int i=2;i<MS;i++)
    {
        if (!isprime[i])
        {
            pri[top++]=i;
            if (i<1000)
            for (int j=i*i;j<MS;j+=i)
            isprime[j]=1;
        }
    }
}

int quick_pow(LL x,int y,int mod)
{
    LL ans=1;
    while (y>0)
    {
        if (y&1) ans=(ans*x)%mod;
        x=(x*x)%mod;
        y>>=1;
    }
    return ans;
}

int main()
{
    prepare();
    scanf("%d",&T);
    while (T--)
    {
        scanf("%d%d%d",&n,&m,&p);
        if ((n-m)<m) m=n-m;
        res=b=1;
        for (int i=1;i<=m;i++)
        {
            res=(res*(n-i+1))%p;
            b=(b*i)%p;
        }
        res=(res*quick_pow(b,p-2,p))%p;
        printf("%I64d\n",res);
    }
    return 0;
}

/**
hdu3398
*/
#include<cstdio>
#include<cstring>

typedef long long LL;
const int NS=2000010;
const int MOD=20100501;

LL res;
int n,m,T;
int c[150000],cnt;
int isprime[NS],pri[150000],top=0;

void prime()
{
    for (int i=2;i<NS;i++)
    {
        if (!isprime[i])
        {
            pri[top++]=i;
            if (i<2000)
            for (int j=i*i;j<NS;j+=i)
            isprime[j]=1;
        }
    }
}

void factor(int k,int f)
{
    int t,r,i;
    for (i=0;i<top&&pri[i]<=k;i++)
    {
        t=k,r=pri[i];
        while (t>=r)
        {
            t/=r;
            c[i]+=t*f;
        }
    }
    cnt=cnt>i?cnt:i;
}

LL quick_pow(int a,int y)
{
    LL ans=1,x=a;
    while (y>0)
    {
        if (y&1) ans=(ans*x)%MOD;
        x=(x*x)%MOD;
        y>>=1;
    }
    return ans;
}

int main()
{
    prime();
    scanf("%d",&T);
    while (T--)
    {
        cnt=0,res=1;
        memset(c,0,sizeof(c));
        scanf("%d %d",&n,&m);
        factor(m,-1);
        factor(n+m,1);
        factor(n+1,-1);
        int t=n-m+1;
        for (int i=0;i<top&&pri[i]<=t;i++)
        {
            int r=pri[i];
            if (t%r==0)
            {
                while (t%r==0)
                c[i]++,t/=r;
            }
        }
        for (int i=0;i<cnt;i++)
        res=(res*quick_pow(pri[i],c[i]))%MOD;
        printf("%I64d\n",res);
    }
    return 0;
}

/**
hdu3037
*/
#include<cstdio>
#include<cstring>

typedef long long LL;

LL res,b;
int n,m,p,n0,m0,T;

LL quick_pow(LL x,int y,int mod)
{
    LL ans=1;
    while (y)
    {
        if (y&1) ans=(ans*x)%mod;
        x=(x*x)%mod;
        y>>=1;
    }
    return ans;
}

int main()
{
    scanf("%d",&T);
    while (T--)
    {
        scanf("%d%d%d",&n,&m,&p);
        res=1,n=n+m;
        while (n>0||m>0)
        {
            n0=n%p,m0=m%p;
            n/=p,m/=p;
            if (n0<m0)
            {
                res=0;
                break;
            }
            m0=(n0-m0)<m0?(n0-m0):m0;
            if (m0==0) continue;
            b=1;
            for (int i=1;i<=m0;i++)
            res=(res*(n0-i+1))%p,b=(b*i)%p;
            res=(res*quick_pow(b,p-2,p))%p;
        }
        printf("%I64d\n",res);
    }
    return 0;
}

/**
hdu3944
*/
#include<cstdio>
#include<cstring>

const int NS=9974;

int fac[NS][1229];
int isprime[NS],pri[NS],top=0;

void prepare()
{
    memset(isprime,0,sizeof(isprime));
    for (int i=2;i<NS;i++)
    {
        if (!isprime[i])
        {
            pri[top]=i,isprime[i]=top++;
            if (i<1000)
            for (int j=i*i;j<NS;j+=i)
            isprime[j]=1;
        }
    }
    for (int i=0;i<top;i++) fac[0][i]=1;
    for (int i=1;i<NS;i++)
    for (int j=0;j<top;j++)
    fac[i][j]=(fac[i-1][j]*i)%pri[j];
}

int quick_pow(int x,int y,int mod)
{
    int ans=1;
    for (;y;y>>=1)
    {
        if (y&1) ans*=x,ans%=mod;
        x*=x,x%=mod;;
    }
    return ans;
}

int combin(int n,int k,int p)
{
    int ans=1,x,y,q=isprime[p];
    while (k>0)
    {
        x=n%p,y=k%p;
        if (x<y) return 0;
        ans*=(fac[x][q]*quick_pow(fac[y][q]*fac[x-y][q]%p,p-2,p))%p;
        ans%=p;
        n/=p,k/=p;
    }
    return ans;
}

int main()
{
    prepare();
    int res,n,k,p,cas=0;
    while (~scanf("%d%d%d",&n,&k,&p))
    {
        if (n-k<k) res=k,n++,k++,k=n-k;
        else res=n-k,n++;
        res+=combin(n,k,p);
        res%=p;
        printf("Case #%d: %d\n",++cas,res);
    }
    return 0;
}

/**
BZOJ2142
*/
#include<cstdio>
#include<cstring>

typedef long long LL;
const int NS=100010;

LL mod,sum,res;
int n,m;
int p[NS],q[NS],w[10],ac[NS];
int a[NS],c[NS],d[NS],cnt;
int isprime[NS],pri[NS],top=0;

void prime_table()
{
    for (int i=2; i<NS; i++)
    {
        if (!isprime[i])
        {
            pri[top++]=i;
            if (i<1000)
                for (int j=i*i; j<NS; j+=i)
                    isprime[j]=1;
        }
    }
}

int eular(int x)
{
    int ret=1,i;
    for(i=2; i*i<=x; i++)
        if (x%i==0)
        {
            x/=i,ret*=i-1;
            while (x%i==0)
                x/=i,ret*=i;
        }
    if(x>1) ret*=x-1;
    return ret;
}


int quick_pow(int x,int y,int mo)
{
    int ans=1;
    while (y>0)
    {
        if (y&1) ans=((LL)ans*x)%mo;
        x=((LL)x*x)%mo;
        y>>=1;
    }
    return ans;
}

void divide(LL k)
{
    int t,r;
    for (int i=0; i<top&&pri[i]<=k; i++)
    {
        r=pri[i];
        if (k%r==0)
        {
            a[cnt]=r,ac[cnt]=1;
            while (k%r==0)
                c[cnt]++,k/=r,ac[cnt]*=r;
            cnt++;
        }
    }
}

int cal(int k,int pos,int flag)
{
    int pm=a[pos],y=ac[pos];
    int res=1,r=k%y,x=k/pm;
    d[pos]+=flag*x;
    for (int i=1; i<y; i++)
    {
        if (i%pm) res=((LL)res*i)%y;
    }
    res=quick_pow(res,k/y,y);
    for (int i=1; i<=r; i++)
    {
        if (i%pm) res=((LL)res*i)%y;
    }
    if (k>pm) res=((LL)res*cal(k/pm,pos,flag))%y;
    return res;
}

int exgcd(int a1,int a2,int &x,int &y)
{
    if (!a2)
    {
        x=1,y=0;
        return a1;
    }
    int r=exgcd(a2,a1%a2,x,y);
    int t=x;
    x=y,y=t-a1/a2*y;
    return r;
}

int main()
{
    int rem,res;
    prime_table();
    while (~scanf("%lld%d%d",&mod,&n,&m))
    {
        sum=cnt=0,res=1;
        memset(c,0,sizeof(c));
        memset(d,0,sizeof(d));
        for (int i=0; i<m; i++)
            scanf("%d",&w[i]),sum+=w[i];
        rem=n-sum;
        if (sum>n)
        {
            printf("Impossible\n");
            continue;
        }
        divide(mod);
        for (int i=0; i<cnt; i++)
        {
            p[i]=cal(n,i,1);
            q[i]=cal(rem,i,-1);
            for (int j=0; j<m; j++)
                q[i]=((LL)q[i]*cal(w[j],i,-1))%ac[i];
            p[i]=(LL)(p[i]*quick_pow(q[i],eular(ac[i])-1,ac[i])%ac[i])*quick_pow(a[i],d[i],ac[i])%ac[i];
        }
        res=0;
        for (int i=0; i<cnt; i++)
        {
            int x,y,r;
            r=exgcd(mod/ac[i],ac[i],x,y);
            res+=(LL)x*mod/ac[i]*p[i]%mod;
            res%=mod;
        }
        printf("%lld\n",((LL)res+mod)%mod);
    }
    return 0;
}

/**
petr2012_B
*/
#include<cstdio>
#include<cstring>
typedef long long LL;

const LL MOD2=1024;
const LL MOD5=9765625;
const LL MOD=10000000000LL;

LL n,m;
LL p[2],q[2],c[9];
int fac5[10000000],fac2[1111];

void prepare()
{
    fac5[0]=fac2[0]=1;
    for (int i=1;i<MOD2;i++)
    if (i&1) fac2[i]=((LL)fac2[i-1]*i)%MOD2;
    else fac2[i]=fac2[i-1];
    for (int i=1;i<MOD5;i++)
    if (i%5) fac5[i]=((LL)fac5[i-1]*i)%MOD5;
    else fac5[i]=fac5[i-1];
}

LL qpow(LL x,LL y,LL md)
{
    LL ans=1;
    for (;y;y>>=1)
    {
        if (y&1) ans*=x,ans%=md;
        x*=x,x%=md;
    }
    return ans;
}

LL cal(LL k,int flag,int z)
{
    LL fd,x,ans;
    int md,y,r;
    md=flag?MOD5:MOD2;
    r=k%md,x=k/md;
    if (flag) y=5,fd=fac5[md-1],ans=fac5[r];
    else y=2,fd=fac2[md-1],ans=fac2[r];
    ans=(ans*qpow(fd,x,md))%md;
    c[flag]+=k/y*z;
    if (k>y) ans=(ans*cal(k/y,flag,z))%md;
    return ans;
}

LL Exgcd(LL a,LL b,LL &x,LL &y)
{
    if (!b)
    {
        x=1,y=0;
        return a;
    }
    LL r=Exgcd(b,a%b,x,y);
    LL t=x;
    x=y,y=t-a/b*y;
    return r;
}

int main()
{
    freopen("combi.in","r",stdin);
    freopen("combi.out","w",stdout);
    prepare();
    while (~scanf("%I64d%I64d",&n,&m))
    {
        c[0]=c[1]=0;
        p[0]=cal(n,0,1);
        q[0]=cal(n-m,0,-1);
        q[0]*=cal(m,0,-1);
        p[1]=cal(n,1,1);
        q[1]=cal(n-m,1,-1);
        q[1]*=cal(m,1,-1);
        q[0]%=MOD2,q[1]%=MOD5;

        LL x,y,res;
        Exgcd(q[0],MOD2,x,y);
        p[0]=(p[0]*x%MOD2+MOD2)*qpow(2,c[0],MOD2)%MOD2;
        Exgcd(q[1],MOD5,x,y);
        p[1]=(p[1]*x%MOD5+MOD5)*qpow(5,c[1],MOD5)%MOD5;
        Exgcd(MOD2,MOD5,x,y);
        res=p[0]*(y*MOD5%MOD)%MOD;
        Exgcd(MOD5,MOD2,x,y);

        res=(res+p[1]*(y*MOD2%MOD)+MOD)%MOD;

        x=1;
        double s=1;
        int flag=0;
        if (n-m<m) m=n-m; //不加这句就Wrong在第4组，被坑的都没脾气了，求解释！！！
        for (LL i=1;i<=m;i++)
        {
            x=x*(n-i+1)/i;
            if (x>=MOD)
            {
                flag=1;
                break;
            }
        }
        if (flag) printf("...%010I64d\n",res);
        else printf("%I64d\n",res);
    }
    return 0;
}