【BZOJ 1951】 [Sdoi2010]古代猪文

 1 #include<cstdio>
 2 #define ll long long
 3 const int P=999911659;
 4 int t[4]={2,3,4679,35617};
 5 int fac[5][40000],M[5],G,N;
 6 int pow(ll a,ll b,ll x){
 7     ll ans=1;
 8     while (b>0){
 9         if (b&1==1) ans=(ans*a)%x;
10         b=b>>1;a=(a*a)%x;
11     }
12     return ans;
13 }
14  
15 int C(int a,int b,int x){
16     if (b>a) return 0;
17     return fac[x][a]*pow(fac[x][a-b]*fac[x][b],t[x]-2,t[x])%t[x];
18 }
19  
20 int lucas(int a,int b,int x){
21     if(!b) return 1;
22     return C(a%t[x],b%t[x],x)*lucas(a/t[x],b/t[x],x)%t[x];
23 }
24  
25 void exgcd(int a,int b,ll &x,ll &y){
26     if (b==0){x=1,y=0;return;}
27     exgcd(b,a%b,x,y);
28     int t=x;x=y,y=t-(a/b)*y;
29 }
30  
31 ll CRT()  
32 {  
33     ll N=1;  
34     ll i,d,x0,y0,ans=0;  
35     for(i=0;i<4;i++)  
36         N*=t[i];
37     for(i=0;i<4;i++)  
38     {  
39         d=N/t[i];  
40         exgcd(d,t[i],x0,y0);  
41         ans=(ans+d*x0*M[i])%N;  
42     }  
43     while(ans<= 0)  
44        ans+=N;  
45     return ans;  
46 }  
47 int main(){
48     scanf("%d%d",&N,&G);
49     if (G==P) {
50         printf("0");
51         return 0;
52     } 
53     G%=P;
54     for (int i=0;i<4;i++){//这里求阶乘 
55         fac[i][0]=1;
56         for (int j=1;j<=t[i];j++) fac[i][j]=(fac[i][j-1]*j)%t[i];
57     }
58     for (int i=1;i*i<=N;i++){
59         if(N%i==0){
60             int tmp=N/i;
61             for (int j=0;j<4;j++){
62                 if (tmp!=i) M[j]=(M[j]+lucas(N,i,j))%t[j];
63                 M[j]=(M[j]+lucas(N,tmp,j))%t[j];
64             }
65         }
66     }
67     printf("%d",pow(G,CRT(),P));
68 }