BZOJ 1951 lucas定理 中国剩余定理

传送门

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

题意：求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))%4679、x4=sigma(C(N,i))%35617，若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;}