1. 程式人生 > >Lucas定理——大組合數取模

Lucas定理——大組合數取模

大組合數取模,求C[n][m]%p

公式:C[n][m]%p == C[n%p][m%p]*C[n/p][m/p]%p

注意,Lucas的要求是n,m<=10^5,如果n,m>=10^5,那麼要求p<=10^5

楊輝三角:

f[0][0]=f[1][0]=f[1][1]=1;
for(int i=2;i<=n;i++){
    f[i][0]=1;
    for(int j=1;j<=m;j++)
        f[i][j]=f[i-1][j-1]+f[i-1][j];
}
//f[p][q]=C[p][q];

程式碼(內含粗略解釋):

#define LL long long
LL mont(LL a,LL b,LL c){//快速冪
LL ans=1; a%=c; while(b){ if(b&1)ans=ans*a%c; b>>=1,a=a*a%c; } return ans%c; } LL C(LL a,LL b,LL p){//p一定為質數 if(a<b)return 0;//不存在 if(a==b)return 1; if(b>a-b)b=a-b;//C[n][m]=C[n][n-m] LL A=1,B=1; for(LL i=0;i<b;i++){ A=(A*(
a-i))%p; //求(a*(a-1)*(a-2)*...*(a-b+1))%p //即求(a!/(a-b)!)%p B=(B*(b-i))%p; //求(b*(b-1)*(b-2)*...*2)%p //即求b!%p } return (A*mont(B,p-2,p))%p; //根據費馬小定理,a^(p-1)%p==1(p為質數) //所以,b^(p-1)%p==1,b*b^(p-2)%p==1 //又因為,b*b^(-1)==1,即b*b^(-1)%p==1(此處指b的逆元) //所以,b^(-1)==b^(p-2),用快速冪解決
} LL Lucas(LL n,LL m,LL p){//C[n][m]%p if(m==0)return 1; return C(n%p,m%p,p)*Lucas(n/p,m/p,p)%p; }

把註釋去掉就是這個樣子↓

#define LL long long
LL mont(LL a,LL b,LL c){
    LL ans=1;
    a%=c; 
    while(b){ 
        if(b&1)ans=ans*a%c; 
        b>>=1,a=a*a%c; 
    } 
    return ans%c; 
}
LL C(LL a,LL b,LL p){
    if(a<b)return 0;
    if(a==b)return 1; 
    if(b>a-b)b=a-b;
    LL A=1,B=1; 
    for(LL i=a-b+1;i<=a;i++)A=A*i%p;
    for(LL i=1;i<=b;i++)B=B*i%p;
    return (A*mont(B,p-2,p))%p; 
}
LL Lucas(LL n,LL m,LL p){//C[n][m]%p 
    if(m==0)return 1; 
    return C(n%p,m%p,p)*Lucas(n/p,m/p,p)%p; 
}

如果p不是很大的話,可以事先求好每個數的階乘和逆元,直接使用。

void pre(){
    Fact[0]=1,Inv[0]=Inv[1]=1;
    for(int i=1;i<=p;i++)Fact[i]=(Fact[i-1]*i)%p;
    for(int i=2;i<=p;i++)Inv[i]=(p-p/i)*Inv[p%i]%p;
    for(int i=1;i<=p;i++)Inv[i]=Inv[i]*Inv[i-1]%p;
}
ll Lucas(ll n,ll m){
    if(m>n)return 0;
    if(!m)return 1;
    if(m<p && n<p)return Fact[n]*Inv[n-m]%p*Inv[m]%p;
    return Lucas(n%p,m%p)*Lucas(n/p,m/p)%p;
}