BZOJ 1951 lucas定理 中国剩余定理

xiaoxiao2021-02-28  122

传送门

照旧记录一下,不然担心以后碰到类似的题依旧不知道怎么做……

题意:G^(sigma(C(N,i)))%MOD,MOD=999911659,i是N的因子(包括1和本身)

思路:首先幂次非常大,所以先降幂次,即sigma(C(N,i))%(MOD-1)+(MOD-1),那么关键就是求sigma(C(N,i))%(MOD-1),我们可以遍历能整除N的数,但999911658是一个合数,它可以分解为四个素数的乘积,即999911658=2*3*4679*35617,那么就可以分别求x1=sigma(C(N,i))%2、x2=sigma(C(N,i))%3、x3=sigma(C(N,i))F79、x4=sigma(C(N,i))5617,若x=sigma(C(N,i))%999911658,那么x%x1=2、x%x2=3、x%x3=4679、x%x4=35617,那么就可以用中国剩余定理求得x了。可以发现素数都比较小,所以可以阶乘预处理。这个题还有个需要注意的地方,就是G=999911659的时候结果就是0了,可是我测试了几组数据,发现不写这个,最后结果也都是0,不知道哪组数据不是0了……

完整代码:

#include <iostream> #include <string.h> #include <stdio.h> using namespace std; typedef long long LL; LL mod[4]={2,3,4679,35617}; const int MOD=999911659; LL fac[4][40000]; LL a[4]; LL quick_mod(LL a,LL b,LL p) { LL ans=1; a%=p; while(b) { if(b&1) { ans=ans*a%p; b--; } b>>=1; a=a*a%p; } return ans; } void init() { for(int i=0;i<4;i++) fac[i][0]=1; for(int i=0;i<4;i++) for(int j=1;j<mod[i];j++) fac[i][j]=fac[i][j-1]*j%mod[i]; } LL C(LL n,LL m,LL p,int cnt) { if(m>n) return 0; LL ans=fac[cnt][n]*quick_mod(fac[cnt][m]*fac[cnt][n-m],p-2,p)%p; return ans; } LL Lucas(LL n, LL m,LL p,int cnt) { if(m==0) return 1; return C(n%p,m%p,p,cnt)*Lucas(n/p,m/p,p,cnt); } LL gcd(LL a,LL b) { if(b==0) return a; else gcd(b,a%b); } void exgcd(LL a,LL b,LL &x,LL &y) { if(b==0) { x=1; y=0; return ; } LL x1,y1; exgcd(b,a%b,x1,y1); x=y1; y=x1-(a/b)*y1; } LL CRT(LL a[]) { LL M = 1; LL ans = 0; for(LL i=0; i<4; i++) M *= mod[i]; for(LL i=0; i<4; i++) { LL x, y; LL Mi=M/mod[i]; exgcd(Mi,mod[i],x,y); ans=(ans+Mi*x*a[i])%M; } if(ans<0) ans += M; return (ans+M)%M; } int main() { init(); LL N,G; while(scanf("%lld%lld",&N,&G)!=-1) { for(int i=0;i<4;i++) a[i]=0; if(G==MOD) { printf("0\n"); continue; } for(int i=1;i*i<=N;i++) { if(N%i==0) { LL x=i; a[0]=(a[0]+Lucas(N,x,mod[0],0))%mod[0]; a[1]=(a[1]+Lucas(N,x,mod[1],1))%mod[1]; a[2]=(a[2]+Lucas(N,x,mod[2],2))%mod[2]; a[3]=(a[3]+Lucas(N,x,mod[3],3))%mod[3]; x=N/i; if(N!=i*i) { a[0]=(a[0]+Lucas(N,x,mod[0],0))%mod[0]; a[1]=(a[1]+Lucas(N,x,mod[1],1))%mod[1]; a[2]=(a[2]+Lucas(N,x,mod[2],2))%mod[2]; a[3]=(a[3]+Lucas(N,x,mod[3],3))%mod[3]; } } } LL temp=CRT(a)+MOD-1; LL ans=quick_mod(G,temp,MOD)%MOD; cout<<ans<<endl; } return 0; }

转载请注明原文地址: https://www.6miu.com/read-34942.html

最新回复(0)