POJ-1845 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). 

题目大意：

给你两个数a和b，让你求a的b次方的因子的和。

解题思路：

蒟蒻是从这篇博客(http://blog.csdn.net/lyy289065406/article/details/6648539)里面得到的思路0.0

这篇博客说的很详细，对于数论新手是非常友好的0.0

以及这道题似乎还可以通过逆元求解。不过我还是不会0.0

这道题的思路其实主要来自于两个定理：

1.对于任何整数，都有且只有一种方式能够写出其素因子的乘积形式。

2.对于已经分解的整数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)

根据这两个定理，就可以解出这个题目了。

代码：

#include <cstdio>
#include <cstring>
using namespace std;
#define mod 9901
typedef long long LL;

const int maxn = 1e4 + 5;
LL prime[maxn], f[maxn], g[maxn];

void init(){
    memset(prime, 0, sizeof(prime));
    for(int i = 2; i < maxn; ++i){
        if(!prime[i]) prime[++prime[0]] = i;
        for(int j = 1; j <= prime[0] && prime[j] < maxn / i; ++j){
            prime[prime[j] * i] = 1;
            if(i % prime[j] == 0) break;
        }
    }
}
LL mod_mul(LL a, LL b){
    LL res = 0;
    while(b){
        if(b & 1) res = (res + a) % mod;
        a = (a + a) % mod;
        b >>= 1;
    }
    return res % mod;
}
LL mod_exp(LL a, LL b){
    LL res = 1;
    while(b){
        if(b & 1) res = mod_mul(res, a);
        a = mod_mul(a, a);
        b >>= 1;
    }
    return res % mod;
}
LL Cal(LL x, LL num){
    if(num == 0) return 1;
    if(num == 1) return (1 + x) % mod;
    
    LL res = 1;
    if(num % 2 == 0){
        res = (res + mod_exp(x, num / 2 + 1)) % mod;
        res = (res * Cal(x, num / 2 - 1)) % mod;
        res = (res + mod_exp(x, num / 2)) % mod;
    }else{
        res = (res + mod_exp(x, num / 2 + 1)) % mod;
        res = (res * Cal(x, num / 2)) % mod;
    }
    return res;
}
int main(){
    init();
    LL a, b, tmp, cnt;
    while(~scanf("%lld%lld", &a, &b)){
        if(a == 1) {puts("1"); continue;}
        tmp = a; cnt = 0;
        // memset(f, 0, sizeof(f));
        // memset(g, 0, sizeof(g));
        for(int i = 1; i <= prime[0]; ++i){
            if(tmp % prime[i] == 0){
                f[cnt] = prime[i];
                g[cnt] = 0;
                while(tmp % prime[i] == 0){
                    tmp /= prime[i];
                    ++g[cnt];
                }
                ++cnt;
            }
            
            if(tmp == 1) break;
        }
        if(tmp != 1){
            f[cnt] = tmp; g[cnt++] = 1;
        }
        LL ans = 1;
        for(int i = 0; i < cnt; ++i) {
            g[i] *= b;
            // printf("f = %lld, g = %lld\n", f[i], g[i]);
            ans = (ans * Cal(f[i], g[i])) % mod;
        }
        printf("%lld\n", ans);
    }
    return 0;
}