题目链接：
POJ 1845 Sumdiv
题意：
给出A和B,求A^B的所有因子和对990取余后的值。0<= A, B <= 50000000.
分析：
特判A=1,A=0,B=0的情况，根据任何一个大于1的自然数若其不为素数则必可唯一的分解为有限个素数的乘积。
假设A = p1^a1 * p2^a2 * p3^a3 * ... * pn^an，其中(p1,p2,p3,...,pn均为素数，a1,a2,...,an是拆分后的指数)，那么A^B = p1^(a1*B) * p2^(a2*B) * p3^(a3*B) * ... * pn^(an*B);.
那么所有因子和的表达式为sum = (1 + p1 + p1^2 + ... + p1^(a1*B)) * (1 + p2 + p2^2 + ... + p2^(a2*B)) *...*(1 + pn + pn^2 + pn^(an*B)).
对于sum的每一项可以二分求等比数列之和。 也可以利用等比数列求和公式，但是需要用的快速乘法。

在使用等比数列求和公式时，慎用乘法逆元。用求和公式可得：(pi^(ei + 1) - 1) / (pi - 1) % mod. 但是(pi - 1）可能是mod的倍数，这样就不能直接使用费马小定理了。因为求分母的逆元要求分母*分母^(-1)%mod=1,但分母%mod=0恒成立，找不出模mod下的逆元。需要特判pi % mod == 0 （直接跳过）,pi % mod == 1（ans = ans * (num * B + 1) % mod）.然后用ex_gcd求逆元，还要考虑求出来的结果为负的情况，详细代码看最后一个。

#include <iostream>
#include <cstdio>
#include <cstring>
#include <string>
#include <algorithm>
#include <climits>
#include <cmath>
#include <ctime>
#include <cassert>
#include <bitset>
#define IOS ios_base::sync_with_stdio(0); cin.tie(0);
using namespace std;
typedef long long ll;
const ll mod = 9901;
const int MAX_N = 10010;

int cnt;
ll A, B;
int prime[MAX_N], vis[MAX_N];

void GetPrime()
{
    cnt = 0;
    memset(vis, 0, sizeof(vis));
    for(int i = 2; i < MAX_N; i++){
        if(vis[i] == 0){
            prime[cnt++] = i;
            for(int j = i * i; j < MAX_N; j += i){
                vis[j] = 1;
            }
        }
    }
}

ll quick_pow(ll n, ll m)
{
    ll res = 1, tmp = n % mod;
    while(m){
        if(m & 1) {
            res = res * tmp % mod;
            m--;
        }
        tmp = tmp * tmp % mod;
        m >>= 1;
    }
    return res;
}

ll sum(ll a, ll n) //计算a^0 + a^1 + a^2 + ... + a^n
{
   //n分奇偶考虑，例如n = 3 和n = 4的情况。
    if(n == 0) return 1;
    ll ans = 1, r, base, extra;
    base = sum(a, (n - 1) / 2); 
    r = quick_pow(a, n / 2 + 1);
    ans = (base + base * r % mod) % mod;
    if(n % 2 == 0){
        ans = (ans + quick_pow(a, n / 2)) % mod;
    }
    //printf("a = %lld n = %lld base = %lld r = %lld\n", a, n, base, r);
    return ans;

}

void solve()
{
    ll ans = 1;
    for(int i = 0; prime[i] * prime[i] <= A; i++){
        if(A % prime[i] == 0){
            int num = 0;
            while(A % prime[i] == 0){
                num ++;
                A /= prime[i];
            }
            //printf("prime[i] = %d num = %d\n", prime[i], num);
            ans = ans * sum((ll)prime[i], (ll)num * B) % mod;
        }
    }
    if(A > 1) ans = ans * sum(A, B) % mod; //别忘了这个！最后A可能是一个特别大的素数！
    printf("%lld\n", ans);
}

int main()
{
    GetPrime();
    while(~scanf("%lld%lld", &A, &B)){
        if(A == 1 || A == 0 || B == 0) printf("1\n");
        else solve();
    }
    return 0;
}

#include <iostream>
#include <cstdio>
#include <cstring>
#include <string>
#include <algorithm>
#include <climits>
#include <cmath>
#include <ctime>
#include <cassert>
#define IOS ios_base::sync_with_stdio(0); cin.tie(0);
using namespace std;
typedef long long ll;
const ll MOD = 9901;
const int MAX_N = 10010;

ll  mod, A, B;
int cnt, prime[MAX_N], vis[MAX_N]; 

void GetPrime()
{
    cnt = 0;
    memset(vis, 0, sizeof(vis));
    for(int i = 2; i < MAX_N; i++){
        if(vis[i] == 0){
            prime[cnt++] = i;
            for(int j = i * i; j < MAX_N; j += i){
                vis[j] = 1;
            }
        }
    }
}

ll mul(ll n, ll m)
{
    ll res = 0, tmp = n % mod;
    while(m){
        if(m & 1) res = (res + tmp) % mod;
        tmp = (tmp + tmp) % mod;
        m >>= 1; 
    }
    return res;
}

ll quick_pow(ll n, ll m)
{
    ll res = 1, tmp = n % mod;
    while(m){
        if(m & 1) res = mul(res, tmp);
        tmp = mul(tmp, tmp);
        m >>= 1;
    }
    return res;
}

void solve()
{
    ll ans = 1;
    for(int i = 0; prime[i] * prime[i] <= A; i++){
        if(A % prime[i] == 0){
            int num = 0;
            while(A % prime[i] == 0){
                num++;
                A /= prime[i];
            }
            mod = (prime[i] - 1) * MOD;
            ans = ans * (quick_pow(prime[i], num * B + 1) + mod - 1) / (prime[i] - 1) % MOD;
        }
    }
    if(A > 1){
        mod = (A - 1) * MOD;
        ans = ans * (quick_pow(A, B + 1) + mod - 1) / (A - 1) % MOD;
    }
    printf("%lld\n", ans);
}

int main()
{
    GetPrime();
    while(~scanf("%lld%lld", &A, &B)){
        if(A == 1 || A == 0 || B == 0) printf("1\n");
        else solve();
    }
    return 0;
}

#include <iostream>
#include <cstdio>
#include <cstring>
#include <string>
#include <algorithm>
#include <climits>
#include <cmath>
#include <ctime>
#include <cassert>
#include <bitset>
#define IOS ios_base::sync_with_stdio(0); cin.tie(0);
using namespace std;
typedef long long ll;
const ll mod = 9901;
const int MAX_N = 500000;

ll A, B;
int prime_cnt;
int prime[MAX_N];
bitset<MAX_N> bs;

void GetPrime()
{
    bs.set();
    prime_cnt = 0;
    for(int i = 2; i < MAX_N; i++){
        if(bs[i] == 1) prime[prime_cnt++] = i;
        for(int j = 0; j < prime_cnt && i * prime[j] < MAX_N; j++){
            bs[i * prime[j]] = 0;
            if(i % prime[j] == 0) break;
        }
    }
    //printf("bs[9901] = %d\n", (int)bs[9901]);
}

ll mul(ll a, ll b, ll m)
{
    ll res = 0, tmp = a % m;
    b %= m;
    while(b){
        if(b & 1) res = (res + tmp) % m;
        tmp = (tmp + tmp) % m;
        b >>= 1;    
    }
    return res;
}

ll quick_pow(ll a, ll b, ll m)
{
    ll res = 1, tmp = a % m;
    while(b){
        if(b & 1) res = mul(res, tmp, m);
        tmp = mul(tmp, tmp, m);
        b >>= 1;
    }
    return res;
}

ll ex_gcd(ll a, ll b, ll& x, ll& y)
{
    if(b == 0) {
        x = 1, y = 0;
        return a;
    }
    ll d = ex_gcd(b, a % b, y, x);
    y -= a / b * x;
    return d;
}

void solve()
{
    ll x = A, ans = 1;
    for(int i = 0; i < prime_cnt && prime[i] * prime[i] <= x; i++){
        if(x % prime[i] == 0){
            int num = 0;
            while(x % prime[i] == 0){
                num++;
                x /= prime[i];
            }
            if(prime[i] % mod == 0) continue;
            else if(prime[i] % mod == 1){
                ans = ans * (num * B + 1) % mod;
                continue;
            }
            ans = ans * (quick_pow(prime[i], num * B + 1, mod) + mod - 1) % mod;
            ll xx, y;
            ex_gcd((ll)prime[i] - 1, mod, xx, y);
            xx = (xx % mod + mod) % mod;
            ans = ans * xx % mod;
        }
    }
    if(x > 1){
        if(x % mod == 0) {
            printf("%lld\n", ans);
            return;
        }
        else if(x % mod == 1) {
            ans = ans * (B + 1) % mod;
            printf("%lld\n", ans);
            return;
        }
        ans = ans * (quick_pow(x, B + 1, mod) + mod - 1) % mod;
        ll xx , y;
        ex_gcd(x - 1, mod, xx, y);
        xx = (xx % mod + mod) % mod;
        ans = ans * xx % mod;
    }
    printf("%lld\n", ans);
}

int main()
{
    GetPrime();
    while(~scanf("%lld%lld", &A, &B)){
        if(A == 0 || A == 1) printf("%lld\n", A);
        else if(B == 0) printf("1\n");
        else solve();
    }
    return 0;
}