卢卡斯算法
原理
用于求大组合数
c(n,m)%mod=C(n%mod,m%mod)*C(n/mod,m/mod)%mod
代码
递归版
LL fac[mod+105],inv[mod+105];
void getInv()
{
fac[0]=fac[1]=inv[1]=1;
for(int i=2;i<mod;i++)
{
fac[i]=fac[i-1]*i%mod;
inv[i]=(mod-mod/i)*inv[mod%i]%mod;
}
}
LL C(LL n,LL m)
{
if(m>n)return 0;
return fac[n]*inv[fac[m]*fac[n-m]%mod]%mod;
}
LL Lucas(LL n,LL m)
{
if(m>n)return 0;
if(n<mod && m<mod)return C(n,m);
return Lucas(n/mod,m/mod)*C(n%mod,m%mod)%mod;
}
非递归版
LL Lucas(LL n,LL m)
{
if(m>n)return 0;
LL ret=1;
for(;m;n/=mod,m/=mod)//m会被先除完,当m为0的时候答案也就不会再变化了
ret=ret*C(n%mod,m%mod)%mod;
return ret;
}
性能分析
- 时间复杂度:O(logn+p)
- 适用范围:mod是不大的素数的大组合数计算
扩展卢卡斯算法
原理
把mod数分解质因数,计算p1^t1 , p2^t2 …. pm^tm,然后用中国剩余定理合并
C(n.m)mod(pi^ti)的计算方法就是先算阶乘然后用逆元搞一搞就出来了,其中单独计算pi的倍数用快速幂算,会出现一个阶乘形式递归处理,还有一个分板块的,以mod数位基础,除去pi的倍数,然后用快速幂计算每一块,剩下的不成块的单独计算。所以可以看作3部分或者4部分
代码
typedef long long LL;
LL exgcd(LL a,LL b,LL &x,LL &y)
{
if(!b)
{
x=1,y=0;
return a;
}
LL d=exgcd(b,a%b,y,x);
y-=a/b*x;
return d;
}
LL inv(LL a,LL mod)
{
LL d,x,y;
d=exgcd(a,mod,x,y);
return d==1?(x%mod+mod)%mod:-1;
}
LL qkpow(LL a,LL p,LL mod)
{
LL t=1,tt=a%mod;
while(p)
{
if(p&1)t=t*tt%mod;
tt=tt*tt%mod;
p>>=1;
}
return t;
}
LL fac(LL n,LL p,LL pt)//计算除去因数p的阶乘
{
if(n==0)return 1;//注意边界
LL ret=1;
if(n/pt)
{
for(LL i=1;i<=pt;i++)//pt为一组
if(i%p)ret=ret*i%pt;//先计算分块部分,注意判断i的性质
ret=qkpow(ret,n/pt,pt);//再将块累加
}
for(LL i=1;i<=n%pt;i++)//因为n可能很大,所以我们这样做
if(i%p)ret=ret*i%pt;//计算单出来的一部分
return ret*fac(n/p,p,pt)%pt;//递归计算
}
LL C(LL n,LL m,LL p,LL pt)//扩展卢卡斯定理计算C(n,m)%(p^t)的值
{
LL a,b,c,d,k=0;//K用来记录上面的p-下面的p有多少个
for(LL t=n;t;t/=p)k+=t/p;
for(LL t=m;t;t/=p)k-=t/p;
for(LL t=n-m;t;t/=p)k-=t/p;//这三排是记录P的个数,这样除的好处是不会炸LL,和分解n!因数是一样的,最终分解出来个数相同
a=qkpow(p,k,pt);//k肯定是正数,因为组合数是个正数
b=fac(n,p,pt);
c=fac(m,p,pt);
d=fac(n-m,p,pt);
return a*b%pt*inv(c,pt)%pt*inv(d,pt)%pt;
}
LL exLucas(LL n,LL m,LL mod)//计算C(n,m)%mod的值+中国剩余定理
{
LL M=mod,pt,ans=0;//中国剩余定理无需和
for(LL i=2;i<=mod;i++)//分解mod
if(mod%i==0)
{
pt=1;//p^t
while(mod%i==0)mod/=i,pt*=i;
ans=(ans+C(n,m,i,pt)*(M/pt)%M*inv(M/pt,pt)%M)%M;//中国剩余定理合并
}
return ans;
}
int main()
{
LL n,m,p;
while(cin>>n>>m>>p)
cout<<exLucas(n,m,p)<<endl;//千万不要写C(n%p,m%p,p)
return 0;
}
注意:
- 所有变量全部开成LL
- 不要随意改代码,每个细节都有可能爆LL
性能分析
- 时间复杂度:O(log^2n*mod)
- 适用范围:对于mod数不是很大但是不是素数的。
- 弱化版:对于给定mod数而且分解质因数后每个质因数都只有一个直接种CRT+Lucas就可以了。