POJ1845 Sumdiv【快速模幂+素因子分解+等比数列+二分法】

/* POJ1845 Sumdiv */

#include <iostream>

using namespace std;

const int MOD = 9901;
const int N = 1000;

int fact[N+1], e[N+1], fcount;

void setfact(int n)
{
fcount = 0;
if(n && n % 2 == 0) {
fact[fcount] = 2;
e[fcount] = 0;
while(n && n % 2 == 0) {
e[fcount]++;
n /= 2;
}
fcount++;
}

for(int i=3; i*i <=n; i+=2) {
if(n % i == 0) {
fact[fcount] = i;
e[fcount] = 0;
while(n && n % i == 0) {
e[fcount]++;
n /= i;
}
fcount++;
}
}

if(n != 1) {
fact[fcount] = n;
e[fcount] = 1;
fcount++;
}
}

// 快速模幂计算函数
long long powermod(long long a, long long n, int m)
{
long long res = 1;
while(n) {
if(n & 1) {        // n % 2 == 1
res *= a;
res %= m;
}
a *= a;
a %= m;
n >>= 1;
}
return res;
}

//递归二分求 (1 + p + p^2 + p^3 +...+ p^n)%mod
long long sum(long long p, long long n)
{
if(n==0)
return 1;
else if(n % 2)
// 奇数：(1 + p + p^2 +...+ p^(n/2)) * (1 + p^(n/2+1))
return (sum(p, n / 2) * (1 + powermod(p, n/2+1, MOD))) % MOD;
else
// 偶数：(1 + p + p^2 +...+ p^(n/2-1)) * (1+p^(n/2+1)) + p^(n/2)
return (sum(p, n / 2 - 1) * (1 + powermod(p, n / 2 + 1, MOD)) + powermod(p, n / 2, MOD)) % MOD;
}

int main()
{
int a, b, ans;
while(cin >> a >> b) {
setfact(a);

ans = 1;
for(int i=0; i<fcount; i++)
ans = (ans * (sum(fact[i], e[i] * b ) % MOD)) % MOD;

cout << ans << endl;
}

return 0;
}