首页 > 代码库 > 扩展lucas定理
扩展lucas定理
i64 POW(i64 a,i64 b,i64 mod) { i64 ans=1; while(b) { if(b&1) ans=ans*a%mod; a=a*a%mod; b>>=1; } return ans; } i64 POW(i64 a,i64 b) { i64 ans=1; while(b) { if(b&1) ans=ans*a; a=a*a; b>>=1; } return ans; } i64 exGcd(i64 a,i64 b,i64 &x,i64 &y) { i64 t,d; if(!b) { x=1; y=0; return a; } d=exGcd(b,a%b,x,y); t=x; x=y; y=t-a/b*y; return d; } bool modular(i64 a[],i64 m[],i64 k) { i64 d,t,c,x,y,i; for(i=2;i<=k;i++) { d=exGcd(m[1],m[i],x,y); c=a[i]-a[1]; if(c%d) return false; t=m[i]/d; x=(c/d*x%t+t)%t; a[1]=m[1]*x+a[1]; m[1]=m[1]*m[i]/d; } return true; } i64 reverse(i64 a,i64 b) { i64 x,y; exGcd(a,b,x,y); return (x%b+b)%b; } i64 C(i64 n,i64 m,i64 mod) { if(m>n) return 0; i64 ans=1,i,a,b; for(i=1;i<=m;i++) { a=(n+1-i)%mod; b=reverse(i%mod,mod); ans=ans*a%mod*b%mod; } return ans; } i64 C1(i64 n,i64 m,i64 mod) { if(m==0) return 1; return C(n%mod,m%mod,mod)*C1(n/mod,m/mod,mod)%mod; } i64 cal(i64 n,i64 p,i64 t) { if(!n) return 1; i64 x=POW(p,t),i,y=n/x,temp=1; for(i=1;i<=x;i++) if(i%p) temp=temp*i%x; i64 ans=POW(temp,y,x); for(i=y*x+1;i<=n;i++) if(i%p) ans=ans*i%x; return ans*cal(n/p,p,t)%x; } i64 C2(i64 n,i64 m,i64 p,i64 t) { i64 x=POW(p,t); i64 a,b,c,ap=0,bp=0,cp=0,temp; for(temp=n;temp;temp/=p) ap+=temp/p; for(temp=m;temp;temp/=p) bp+=temp/p; for(temp=n-m;temp;temp/=p) cp+=temp/p; ap=ap-bp-cp; i64 ans=POW(p,ap,x); a=cal(n,p,t); b=cal(m,p,t); c=cal(n-m,p,t); ans=ans*a%x*reverse(b,x)%x*reverse(c,x)%x; return ans; } //计算C(n,m)%mod i64 Lucas(i64 n,i64 m,i64 mod) { i64 i,t,cnt=0; i64 A[205],M[205]; for(i=2;i*i<=mod;i++) if(mod%i==0) { t=0; while(mod%i==0) { t++; mod/=i; } M[++cnt]=POW(i,t); if(t==1) A[cnt]=C1(n,m,i); else A[cnt]=C2(n,m,i,t); } if(mod>1) { M[++cnt]=mod; A[cnt]=C1(n,m,mod); } modular(A,M,cnt); return A[1]; }
扩展lucas定理
声明:以上内容来自用户投稿及互联网公开渠道收集整理发布,本网站不拥有所有权,未作人工编辑处理,也不承担相关法律责任,若内容有误或涉及侵权可进行投诉: 投诉/举报 工作人员会在5个工作日内联系你,一经查实,本站将立刻删除涉嫌侵权内容。