解法1：0x03 使用递归和分治算法解决取模乘法下等比数列的求和问题

#include <bits/stdc++.h>

using namespace std;

typedef long long ll;

ll a, b;
ll p[105], cnt[105];
ll mod = 9901;
int num = 0;

ll ksm(ll a, ll b){
    ll ans = 1;
    while(b){
        if(b & 1)
            ans = ans * a % mod;
        a = a * a % mod;
        b >>= 1;
    }
    return ans % mod;
}

ll sum(ll x, ll y){
    if(y == 1) return 1 + x;

    if(y % 2 == 1) return (1 + ksm(x, (y + 1)/2)) % mod * sum(x, y/2) % mod;
    else return ksm(x, y) + sum(x, y - 1) % mod;
}

int main(){
    cin >> a >> b;
    if(a == 0){
        cout << 0 << endl;
        return 0;
    }

    if(b == 0){
        cout << 1 << endl;
        return 0;
    }

    ll aa = a;
    for(int i = 2; i <= sqrt(a); ++ i){
        if(aa % i == 0){
            p[num] = i;
            int c = 0;
            while(aa % i == 0){
                c ++;
                aa /= i;
            }
            cnt[num ++] = c;
        }
    }
    if(aa != 1) p[num] = aa, cnt[num ++] = 1;

//    for(int i = 0; i < num; ++ i)
//        cout << p[i] << " " << cnt[i] << endl;

    ll res = 1;
    for(int i = 0; i < num; ++ i){
        res = res * sum(p[i], cnt[i] * b) % mod;
    }

    cout << res << endl;
    return 0;
}

解法2：0x33 使用递归和分治算法解决取模乘法下等比数列的求和问题

#include <bits/stdc++.h>

using namespace std;

int a, b, p[20], c[20], m, ans = 1, mod = 9901;

void divide(int n) {
    m = 0;
    for(int i = 2; i * i <= n; ++ i) {
        if(n % i == 0) {
            p[++ m] = i, c[m] = 0;
            while(n % i == 0) {
                c[m] ++;
                n /= i;
            }
        }
    }
    if(n > 1) p[++ m] = n, c[m] = 1;
}

int ksm(int a, int b) {
    int ret = 1;
    while(b) {
        if(b & 1) ret = 1ll * ret * a % mod;
        a = 1ll * a * a % mod;
        b >>= 1;
    }
    return ret;
}

int main(){
    cin >> a >> b;
    if(a == 0){
        cout << 0 << endl;
        return 0;
    }
    divide(a);
    for(int i = 1; i <= m; ++ i) {
        if((p[i] - 1) % mod == 0) {
            ans = (1ll * b * c[i] + 1) % mod * ans % mod;
            continue;
        }
        int x = ksm(p[i], b * c[i] + 1);
        x = (x - 1 + mod) % mod;
        int y = ksm(p[i] - 1, mod - 2);
        ans = 1ll * ans * x % mod * y % mod;
    }
    cout << ans << endl;
}