学步园

题目链接 ： 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;
}