[POJ1845]Sumdiv

 1 #include <iostream>
 2 #include <cstdio>
 3 #include <cstring>
 4 #include <vector>
 5 #include <algorithm>
 6 #include <cmath>
 7 
 8 using namespace std;
 9 
10 typedef long long LL;
11 const int mod = 9901;
12 vector<LL> p;
13 vector<LL> k;
14 LL n, a, b;
15 
16 //快速幂
17 LL quickmul(LL x, LL n) {
18     LL ans = 1;
19     LL t = x;
20     while(n) {
21         if(n & 1) {
22             ans = (ans * t) % mod;
23         }
24         t = t * t % mod;
25         n >>= 1;
26     }
27     return ans;
28 }
29 
30 //分解质因数
31 void factor(LL n) {
32     p.clear();
33     k.clear();
34     int nn = (int)sqrt(n*1.0);
35     for(int i = 2; i <= nn; i+=2) {
36         if(n % i == 0) {    //分解n,使n变为p1^k1*p2^k2+...pn^kn
37             p.push_back(i);
38             k.push_back(0);
39             while(n % i == 0) {
40                 k.back()++;
41                 n = n / i;
42             }
43         }
44         if(i == 2) {
45             i--;
46         }
47         if(n == 1){
48             break;
49         }
50     }
51     if(n != 1) {    //特判质数
52         p.push_back(n);
53         k.push_back(1);
54     }
55 }
56 
57 //同余模公式
58 //(a+b)%m=(a%m+b%m)%m
59 //(a*b)%m=(a%m*b%m)%m
60 //二分求1+pi+pi^2+...+pi^n的和
61 LL modu(LL x,LL y) {  
62     if(y == 0) {
63         return 1;  
64     }
65     if(y % 2 == 0) {
66         return (((modu(x,y/2-1)%mod)*((1+quickmul(x,y/2+1))%mod))%mod+quickmul(x,y/2)%mod)%mod;  
67     }
68     if(y % 2 != 0) {
69         return (modu(x,y/2)%mod)*((1+quickmul(x,y/2+1))%mod)%mod;  
70     }
71 }
72 
73 void solve() {
74     factor(a);  //分解a,使a = p1^k1*p2^k2+...pn^kn
75     for(int i = 0; i < k.size(); i++) {
76         k[i] = k[i] * b; //所以a^b = p1^(k1*b)*p2^(k2*b)...pn^(kn*b)
77     }
78     LL sum = 1;
79     for(int i = 0; i < k.size(); i++) {
80         sum = sum * modu(p[i], k[i]) % mod;
81     }
82     printf("%I64d\n", sum);
83 }
84 
85 int main() {
86     while(~scanf("%d %d", &a, &b)) {
87         solve();
88     }
89     return 0;
90 }