卢卡斯定理大组合数取模及扩展

orz hhw posted @ 2015年6月16日 12:21 in 算法学习 with tags 数论 算法学习 卢卡斯 bzoj , 6397 阅读

卢卡斯定理内容：Lucas(n,m,p)=c(n%p,m%p)*Lucas(n/p,m/p,p)=n!/(m!*(n-m)!)*Lucas(n/p,m/p,p)

当p<10^5次时可以用卢卡斯定理进行组合数取模

HDOJ3037费马小定理版本:

#include<cstdio>
#define N 111111 
typedef long long LL;
LL n,m,p,T,i,fac[N];
LL pow(LL a,LL b){
	for(LL ans=1,tmp=a%p;b;b>>=1){
		if(b&1)ans=ans*tmp%p;
		tmp=tmp*tmp%p;
	}
	return ans;
}
LL C(LL n,LL m){
	if(m>n)return 0;
	return fac[n]*pow(fac[m]*fac[n-m],p-2)%p;//费马小定理求逆元
}
LL lucas(LL n,LL m){
	if(!m)return 1;
	return(C(n%p,m%p)*lucas(n/p,m/p))%p;//卢卡斯定理
}
int main(){
	scanf("%lld",&T);
	while(T--){
		scanf("%lld%lld%lld",&n,&m,&p);
		fac[0]=1;for(i=1;i<=p;i++)fac[i]=fac[i-1]*i%p;//阶乘预处理
		printf("%lld\n",lucas(n+m,m));
	}
}

不用费马小定理的模板：

#include<cstdio>
#define N 111111 
typedef long long LL;
LL n,m,p,T,i,x,y,fac[N];
void exgcd(LL a,LL b,LL &x,LL &y){
	if(!b){x=1;y=0;return;}
	exgcd(b,a%b,x,y);
	LL t=x;x=y;y=t-a/b*y;
}
LL C(LL n,LL m){
	if(m>n)return 0;
	exgcd(fac[m]*fac[n-m],p,x,y);
	x=(x+p)%p;
	return fac[n]*x%p;
}
LL lucas(LL n,LL m){
	if(!m)return 1;
	else return(C(n%p,m%p)*lucas(n/p,m/p))%p;
}
int main(){
	scanf("%lld",&T);
	while(T--){
		scanf("%lld%lld%lld",&n,&m,&p);
		fac[0]=1;for(i=1;i<=p;i++)fac[i]=fac[i-1]*i%p;
		printf("%lld\n",lucas(n+m,m));
	}
}

BZOJ1954：卢卡斯定理+中国剩余定理

Ans=G^(Σ(C(n,i)|i是n的约数))%999911659，然后组合数可以卢卡斯求出

然后这个指数太大了，快速幂都无法承受，有公式a^b mod c=a^(b mod phi(c) + phi(c))

phi(999911659)=999911658=2*3*4679*35617

求的时候模上面四个数，然后用中国剩余定理得到最后答案 中国剩余定理都忘了，没希望了

#include<cstdio>
#define M 999911659
typedef long long LL;
LL w[4]{2,3,4679,35617},i,j,x,y,ans,g,sum,d,n,a[4],fac[4][36666];
void exgcd(LL a,LL b,LL &x,LL &y){if(!b){x=1;y=0;return;}exgcd(b,a%b,y,x);y-=x*(a/b);}//扩展欧几里得
LL pow(LL a,LL b,LL p){for(sum=1,a%=p;b;a=a*a%p,b>>=1)if(b&1)sum=sum*a%p;return sum;}//快速幂
LL C(LL n,LL m,LL p){ return n<m?0:fac[p][n]*pow((fac[p][n-m]*fac[p][m]),w[p]-2,w[p])%w[p];}//求组合数
LL lucas(LL n,LL m,LL p){return m==0?1:lucas(n/w[p],m/w[p],p)*C(n%w[p],m%w[p],p)%w[p];}//卢卡斯
LL CRT(){//中国剩余定理
	for(ans=0,sum=M-1,i=0;i<4;i++){
		d=sum/w[i];
		exgcd(d,w[i],x,y);
		ans=(ans+d*x*a[i])%sum;
	}
	while(ans<=0)ans+=sum;
	return ans;
}
int main(){
	for(i=0;i<4;i++)for(fac[i][0]=1,j=1;j<=w[i];j++)fac[i][j]=(fac[i][j-1]*j)%w[i];//预处理阶乘
	scanf("%lld%lld",&n,&g);g%=M;
	for(i=1;i*i<=n;i++)if(n%i==0)
		for(j=0;j<4;j++){
			a[j]=(a[j]+lucas(n,i,j))%w[j];
			if(i*i!=n)a[j]=(a[j]+lucas(n,n/i,j))%w[j];
		}
	printf("%lld",pow(g,CRT(),M));
}

三道题都是27行，是不是觉得卢卡斯定理很简单呢。。

BZOJ2111更简单的题，一眼秒了，只要二叉树上统计排列组合就行了，然后奇怪的WA了，而且改成I64d就PE了，但是本地过了，应该是数据youginzi

#include<cstdio>
typedef long long LL;
int n,i;LL x,y,p,ans,size[2000002],fac[2000002];
void exgcd(LL a,LL b,LL &x,LL &y){
	if(!b){x=1;y=0;return;}
	exgcd(b,a%b,y,x);y-=x*(a/b);
}
LL C(LL a,LL b){
	if(!a||!b)return 1;
	exgcd(fac[a]*fac[b],p,x,y);if(x<0)x+=p;
	return fac[a+b]*x%p;
}
int main(){
	scanf("%d%lld",&n,&p);
	for(fac[1]=1,i=2;i<=n*2;i++)fac[i]=fac[i-1]*i%p;
	for(ans=1,i=n;i;i--){
		ans=(ans*C(size[i*2],size[i*2+1]))%p;
		size[i]=size[i*2]+size[i*2+1]+1;
	}
	printf("%lld",ans);
}

扩展lucas一直是个坑，今天看到BZOJ上有道题只能这么做，就学习了下

大概就是和求(n!)%p^k时，把所有p的倍数的放到一边，剩余的以p^k为循环节用快速面求答案，然后把p的倍数继续递归处理，这样的时间复杂度是线性的，最后合并的时候CRT合并一下就行了

#include<bits/stdc++.h>
#define N 100100
typedef long long LL;LL ans,a[N];int M,n,m,i,jj,l,tot,p[N],P[N],w[9];struct U{LL x,z;};
void exgcd(LL a,LL b,LL &x,LL &y){if(!b){x=1;y=0;return;}exgcd(b,a%b,y,x);y-=x*(a/b);}
LL pow(LL a,LL b,LL p){LL sum=1;for(a%=p;b;a=a*a%p,b>>=1)if(b&1)sum=sum*a%p;return sum;}
LL inv(LL a,LL p){LL x,y;exgcd(a,p,x,y);return(x+p)%p;}
U fac(int k,LL n){
	if (!n)return U{1,0};
	LL x=n/p[k],y=n/P[k],ans=1;int i;
	if(y){
		for(i=2;i<P[k];i++)if(i%p[k])ans=ans*i%P[k];
	    ans=pow(ans,y,P[k]);
	}
	for(i=y*P[k];i<=n;i++)if(i%p[k])ans=ans*i%M;
	U z=fac(k,x);return U{ans*z.x%M,x+z.z};
}
LL get(int k,LL n,LL m){
	U a=fac(k,n),b=fac(k,m),c=fac(k,n-m);
	return pow(p[k],a.z-b.z-c.z,P[k])*a.x%P[k]*inv(b.x,P[k])%P[k]*inv(c.x,P[k])%P[k];
}
LL CRT(){
	LL d,w,y,x,ans=0;
	for(int i=1;i<=l;i++)w=M/P[i],exgcd(w,P[i],x,y),ans=(ans+w*x*a[i])%M;
	return(ans+M)%M;
}
LL C(LL n,LL m){
	for(int i=1;i<=l;i++)a[i]=get(i,n,m);
	return CRT();
}
int main(){
	for(scanf("%d%d%d",&M,&n,&m),jj=M,i=2;i*i<=jj;i++)if(jj%i==0)
		for(p[++l]=i,P[l]=1;jj%i==0;P[l]*=p[l])jj/=i;
	if(jj>1)l++,p[l]=P[l]=jj;
	for(i=1;i<=m;i++)scanf("%d",&w[i]),tot+=w[i];
	if(tot>n)return puts("Impossible"),0;
	for(ans=C(n,tot),i=1;i<=m;i++)ans=ans*C(tot,w[i])%M,tot-=w[i];
	printf("%lld",ans);
}

扩展lucas求卡特兰数：

#include<bits/stdc++.h>
#define N 1001000
typedef long long LL;LL n,w,ans,M,i,jj,l,a[N],p[N],P[N];struct U{LL x,z;};
LL pow(LL a,LL b,LL p){LL sum=1;for(a%=p;b;a=a*a%p,b>>=1)if(b&1)sum=sum*a%p;return sum;}
U fac(LL k,LL n){
	if(!n)return U{1,0};
	LL x=n/p[k],y=n/P[k],ans=1;LL i;
	if(y){
		for(i=2;i<P[k];i++)if(i%p[k])ans=ans*(i%P[k])%P[k];
	    ans=pow(ans,y,P[k]);
	}
	for(i=y*P[k];i<=n;i++)if(i%p[k])ans=ans*(i%M)%M;
	U z=fac(k,x);return U{ans*z.x%M,x+z.z};
}
LL get(LL k,LL n){
	U a=fac(k,n*2),b=fac(k,n),c=fac(k,n+1);
	return pow(p[k],a.z-b.z-c.z,P[k])*a.x%P[k]*pow(b.x,P[k]/p[k]*(p[k]-1)-1,P[k])%P[k]*pow(c.x,P[k]/p[k]*(p[k]-1)-1,P[k])%P[k];
}
void C(LL n){for(LL i=1;i<=l;i++)w=M/P[i],ans=(ans+w*get(i,n)%M*pow(w,P[i]/p[i]*(p[i]-1)-1,P[i]))%M;}
int main(){
	for(scanf("%lld%lld",&n,&M),jj=M,i=2;i*i<=jj;i++)if(jj%i==0)
		for(p[++l]=i,P[l]=1;jj%i==0;P[l]*=p[l])jj/=i;
	if(jj>1)l++,p[l]=P[l]=jj;C(n-1);
	printf("%lld",ans);
}

BZOJ3129 容斥把所有符号转成≥，然后用隔板法，膜数不是质数，用扩展lucas求组合数

#include<bits/stdc++.h>
using namespace std;
typedef long long LL;LL ans,a[7];int T,M,x,n,m,i,j,l,n1,n2,tot,o[9],p[9],P[9],w[9];struct U{LL x,z;};
LL pow(LL a,int b,int p){LL sum=1;for(a%=p;b;a=a*a%p,b>>=1)if(b&1)sum=sum*a%p;return sum;}
U fac(int k,int n){
	if (!n)return U{1,0};
	LL x=n/p[k],y=n/P[k],ans=1,i;
	for(i=2;i<P[k];i++)if(i%p[k])ans=ans*i%P[k];
	for(ans=pow(ans,y,P[k]),i=y*P[k];i<=n;i++)if(i%p[k])ans=ans*i%M;
	U z=fac(k,x);return U{ans*z.x%M,x+z.z};
}
LL get(int k,int n,int m){
	U a=fac(k,m),b=fac(k,n),c=fac(k,m-n);
	return pow(p[k],a.z-b.z-c.z,P[k])*a.x%P[k]*pow(b.x,P[k]/p[k]*(p[k]-1)-1,P[k])%P[k]*pow(c.x,P[k]/p[k]*(p[k]-1)-1,P[k])%P[k];
}
LL CRT(){
	LL d,w,y,x,ans=0,i;
	for(i=1;i<=l;i++)w=M/P[i],ans=(ans+w*(a[i]+M)%M*pow(w,P[i]/p[i]*(p[i]-1)-1,P[i]))%M;
	return(ans+M)%M;
}
LL C(int n,int m){for(int i=1;i<=l;i++)a[i]=get(i,n,m);return CRT();}
void dfs(int x,int y,int m){
	if(m<n)return;
	if(x>n1){
		if(!y)ans=(ans+C(n-1,m-1))%M;else ans=(ans-C(n-1,m-1)+M)%M;
		return;
	}
	dfs(x+1,y,m);dfs(x+1,y^1,m-o[x]);
}
int main(){
	for(scanf("%d%d",&T,&M),j=M,i=2;i*i<=j;i++)if(j%i==0)for(p[++l]=i,P[l]=1;j%i==0;P[l]*=p[l])j/=i;
	if(j>1)l++,p[l]=P[l]=j;
	for(;T--;ans=0){
		scanf("%d%d%d%d",&n,&n1,&n2,&m);
		for(i=1;i<=n1;i++)scanf("%d",&o[i]);
		for(;n2--;m-=x-1)scanf("%d",&x);
		dfs(1,0,m);printf("%lld\n",ans);
	}
}
Small_feed 评论 (4)
1
[回复]
Avatar_small
nbdhhzh 说:
2015年6月17日 06:47

每天扩展欧几里得用用不知道想些什么。。费马小定理不知道高明到哪里去了

2
[回复]
Avatar_small
nbdhhzh 说:
2015年6月17日 06:47

直接学扩展卢卡斯啊。。卢卡斯有卵用。。

3
[回复]
Avatar_small
orzhzh 说:
2015年6月18日 15:27

不学卢卡斯怎么学扩展卢卡斯。。。

4
[回复]
Avatar_small
orzorzhzh 说:
2016年4月16日 04:43

不学卢卡斯怎么学扩展卢卡斯+1


登录 *


loading captcha image...
(输入验证码)
or Ctrl+Enter