[SDOI2010]古代猪文
题目背景
“在那山的那边海的那边有一群小肥猪。他们活泼又聪明，他们调皮又灵敏。他们自由自在生活在那绿色的大草坪，他们善良勇敢相互都关心……”
——选自猪王国民歌
很久很久以前，在山的那边海的那边的某片风水宝地曾经存在过一个猪王国。猪王国地理位置偏僻，实施的是适应当时社会的自给自足的庄园经济，很少与外界联系，商贸活动就更少了。因此也很少有其他动物知道这样一个王国。
猪王国虽然不大，但是土地肥沃，屋舍俨然。如果一定要拿什么与之相比的话，那就只能是东晋陶渊明笔下的大家想象中的桃花源了。猪王勤政爱民，猪民安 居乐业，邻里和睦相处，国家秩序井然，经济欣欣向荣，社会和谐稳定。和谐的社会带给猪民们对工作火红的热情和对未来的粉色的憧憬。
小猪iPig是猪王国的一个很普通的公民。小猪今年10岁了，在大肥猪学校上小学三年级。和大多数猪一样，他不是很聪明，因此经常遇到很多或者稀奇 古怪或者旁人看来轻而易举的事情令他大伤脑筋。小猪后来参加了全猪信息学奥林匹克竞赛(Pig Olympiad in Informatics, POI)，取得了不错的名次，最终保送进入了猪王国大学(Pig Kingdom University, PKU)深造。
现在的小猪已经能用计算机解决简单的问题了，比如能用P++语言编写程序计算出A + B的值。这个“成就”已经成为了他津津乐道的话题。当然，不明真相的同学们也开始对他刮目相看啦~
小猪的故事就将从此展开，伴随大家两天时间，希望大家能够喜欢小猪。
题目描述
猪王国的文明源远流长，博大精深。
iPig在大肥猪学校图书馆中查阅资料，得知远古时期猪文文字总个数为N。当然，一种语言如果字数很多，字典也相应会很大。当时的猪王国国王考虑到 如果修一本字典，规模有可能远远超过康熙字典，花费的猪力、物力将难以估量。故考虑再三没有进行这一项劳猪伤财之举。当然，猪王国的文字后来随着历史变迁 逐渐进行了简化，去掉了一些不常用的字。
iPig打算研究古时某个朝代的猪文文字。根据相关文献记载，那个朝代流传的猪文文字恰好为远古时期的k分之一，其中k是N的一个正约数（可以是1和N）。不过具体是哪k分之一，以及k是多少，由于历史过于久远，已经无从考证了。
iPig觉得只要符合文献，每一种能整除N的k都是有可能的。他打算考虑到所有可能的k。显然当k等于某个定值时，该朝的猪文文字个数为N / k。然而从N个文字中保留下N / k个的情况也是相当多的。iPig预计，如果所有可能的k的所有情况数加起来为P的话，那么他研究古代文字的代价将会是G的P次方。
现在他想知道猪王国研究古代文字的代价是多少。由于iPig觉得这个数字可能是天文数字，所以你只需要告诉他答案除以999911659的余数就可以了。
输入输出格式
输入格式：输入文件ancient.in有且仅有一行：两个数N、G，用一个空格分开。
输出格式：输出文件ancient.out有且仅有一行：一个数，表示答案除以999911659的余数。
输入输出样例
4 2
2048
说明
数据规模
10%的数据中，1 <= N <= 50；
20%的数据中，1 <= N <= 1000；
40%的数据中，1 <= N <= 100000；
100%的数据中，1 <= G <= 1000000000，1 <= N <= 1000000000。
题目大意：给定N,G，求
但是的x会很大，以至于无法用存储
这时要用到指数取模算法
对于G^x，我们有
下面给出证明
因为Gμ(p)≡1 (mod p)
我们令x=k*μ(p)+b
那么Gx=(Gμ(p))^k*Gb
因为(Gμ(p))^k≡1 (mod p)
因为b=x%μ(p)=x%(p-1)
所以Gx≡Gx%(p-1) (mod p)
所以对指数取模，模p-1就行了
组合数取模用Lucas
这时有一个问题：p-1不是素数，不能用lucas和求逆元
所以将p-1分解为4个素数，分别算出同余方程，再用中国剩余定理合并
这里写出p-1的4个素数：
999911658=2*3*4679*35617
1 #include<iostream> 2 #include<cstdio> 3 #include<cstring> 4 #include<algorithm> 5 using namespace std; 6 typedef long long ll; 7 ll Mod=999911658; 8 ll fac[100001],A[100001],cnt,ys[100001],n,g,pri[5],b[5],c[5],a[5]; 9 ll exgcd(ll a,ll b,ll &x,ll &y) 10 { 11 if (b==0) 12 { 13 x=1;y=0; 14 return a; 15 } 16 ll d=exgcd(b,a%b,x,y); 17 ll t=x;x=y;y=t-(a/b)*y; 18 return d; 19 } 20 ll rev(ll a,ll p) 21 { 22 ll x,y; 23 exgcd(a,p,x,y); 24 return (x%p+p)%p; 25 } 26 ll C(ll x,ll y,ll p) 27 { 28 if (x>y) return 0; 29 if (x<p&&y<p) return (fac[y]*A[x]%p)*A[y-x]%p; 30 return C(x%p,y%p,p)*C(x/p,y/p,p)%p; 31 } 32 ll qpow(ll x,ll y,ll p) 33 { 34 ll res=1; 35 while (y) 36 { 37 if (y&1) res=res*x%p; 38 x=x*x%p; 39 y/=2; 40 } 41 return res; 42 } 43 ll cal(int p) 44 {int i; 45 fac[0]=1; 46 for (i=1;i<p;i++) 47 fac[i]=fac[i-1]*i%p; 48 A[1]=1;A[0]=1; 49 for (i=2;i<p;i++) 50 A[i]=(p-p/i)*A[p%i]%p; 51 for (i=1;i<p;i++) 52 A[i]=A[i]*A[i-1]%p; 53 ll s=0; 54 for (i=1;i<=cnt;i++) 55 { 56 s+=C(ys[i],n,p); 57 s%=p; 58 } 59 return s; 60 } 61 int main() 62 {int i; 63 cin>>n>>g; 64 if (g==Mod+1) 65 { 66 cout<<0; 67 return 0; 68 } 69 for (i=1;i*i<=n;i++) 70 if (n%i==0) 71 { 72 if (i*i==n) 73 ys[++cnt]=i; 74 else 75 { 76 ys[++cnt]=i; 77 ys[++cnt]=n/i; 78 } 79 } 80 pri[1]=2;pri[2]=3;pri[3]=4679;pri[4]=35617; 81 for (i=1;i<=4;i++) 82 b[i]=cal(pri[i]); 83 for (i=1;i<=4;i++) 84 c[i]=Mod/pri[i]; 85 for (i=1;i<=4;i++) 86 a[i]=rev(c[i],pri[i]); 87 ll s=0; 88 for (i=1;i<=4;i++) 89 { 90 s+=((c[i]*a[i]%Mod)*b[i])%Mod; 91 s%=Mod; 92 } 93 printf("%lld\n",qpow(g,s,Mod+1)); 94 }
