POJ1845Sumdiv（求所有因子和 + 唯一分解定理）

Sumdiv

Description

Consider two natural numbers A and B. Let S be the sum of all natural divisors of A^B. Determine S modulo 9901 (the rest of the division of S by 9901).

Input

The only line contains the two natural numbers A and B, (0 <= A,B <= 50000000)separated by blanks.

Output

The only line of the output will contain S modulo 9901.

Sample Input

2 3

Sample Output

15

Hint

2^3 = 8.
The natural divisors of 8 are: 1,2,4,8.
Their sum is 15.
15 modulo 9901 is 15 (that should be output).

Source

Romania OI
2002

题意：求A ^ B的所有因子的和；

分析：对A用唯一分解定理分解 A = p1 ^ a1 * p2 ^ a2 * p3 ^ a3 ... * pn ^ an

其中A的因子的个数为 （ 1 + a1) * ( 1 + a2 ) * ( 1 + a3 ) * ... * ( 1 + an)

则A的所有因子的和为 （1 + p1 + p1 ^ 2 + p1 ^ 3 ... + p1 ^ a1) * ( 1 + p2 ^ 1 + p2 ^ 2 + p3 ^ 3 + ... + p2 ^ a2 ) * ... * ( 1 + pn + pn ^ 2 + ... + pn ^ an)

求  a + a ^ 2 + a ^ 3 + ... + a ^ n

如果n为奇数：  a + a ^ 2 + ... + a ^ (n / 2 )  + a ^ (n / 2 + 1) + ( a + a ^ 2 + ... + a ^ ( n / 2 ) ) *  a ^ ( n / 2 + 1）

如果n为偶数：  a + a ^ 2 + ... + a ^ (n / 2 ） + （ a + a ^ 2 + ... + a ^ ( n /. 2) ) * a ^ (n / 2)

 1 #include <iostream>
 2 #include <cstring>
 3 #include <cstdio>
 4 #include <algorithm>
 5 using namespace std;
 6 typedef long long LL;
 7 const int Max = 10000;
 8 const int Mod = 9901;
 9 int prime[Max + 10],flag[Max],cnt;
10 void get_prime()
11 {
12     cnt = 0;
13     memset(flag, 0, sizeof(flag));
14     for(int i = 2; i <= Max; i++)
15     {
16         if(flag[i] == 0)
17         {
18             flag[i] = 1;
19             prime[++cnt] = i;
20             for(int j = i; j <= Max / i; j++)
21                 flag[i * j] = 1;
22         }
23     }
24 }
25 LL pow_mod(LL n, LL k)
26 {
27     LL res = 1;
28     while(k)
29     {
30         if(k & 1)
31             res = res * n % Mod;
32         n = n * n % Mod;
33         k >>= 1;
34     }
35     return res;
36 }
37 LL get_sum(LL n, LL m)
38 {
39     if(m == 0)
40         return 1;
41     if(m & 1)
42     {
43         return get_sum(n, m / 2) *( 1 + pow_mod(n, m / 2 + 1) ) % Mod;
44     }
45     else
46     {
47         return ( get_sum(n, m / 2 - 1) * (1 + pow_mod(n, m / 2 + 1)) % Mod + pow_mod(n, m / 2) ) % Mod;
48     }
49 }
50 int main()
51 {
52     LL a,b;
53     get_prime();
54     while(scanf("%I64d%I64d", &a, &b) != EOF)
55     {
56         LL ans = 1;
57         if(a == 0 && b) //特殊情况
58             ans = 0;
59         LL m;
60         for(int i = 1; i <= cnt; i++)
61         {
62             if(prime[i] > a)
63                 break;
64             m = 0;
65             if(a % prime[i] == 0)
66             {
67                 while(a % prime[i] == 0)
68                 {
69                     a = a / prime[i];
70                     m++;
71                 }
72                 m = m * b;  // m要设成LL，否则这里会溢出
73                 ans = ans * get_sum((LL)prime[i], m) % Mod;
74             }
75         }
76         if(a > 1)
77             ans = ans * get_sum(a, b) % Mod;
78         printf("%I64d\n", ans);
79     }
80     return 0;
81 }