Lucas定理——大組合數取模
阿新 • • 發佈:2019-02-07
大組合數取模,求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;
}