Sumdiv POJ - 1845

题目传送门

题意：给你一个A，B，问你A^B的所有正因子的和模9901的结果。

思路：这一个题目学了好多东西啊，先学习了逆元，发现如果是出发的话要用逆元来取模
a / b % MOD = (a % (MOD * b)) / b
然后这个题目是可以用唯一分解定理来做的
N=P1^a1 * P2^a2 * P3a3 * …… * Pn^an
唯一分解定理有一个结论，一个数字的所有的正因子的和为
(1 + P1 + P1^2 + … + P1^a1) * (1 + P2 + P2^2 + … + P2^a2) * … *
(1 + Pn + Pn^2 + … + Pn^an)
然后就可以用等比数列来做了，但是又遇到了一个新的问题，好像只用快速幂过不了啊，用了一下二分乘法终于过了这个题目。

#include <algorithm>
#include <cmath>
#include <cstdio>
#include <cstring>
#include <fstream>
#include <iostream>
#include <list>
#include <map>
#include <queue>
#include <set>
#include <sstream>
#include <stack>
#include <string>
#include <vector>

#define MAXN 20100
#define MAXE 5
#define INF 1000000000
#define MOD 9901
#define LL long long
#define pi acos(-1.0)

using namespace std;

LL p[MAXN];
bool prime[MAXN];
int cnt;

struct Node {
    LL x;
    LL y;
};

vector<Node> vec;

void is_prime() {
    cnt = 0;
    for (int i = 0; i < MAXN; ++i) {
        prime[i] = true;
    }
    for (int i = 2; i < MAXN - 10; i++) {
        if (prime[i]) {
            p[cnt++] = i;
            for (int j = i + i; j < MAXN - 10; j += i)
                prime[j] = false;
        }
    }
}

LL multi(LL a, LL b, LL m) {
    LL ans = 0;
    a %= m;
    while (b) {
        if (b & 1) {
            ans = (ans + a) % m;
            b--;
        }
        b >>= 1;
        a = (a + a) % m;
    }
    return ans;
}

LL quick_mod(LL a, LL b, LL m) {
    LL ans = 1;
    a %= m;
    while (b) {
        if (b & 1) {
            ans = multi(ans, a, m);
            b--;
        }
        b >>= 1;
        a = multi(a, a, m);
    }
    return ans;
}

LL solve(LL a, LL b) {
    for (int i = 0; i < cnt; ++i) {
        LL ans = 0;
        while (a % p[i] == 0) {
            ans++;
            a /= p[i];
        }
        if (ans > 0) {
            vec.push_back((Node) {p[i], ans});
        }
    }
    LL sum = 1;
    for (int i = 0; i < vec.size(); ++i) {
        LL M = (vec[i].x - 1) * MOD;
        sum = sum * (quick_mod(vec[i].x, vec[i].y * b + 1, M) + M - 1) / (vec[i].x - 1);
        sum %= MOD;
    }
    if (a > 1) {
        LL M = MOD * (a - 1);
        sum *= (quick_mod(a, b + 1, M) + M - 1) / (a - 1);
        sum %= MOD;
    }
    return sum % MOD;
}

int main() {
    std::ios::sync_with_stdio(false);
    is_prime();
    LL a, b;
    cin >> a >> b;
    cout << solve(a, b) << endl;
    return 0;
}