POJ 1845 - Sumdiv（唯一分解+递推）

题目链接 https://vjudge.net/problem/POJ-1845

【题目描述】
Consider two natural numbers A and B. Let S be the sum of all natural divisors of A^B. Determine S modulo 9901 (the rest of the division of S by 9901).
Input
The only line contains the two natural numbers A and B, (0 <= A,B <= 50000000)separated by blanks.
Output
The only line of the output will contain S modulo 9901.
Sample Input
2 3
Sample Output
15
Hint
2^3 = 8.
The natural divisors of 8 are: 1,2,4,8. Their sum is 15.
15 modulo 9901 is 15 (that should be output).

【题意】
给定a,b，要求a^b的所有因子之和对9901取模后的结果（1 <= a, b <= 5e7）
bb
【思路】
先对a进行唯一分解，a=p1^a1×p2^a2×…×pn^an，那么a^b=p1^(a1×b)×p2^(a2×b)×…×pn^(an×b)，所以这个数的所有因子之和sum=(1+p1+p1^2+…+p1^(a1×b))×(1+p2+p2^2+…+p2^(a2×b))×(1+pn+pn^2+…+pn^(an×b))，每个括号里都是枚举选择每种素因子的个数，每个括号里的式子都是一个等比数列，可以递归地求解sum(x,n)=1+x+x^2+…+x^n，当n为奇数时y=(1+x+x^2+…+x^(n/2))+(x^(n/2+1)+x^(n/2+2)+…+x^n)，提公因式就能得到
sum(x,n)=(1+x+x^2+…+x^(n/2))×(1+x^(n/2+1))=sum(x,n/2)×(1+x^(n/2+1))，
同理当n为偶数时可以同同样的方法推出
sum(x,n)=sum(x,n/2-1)×(1+x^(n/2+1))+x^(n/2)，这里的除法’/’指的都是整除。

#include<cmath>
#include<cstdio>
#include<cstring>
#include<iostream>
#include<algorithm>
using namespace std;
typedef long long ll;
const int mod = 9901;

ll a, b;
int fac[10050], cnt[10050], num;

void init() {
    num = 0;
    memset(cnt, 0, sizeof(cnt));
    int m = (int)sqrt(a + 0.5);
    ll cpy = a;
    for (int i = 2; i <= m; ++i) {
        if (cpy % i == 0) {
            fac[num] = i;
            while (cpy%i == 0) {
                cpy /= i;
                ++cnt[num];
            }
            ++num;
        }
    }
    if (cpy > 1) {
        fac[num] = cpy;
        ++cnt[num];
        ++num;
    }
}

ll pow_mod(ll x, ll n) {
    ll ans = 1;
    while (n) {
        if (n & 1) ans = (ans * x) % mod;
        x = (x * x) % mod;
        n >>= 1;
    }
    return ans;
}

ll sum(ll x, ll n) {//计算1+x+x^2+...+x^n
    if (0 == n) return 1;
    if (n & 1) {
        return sum(x,n/2) * (1+pow_mod(x,n/2+1)) % mod;
    }
    else {
        return (sum(x, n/2-1) * (1+pow_mod(x,n/2+1)) + pow_mod(x, n / 2)) % mod;
    }
}

int main() {
    while (scanf("%lld%lld", &a, &b) == 2) {
        init();
        ll ans = 1;
        for (int i = 0; i < num; ++i) {
            ans = (ans * sum(fac[i], cnt[i] * b)) % mod;
        }
        printf("%lld\n", ans);
    }
    return 0;
}

这道题容易想到的还是用等比数列求和公式去做，1+x+x^2+…+x^n=(x^(n+1)-1)/(x-1)，但是这涉及到了除法，虽然9901是素数，正常情况下可以用逆元去做，但是如果我们要知道一个数a模9901的逆元，而a恰好是9901的整数倍，那么逆元就不存在了，所以不能直接用逆元。但我们在计算一个分式A/B(mod p)的时候，有A/B(mod p)=Amod(B*p)/B，可以通过这个式子来求解，注意快速幂要用大数乘法，这道题的数据会爆long long。

#include<cmath>
#include<cstdio>
#include<cstring>
#include<iostream>
#include<algorithm>
using namespace std;
typedef long long ll;
const int Mod = 9901;

ll fac[10050], cnt[10050], num;

void init(ll a) {
    num = 0;
    memset(cnt, 0, sizeof(cnt));
    int m = (int)sqrt(a + 0.5);
    ll cpy = a;
    for (int i = 2; i <= m; ++i) {
        if (cpy % i == 0) {
            fac[num] = i;
            while (cpy%i == 0) {
                cpy /= i;
                ++cnt[num];
            }
            ++num;
        }
    }
    if (cpy > 1) {
        fac[num] = cpy;
        ++cnt[num];
        ++num;
    }
}

ll mul_mod(ll a, ll b, ll mod) {//大数乘法取模
    ll ans = 0;
    while (b) {
        if (b & 1) ans = (ans + a) % mod;
        a = (a + a) % mod;
        b >>= 1;
    }
    return ans;
}

ll pow_mod(ll x, ll n, ll mod) {
    ll ans = 1;
    while (n) {
        if (n & 1) ans = mul_mod(ans, x, mod);
        x = mul_mod(x, x, mod);
        n >>= 1;
    }
    return ans;
}

int main() {
    ll a, b;
    while (scanf("%lld%lld", &a, &b) == 2) {
        init(a);
        ll ans = 1;
        for (int i = 0; i < num; ++i) {
            ll M = Mod * (fac[i] - 1);
            ans *= (pow_mod(fac[i], cnt[i] * b + 1, M) - 1 + M) / (fac[i] - 1);
            ans %= Mod;
        }
        printf("%lld\n", ans);
    }
    return 0;
}