題目鏈接：http://poj.org/problem?id=1845

思路：

１．整數唯一分解定理：

任意正整數都有且只有一種方式寫出其素因子的乘積表達式。

a=(p1^k1)*(p2^k2)*(p3^k3)*....*(pn^kn) 其中pi均為素數

２．約數和公式：

　　對於已經分解的整數a=(p1^k1)*(p2^k2)*(p3^k3)*....*(pn^kn)

　　有a的所有因子之和為

　S` = (1+p1+p1^2+p1^3+...p1^k1) * (1+p2+p2^2+p2^3+….p2^k2) * (1+p3+ p3^3+…+ p3^k3) * .... * (1+pn+pn^2+pn^3+...pn^kn)

那麽 a^b 的所有因子和為

　　S = (1+p1+p1^2+p1^3+...p1^(k1*b)) * (1+p2+p2^2+p2^3+….p2^(k2*b)) * (1+p3+ p3^3+…+ p3^(k3*b)) * .... * (1+pn+pn^2+pn^3+...pn^(kn*b))

對於數列 1, p, p^2, p^3 ... p^n % mod，其中 mod 為質數，打個表可以發現該數列是一個循環數列，其中存在一個循環節為 1, p, p^1, ... p^(mod-2)．其實這點在費馬小定理中是有體現的，a^(mod-1) = 1 (% mod)．

那麽對於求 cnt = 1 + p + p^2 + ... + p^n % mod，可以令

　　cc1 = (n + 1) / (mod - 1)

　　cc2 = (n + 1) % (mod - 1)

　　cnt1 = 1 + p + p^2 + ... + p^(mod - 2)

　　cnt2 = 1 + p + p^2 + ... + p^(cc2 - 1)

那麽 cnt = cc1 * cnt1 + cnt2 % mod

對於求 a^b，可以先將 a 質因分解，得到 S 的表達式，對於 S 表達式，只需要按照上面的方法求出其中每個乘數項即可．時間復雜度為 O(loga * mod), 本題中 mod = 9901, 時間上是允許的．

代碼：

 1 #include <iostream>
 2
 
 #include <stdio.h>
 3 #include <map>
 4 #define ll long long
 5 using namespace std;
 6 
 7 const int mod = 9901;
 8 const int MAXN = 1e5 + 10;
 9 ll prime[MAXN], indx = 0;
10 map<int, int> num;
11 
12 int get(ll a, ll b){//計算sigma(a^i),其中0<=i<=b
13     ll sol1 = 0, sol2 = -1, cnt = 1;
14     ll cc1 = (b + 1) / (mod - 1);
15     ll cc2 = (b + 1) % (mod - 1);
16     for(int i = 0; i < mod - 1; i++){
17         sol1 = (sol1 + cnt) % mod;
18         if(sol2 == -1 && i + 1 == cc2) sol2 = sol1;
19         cnt = (cnt * a) % mod;
20     }
21     return (sol1 * cc1 % mod + sol2) % mod;
22 }
23 
24 int main(void){
25     ll a, b, sol = 1;
26     cin >> a >> b;
27     for(int i = 2; i * i <= a; i++){
28         if(a % i == 0){
29             prime[indx] = i;
30             while(a % i == 0){
31                 num[i]++;
32                 a /= i;
33             }
34             indx++;
35         }
36     }
37     if(a > 1) prime[indx++] = a, num[a]++;
38     for(int i = 0; i < indx; i++){
39         sol = (sol * get(prime[i], b * num[prime[i]])) % mod;
40     }
41     cout << sol << endl;
42     return 0;
43 }

但是當 mod 比較大時這個方法就不行了，mod 比較大時可以用下面這個代碼，貼的 Kuangbin 模板

代碼：

 1 #include <stdio.h>
 2 #include <math.h>
 3 #include <iostream>
 4 #include <algorithm>
 5 #include <string.h>
 6 #define ll long long
 7 using namespace std;
 8 
 9 //******************************************
10 //素數篩選和合數分解
11 const int MOD = 9901;
12 const int MAXN=10000;
13 int prime[MAXN+1];
14 
15 void getPrime(void){
16     memset(prime, 0, sizeof(prime));
17     for(int i = 2; i <= MAXN; i++)
18     {
19         if(!prime[i]) prime[++prime[0]] = i;
20         for(int j=1; j<= prime[0]&&prime[j] <= MAXN/i; j++)
21         {
22             prime[prime[j]*i] = 1;
23             if(i % prime[j] == 0) break;
24         }
25     }
26 }
27 
28 ll factor[100][2];
29 int fatCnt;
30 
31 int getFactors(ll x){
32     fatCnt = 0;
33     ll tmp = x;
34     for(int i=1; prime[i] <= tmp/prime[i]; i++){
35         factor[fatCnt][1] = 0;
36         if(tmp % prime[i] == 0){
37             factor[fatCnt][0] = prime[i];
38             while(tmp % prime[i] == 0){
39                 factor[fatCnt][1]++;
40                 tmp /= prime[i];
41             }
42             fatCnt++;
43         }
44     }
45     if(tmp!=1)
46     {
47         factor[fatCnt][0]=tmp;
48         factor[fatCnt++][1]=1;
49     }
50     return fatCnt;
51 }
52 
53 //******************************************
54 ll pow_m(ll a, ll n)//快速模冪運算
55 {
56     ll res = 1;
57     ll tmp = a % MOD;
58     while(n){
59         if(n & 1){ 
60             res *= tmp;
61             res%=MOD;
62         }
63         n >>= 1;
64         tmp *= tmp;
65         tmp %= MOD;
66     }
67     return res;
68 }
69 
70 ll sum(ll p, ll n){//計算1+p+p^2+...+p^n
71     if(p == 0)return 0;
72     if(n == 0)return 1;
73     if(n & 1){//奇數
74         return ((1 + pow_m(p, n/2 + 1)) % MOD * sum(p, n / 2) % MOD) % MOD;
75     }else return ((1 + pow_m(p, n / 2 + 1)) % MOD * sum(p, n / 2 - 1) + pow_m(p, n / 2) % MOD) % MOD;
76 
77 }
78 
79 int main(void){
80     int A, B;
81     getPrime();
82     while(scanf("%d%d", &A, &B) != EOF){
83         getFactors(A);
84         ll ans = 1;
85         for(int i = 0; i < fatCnt; i++){
86             ans *= (sum(factor[i][0], B * factor[i][1]) % MOD);
87             ans %= MOD;
88         }
89         printf("%lld\n",ans);
90     }
91     return 0;
92 }

hdu1845(a^b的因子和%p)