\(Sumdiv\)

  题意: 求\(A^B\)的所有约数之和\(\mod{9901}\left(1\leq A,B \leq5*10^7\right)\)

  分析: A分解质因数为\(p_1^{c_1}\times p_2^{c_2}\times\cdots\times p_n^{c_n}\)。那么\(A^B\)表示为\(p_1^{Bc_1}\times p_2^{Bc_2}\times\cdots\times p_n^{Bc_n}\)。\(A^B\)的所有约数表示为集合\(\{p_1^{k_1}\times p_2^{k_2}\times\cdots\times p_n^{k_n}\}\)，其中\(0\leq k_i \leq B\times c_i\left(1\leq i \leq n\right)\)

  由乘法分配律, \(A^B\)的所有约数之和为:

\[\left(1+p_1+\cdots+p_1^{B*c_1}\right)*\left(1+p_2+\cdots+p_2^{B*c_2}\right)*\cdots\\*\left(1+p_n+\cdots+p_n^{B*c_n}\right)=\prod_{i=1}^n\left(\sum_{j=0}^{B*c_i}\left(p_i\right)^j\right)\]

  括号内为等比数列，直接使用求和公式，需要做除法。答案还需对9901取模，mod运算只对加、减、乘有分配律，不能对分子分母取模后做除法，换一种思路，用分治法对等比数列求和。

  用分治法求\(\displaystyle\sum_{i=1}^cp^i\quad\text{若c为奇数:}\) \[\sum_{i=1}^cp^i=\left(1+p+\cdots+p^\frac{c-1}{2}\right)+\left(p^\frac{c+1}{2}+\cdots+p^c\right)\\=\left(1+p+\cdots+p^\frac{c-1}{2}\right)+p^\frac{c+1}{2}*\left(1+p+\cdots+p^\frac{c-1}{2}\right)\\=(1+p^\frac{c+1}{2})*\sum_{i=1}^\frac{c+1}{2}p^i\]

  若c为偶数，类似有:

\[\sum_{i=1}^cp^i=(1+p^\frac{c}{2})*\sum_{i=1}^{\frac{c}{2}-1}p^i+p^c\]

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
#include <iostream>
#include <cstdio>
using namespace std;
#define ll long long
const ll N = 1005, MOD = 9901;
pair<ll, ll> fac[N];
ll cnt = 0;

ll qpow(ll a, ll b)
{
    ll res = 1;
    for (; b; b >>= 1, a = a * a % MOD)
        if (b & 1) res = res * a % MOD;
    return res;
}

ll sum(ll p, ll c)
{
    if (!c) return 1;
    if (c & 1) return (qpow(p, (c + 1) / 2) + 1) * sum(p, c / 2) % MOD;
    return ((qpow(p, c / 2) + 1) * sum(p, c / 2 - 1) + qpow(p, c)) % MOD;
}

void divide(ll n)
{
    for (ll i = 2; i * i <= n; ++i) {
        if (n % i == 0) {
            ll num = 0;
            while (n % i == 0) {
                num++;
                n /= i;
            }
            fac[++cnt] = make_pair(i, num);
        }
    }
    if (n > 1) fac[++cnt] = make_pair(n, 1);
}

int main()
{
    ll a, b;
    scanf("%lld%lld", &a, &b);
    divide(a);
    ll ans = 1;
    for (ll i = 1; i <= cnt; ++i) {
        ll p = fac[i].first, c = fac[i].second;
        ans = ans * sum(p, b * c) % MOD;
    }
    printf("%lld\n", a == 0 ? 0 : ans);
    return 0;
}