数学的东西(BZOJ1951)

#include <cstdio>
#define LL long long 

  LL finmo=999911659;
  LL fac[4][40001],inv[4][40001];
  LL tmp[4],rev[4];
  LL n,g,x,y;
  const LL mo[4]={2,3,4679,35617};

  LL qpow(LL bas,LL pow,LL mo){
      LL ret=1;
      for (;pow;bas*=bas,bas%=mo,pow=pow>>1){
        if (pow&1) ret*=bas,ret%=mo;
    }
    return(ret);
  }//快速幂
  
  LL C(LL t1,LL t2,LL mopo){
      if (t2>t1) return(0);
      return((fac[mopo][t1]*inv[mopo][t2])%mo[mopo]*inv[mopo][t1-t2]%mo[mopo]);
  }//组合数
  
  LL lucas(int t1,int t2,int mopo){
      LL ret=1;
      while(t1||t2){
        ret*=C(t1%mo[mopo],t2%mo[mopo],mopo);ret%=mo[mopo];
      t1/=mo[mopo];t2/=mo[mopo];    
    }
    return(ret);
  }//lucas定理
  
  LL solve(LL num){
      for (int i=0;i<=3;i++) tmp[i]=lucas(n,num,i);
    LL ret=0;
    for (int i=0;i<=3;i++) 
      ret+=tmp[i]*(finmo/mo[i])*rev[i],ret%=finmo-1;
    return(ret);    
  }//线性同余方程的解
  
  void exgcd(LL a,LL b,LL &x,LL &y){
      if (b==0){
        x=1LL;y=0LL;return;    
    }
    exgcd(b,a%b,x,y);
    LL t=x;x=y;y=t-(a/b)*y;
  }//扩展欧几里得

  int main(){      
      scanf("%lld%lld",&n,&g);
      for (int i=0;i<=3;i++){
        fac[i][0]=1;inv[i][0]=1;
      for (int j=1;j<mo[i];j++) {fac[i][j]=(fac[i][j-1]*j)%mo[i];inv[i][j]=qpow(fac[i][j],mo[i]-2,mo[i]);}    
    }
    for (int i=0;i<=3;i++){
      exgcd(finmo/mo[i],mo[i],x,y);
      rev[i]=(x%mo[i]+mo[i])%mo[i];    
    }
    
    LL ans=0;
    for (int i=1;i*i<=n;i++)
     if (n%i==0){
      ans+=solve(i);ans%=finmo-1;    
      if (n/i!=i) ans+=solve(n/i),ans%=finmo-1;
    }
    
    LL finans=qpow(g,ans,finmo);
    if (g==999911659)printf("0
");else printf("%lld
",finans);
  }