bzoj 1951: [Sdoi2010]古代猪文

 1 #include<cstdio>
 2 #include<iostream>
 3 #include<cstring>
 4 #include<cmath>
 5 #define ll long long
 6 #define mul 999911659
 7 using namespace std;
 8 int n,g,a[4];
 9 int sh[4]={2,3,4679,35617},C1[4][35618];
10 void exgcd(int a1,int a2,int &x,int &y)
11 {
12     if(!a2)
13       {
14         x=1;
15         y=0;
16         return;
17       }
18     exgcd(a2,a1%a2,x,y);
19     int t=x;
20     x=y;
21     y=t-a1/a2*y;
22 }
23 int kuai(ll n,int k,int p)
24 {
25     int ans=1;
26     for(;k;)
27       {
28         if(k%2)
29           ans=(ans*n)%p;
30         n=(n*n)%p;
31         k/=2;
32       }
33     return ans;
34 }
35 int C(int n,int m,int p)
36 {
37     if(n<m)
38       return 0;
39     return C1[p][n]*kuai(C1[p][m]*C1[p][n-m],sh[p]-2,sh[p])%sh[p];
40 }
41 int lucas(int n,int m,int p)
42 {
43     if(!m)
44       return 1;
45     return (C(n%sh[p],m%sh[p],p)*lucas(n/sh[p],m/sh[p],p))%sh[p];
46 }
47 int solve()
48 {
49     int x,y,a1,b1,a2,b2;
50     a1=sh[0];
51     b1=a[0];
52     for(int i=1;i<4;i++)
53       {
54         a2=sh[i];
55         b2=a[i];
56         exgcd(a1,a2,x,y);
57         x=((b2-b1)*x%a2+a2)%a2;
58         b1=b1+x*a1;
59         a1=a1*a2;
60       }
61     return b1;
62 }
63 int main()
64 {
65     scanf("%d%d",&n,&g);
66     if(g==mul)
67       {
68         printf("0\n");
69         return 0;
70       }
71     for(int i=0;i<4;i++)
72       {
73         C1[i][0]=1;
74         for(int j=1;j<=sh[i];j++)
75           C1[i][j]=(C1[i][j-1]*j)%sh[i];
76       }
77     for(int i=1;i<=sqrt(n);i++)
78       if(n%i==0)
79         {
80             for(int j=0;j<4;j++)
81               {
82                 a[j]=(a[j]+lucas(n,i,j))%sh[j];
83                 if(i!=n/i)
84                   a[j]=(a[j]+lucas(n,n/i,j))%sh[j];
85               }
86         }
87     printf("%d\n",kuai(g,solve(),mul));
88     return 0;
89 }