POJ 1845 Sumdiv (快速幂)

题目链接 ： http://poj.org/problem?id=1845

题意 ： 求∑(x^n的所有因子) mod 9901;

思路 ： 如果 x = p1^a1 * p2^a2 * p3^a3...pn^an， (p1, p2, ...都是x的质因子)， 那么所求解S = (1 + p1 + p1^2 ...+p1^(a1*n)) * (1 + p2 +....).....

所以其实只要解决一个等比数列求和问题就行了， 坑的是这题给的mod 有点小， 如果用公式 (x ^(n+1) - 1) / (x - 1)去做得话， 需要求逆元， 但是存在逆元需要满足

gcd(a, b) == 1

这一性质。  我是用快速幂来做的。 



CODE ：

#include <cstdio>
#include <cstring>
#include <algorithm>

using namespace std;
typedef __int64 lld;
const int mod = 9901;
struct mat{
lld a[4][4];
};

void init(mat & a){
for (int i = 1; i <= 3; i++)
for (int j = 1; j <= 3; j++)
a.a[i][j] = 0;
}

mat mul(mat a, mat b){
mat ans;
for (int i = 1; i <= 3; i++){
for (int j = 1; j <= 3; j++){
ans.a[i][j] = 0;
for (int k = 1; k <= 3; k++)
ans.a[i][j] = (ans.a[i][j] + (a.a[i][k] % mod) * (b.a[k][j] % mod)) % mod;
}
}
return ans;
}

void debug(mat & a){
for (int i = 1; i <= 3; i++){
for (int j = 1; j <= 3; j++)printf("%d ", a.a[i][j]);
printf("\n");
}
}

lld calc(lld x, lld n){
mat ans, e;
init(e); init(ans);
for (int i = 1; i <= 3; i++)ans.a[i][i] = 1;
e.a[1][1] = 1; e.a[2][2] = x; e.a[3][1] = 1; e.a[3][3] = x;
while (n){
if (n & 1){
ans = mul(ans, e);
}
e = mul(e, e);
n >>= 1;
}
//debug(ans);
return (ans.a[3][1] + ans.a[3][3]) % mod;
}
const int maxn = 10005;
int cnt[maxn], p[maxn];

lld solve(lld x, lld n){
memset(p, 0, sizeof(p));
memset(cnt, 0, sizeof(cnt));
int m = 0;
for (int i = 2; i <= mod && x > 1; i++){
if (x % i == 0 && x > 1){
p[++m] = i;
while (x % i == 0 && x > 1){
cnt[m]++;
x /= i;
}
}
}
if (x > 1){
m++;
p[m] = x;
cnt[m] = 1;
}
lld res = 1;
for (int i = 1; i <= m; i++){
res = (res * calc(p[i], cnt[i] * n)) % mod;
}
return res;
}

lld x, n;

int main(){
while(scanf("%I64d%I64d", &x, &n) != EOF){
printf("%I64d\n", solve(x, n));
}
return 0;
}