卢卡斯定理大组合数取模及扩展
卢卡斯定理内容: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); } }
2015年6月17日 06:47
每天扩展欧几里得用用不知道想些什么。。费马小定理不知道高明到哪里去了
2015年6月17日 06:47
直接学扩展卢卡斯啊。。卢卡斯有卵用。。
2015年6月18日 15:27
不学卢卡斯怎么学扩展卢卡斯。。。
2016年4月16日 04:43
不学卢卡斯怎么学扩展卢卡斯+1