题目链接

一个数 A 的约数之和可以由分解质因数求得，一个质因数 p 在 A 中出现 c 次，那么对于一个约数 X ，X 中 p 出现的次数一定为 0次 ~ c次，所以对于每个 p 可以构造一个等比数列 1 + p + p^2 + p^3 + ...... +p^c，最后将所有的等比数列之和乘起来即可。即 ans = ans * sum(p[i], c[p[i]]) %?mod; 其中 p[i] 为 A 的质因数， c[p[i]] 为 p[i] 在 A 中出现的次数， sum为首项为 1 ，尾项为 (p[i] ^ c[p[i]]) 的等比数列求和函数。

那么 A的B次方 求约数之和就是在原来每个质因数出现的次数乘上 B 即可。

接下来是三种作法：

1. 费马小定理。

#include <bits/stdc++.h>
using namespace std;

typedef long long ll;
typedef unsigned long long ull;
#define int ll
#define PI acos(-1.0)
#define INF 0x3f3f3f3f3f3f3f3f
#define P pair<int, int>
#define fir first
#define sec second
#define fastio ios::sync_with_stdio(false), cin.tie(0)
const int mod = 9901;
const int MAXM = 50000 + 10;
const int MAXN = 1000000 + 10;

int a, b;
int p[MAXN], c[MAXN], cnt;

int POW(int n, int m)
{
    int ans = 1;
    while(m) {
        if(m & 1) ans = ans * n % mod;
        n = n * n % mod;
        m >>= 1;
    }
    return ans % mod;
}

void deal()
{
    for(int i = 2; i <= a; i ++) {
        if(a % i == 0) {
            p[cnt++] = i;
            while(a % i == 0) {
                a /= i;
                c[i] ++;
            }
        }
    }
}

signed main()
{
    cin >> a >> b;
    if(a == 0) {
        cout << 0 << endl;
        return 0;
    }
    deal();
    int ans = 1;
    for(int i = 0; i < cnt; i ++) {
        if(p[i] - 1 % mod == 0) {
            ans = ans * (b * c[p[i]] + 1) % mod;
        }
        else {
            ans = ans * ((POW(p[i], b * c[p[i]] + 1) - 1 + mod) % mod) * (POW(p[i] - 1, mod - 2)) % mod;
        }
    }
    cout << ans << endl;
}

2. 矩阵快速幂

#include <bits/stdc++.h>
using namespace std;

typedef long long ll;
typedef unsigned long long ull;
#define int ll
#define PI acos(-1.0)
#define INF 0x3f3f3f3f3f3f3f3f
#define P pair<int, int>
#define fir first
#define sec second
#define fastio ios::sync_with_stdio(false), cin.tie(0)
const int mod = 9901;
const int M = 2 + 10;
const int N = 1000000 + 10;

int a, b;
int p[N], c[N], cnt;

void deal()
{
    for(int i = 2; i <= a; i ++) {
        if(a % i == 0) {
            p[cnt++] = i;
            while(a % i == 0) {
                a /= i;
                c[i] ++;
            }
        }
    }
}

struct Mat {
    int n;
    int m[M][M];
};
Mat operator *(Mat a,Mat b)
{
    int n = 2;
    Mat c;
    memset(c.m, 0, sizeof c.m);
    for(int k = 0; k < n; k ++){
        for(int i = 0; i < n; i ++){
            if(a.m[i][k])
                for(int j = 0; j < n; j ++)
                    if(b.m[k][j])
                        c.m[i][j] = (c.m[i][j] + a.m[i][k] * b.m[k][j]) % mod;
        }
    }
    return c;
}
Mat operator ^(Mat a, int k)
{
    int n = 2;
    Mat c;
    memset(c.m, 0, sizeof c.m);
    for(int i = 0; i < n; i ++){
        c.m[i][i] = 1;
    }
    for(; k; k >>= 1)
    {
        if(k&1) c = c * a;
        a = a * a;
    }
    return c;
}

signed main()
{
    cin >> a >> b;
    if(a == 0) {
        cout << 0 << endl;
        return 0;
    }
    deal();
    int ans = 1;
    for(int i = 0; i < cnt; i ++) {
        Mat start, ed, mid;
        
        mid.m[0][0] = 1; mid.m[0][1] = 0;
        mid.m[1][0] = p[i]; mid.m[1][1] = p[i];
        
        start.m[0][0] = 1; start.m[0][1] = 1;
        
        ed = start * (mid ^ (c[p[i]] * b));
        ans = ans * ed.m[0][0] % mod;
    }
    cout << ans << endl;
}

?

3. 分治

#include <bits/stdc++.h>
using namespace std;

typedef long long ll;
typedef unsigned long long ull;
#define int ll
#define PI acos(-1.0)
#define INF 0x3f3f3f3f3f3f3f3f
#define P pair<int, int>
#define fir first
#define sec second
#define fastio ios::sync_with_stdio(false), cin.tie(0)
const int mod = 9901;
const int M = 2 + 10;
const int N = 1000000 + 10;

int a, b;
int p[N], c[N], cnt;

void deal()
{
    for(int i = 2; i <= a; i ++) {
        if(a % i == 0) {
            p[cnt++] = i;
            while(a % i == 0) {
                a /= i;
                c[i] ++;
            }
        }
    }
}

int POW(int n, int m)
{
    int ans = 1;
    while(m) {
        if(m & 1) ans = ans * n % mod;
        n = n * n % mod;
        m >>= 1;
    }
    return ans % mod;
}

int divide(int q, int n)
{
    if(n == 0) return 1;
    if(n & 1) {
        return (1 + POW(q, (n + 1) / 2)) * divide(q, (n - 1) / 2) % mod;
    }
    else {
        return (POW(q, n) + (1 + POW(q, n / 2)) * divide(q, (n - 2) / 2)) % mod;
    }
}

signed main()
{
    cin >> a >> b;
    if(a == 0) {
        cout << 0 << endl;
        return 0;
    }
    deal();
    int ans = 1;
    for(int i = 0; i < cnt; i ++) {
        ans = ans * divide(p[i], c[p[i]] * b) % mod;
    }
    cout << ans << endl;
}

?