1951: [Sdoi2010]古代猪文

#include<cstdio>
using namespace std;
typedef long long ll;
const int P=999911659;
int M[5],t[5]={2,3,4679,35617};
int N,G,fac[5][(int)4e4];
int fpow(ll a,ll p,ll mod){
    ll res=1;
    for(;p;p>>=1,a=a*a%mod) if(p&1) res=res*a%mod;
    return res;
}
int C(int n,int m,int x){
    if(m>n) return 0;
    return fac[x][n]*fpow(fac[x][n-m]*fac[x][m],t[x]-2,t[x])%t[x];
}
void exgcd(int a,int 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;
}
int lucas(int a,int b,int x){
    if(!b) return 1;
    return C(a%t[x],b%t[x],x)*lucas(a/t[x],b/t[x],x)%t[x];
} 
ll CRT(){
    ll x0,y0,ans=0,Z=P-1;
    for(int i=0;i<4;i++){
        int d=Z/t[i];
        exgcd(d,t[i],x0,y0);
        ans=(ans+d*x0*M[i])%Z;
    }
    while(ans<=0) ans+=Z;
    return ans;
}
int main(){
    scanf("%d%d",&N,&G);
    if(N==G){puts("0");return 0;}
    G%=P;
    for(int i=0;i<4;i++){
        fac[i][0]=1;
        for(int j=1;j<=t[i];j++){
            fac[i][j]=(fac[i][j-1]*j)%t[i];
        }
    }
    for(int i=1,tmp;i*i<=N;i++){
        if(N%i==0){
            tmp=N/i;
            for(int j=0;j<4;j++){
                if(tmp!=i) M[j]=(M[j]+lucas(N,i,j))%t[j];
                M[j]=(M[j]+lucas(N,tmp,j))%t[j];
            }
        }
    }
    printf("%d",fpow(G,CRT(),P));
    return 0;
}