题目来源：AcWing 97. 约数之和

题目描述

假设现在有两个自然数 $A$ 和 $B$ ， $S$ 是 $A^B$ 的所有约数之和。

请你求出 $S$ mod $9901$ 的值是多少。

输入格式

在一行中输入用空格隔开的两个整数 $A$ 和 $B$ 。

输出格式

输出一个整数，代表 $S$ mod $9901$ 的值。

数据范围

$0≤A,B≤5×10^7$
$A$ 和 $B$ 不会同时为0

输入样例

2 3

输出样例

解题思路

计算约数的和

为了简化问题，我们先来看如何对 $A$ 求约数之和。设存在质数 $p_i$ 为A的质因数，则对于 $A$ 有:

A = {p_1}^{k_1} + {p_2}^{k_2} + ... + {p_n}^{k_n}

以质因数 ${p_1}$ 为例， $A$ 拥有 ${k_1}$ 个 ${p_1}$ ，故对于质因数 ${p_1}$ 存在 ${k_1} + 1$ 种选法，组成不同的因数。因此， $A$ 的因数个数 $Count_{A}$ 的值为（约数个数定理）：

Count_{A} = ({k_1} + 1)({k_2} + 1)...({k_n} + 1)

由 $A$ 的因数个数 $Count_{A}$ 和乘法原理，则可以得到 $A$ 的约数之和 $Sum_{A}$ 的值（约数和定理）:

Sum_{A} = ({p_1}^0 + {p_1}^1 + ... + {p_{1}^{k_1}})({p_2}^0 + {p_2}^1 + ... + {p_{2}^{k_2}})...({p_n}^0 + {p_n}^1 + ... + {p_{n}^{k_n}})

其中，任意的 ${p_{i}^{k_j}}$ 均为 $A$ 的约数。将该式推广，则可以得到 $Sum_{A^B}$ 的值:

Sum_{A^B} = ({p_1}^0 + {p_1}^1 + ... + {p_{1}^{B*k_1}})({p_2}^0 + {p_2}^1 + ... + {p_{2}^{B*k_2}})...({p_n}^0 + {p_n}^1 + ... + {p_{n}^{B*k_n}})

其中，任意的 ${p_{i}^{B*k_j}}$ 均为 $A^B$ 的约数，每一项对应一个质因数 ${p_i}$ ，分别对每一项求和并累乘即可得到结果。根据上述思路，我们给出相应的代码。

int main() {
    int a, b;
    cin>>a>>b;
    int res = 1;
    for(int x = 2; x <= a; x++) { // 求质因数
        int cnt = 0;
        while(a % x == 0) {
            cnt += 1;
            a /= x;
        }
        res = res * sum(x, cnt * b) % mod; // 拥有b*cnt个质因数x
    }
    if(a == 0) res = 0;
    cout<<res<<endl;
    return 0;
}

分治求解每一项的和

对于质因数 $p_i$ 的项 ${p_i}^0 + {p_i}^1 + ... + {p_{i}^N}$ ，使用分治法求和能够将时间复杂度压缩到 $O(logN)$ 。

每一项从 $N/2$ 处分开，则原式被分为两部分，分别为前式 ${p_i}^0 + {p_i}^1 + ... + {p_{i}^{N/2}}$ 与后式 ${p_i}^{N/2 + 1} + {p_i}^{N/2 + 2} + ... + {p_{i}^{N}}$ 。

当 $N$ 为偶数时，从后式中提出公因子 ${p_i}^{N/2 + 1}$ ，则后式变为 ${p_i}^{N/2 + 1} * ({p_i}^0 + {p_i}^1 + ... + {p_{i}^{N/2}})$ ，整理上述式子。

{p_i}^0 + {p_i}^1 + ... + {p_{i}^N} = (1 + {p_i}^{N/2 + 1}) * ({p_i}^0 + {p_i}^1 + ... + {p_{i}^{N/2}})

设 $sum_{p_i}(N)={p_i}^0 + {p_i}^1 + ... + {p_{i}^N}$ ，则当 $N$ 为偶数时，有：

sum_{p_i}(N) = (1 + {p_i}^{N/2 + 1}) * sum_{p_i}(N / 2)

当 $N$ 为奇数时，则将 $N$ 减去1变为偶数，有：

sum_{p_i}(N) = {p_i} * sum_{p_i}(N-1) + {p_i}^0 = {p_i} * sum_{p_i}(N-1) + 1

我们给出对应的代码。

int sum(int p, int k) {
    if(k == 0) return 1;
    if(k % 2 == 0) return (p % mod * sum(p, k - 1) + 1) % mod;
    return (1 + qmi(p, k / 2 + 1)) * sum(p, k / 2) % mod;
}

快速幂求公因子 ${p_i}^{N/2 + 1}$

为了方便描述，我们将公因子简化为 $p^N$ 。对任意的整数 $N$ ，其一定可以转化为对应的二进制表示。设对应的二进制表示中的各个位都为1，则每一位二进制位都有值 $p^{2 * k}=p^k * p^k$ 。其中, $p^{2 * k}$ 为当前二进制位对应的值， $p^k$ 为前一位二进制位对应的值。

当二进制位为1时，累加当前值，当二进制位为0时，不加当前值。我们给出对应的代码。

int qmi(int p, int k) {
    p %= mod;
    int res = 1;
    while(k) {
        if(k & 1) res = res * p % mod;
        p = p * p % mod;
        k >>= 1;
    }
    return res;
}

参考代码

#include<iostream>
using namespace std;
const int mod = 9901;

int qmi(int p, int k) {
    p %= mod;
    int res = 1;
    while(k) {
        if(k & 1) res = res * p % mod;
        p = p * p % mod;
        k >>= 1;
    }
    return res;
}

int sum(int p, int k) {
    if(k == 0) return 1;
    if(k % 2 == 0) return (p % mod * sum(p, k - 1) + 1) % mod;
    return (1 + qmi(p, k / 2 + 1)) * sum(p, k / 2) % mod;
}

int main() {
    int a, b;
    cin>>a>>b;
    int res = 1;
    for(int x = 2; x <= a; x++) {
        int cnt = 0;
        while(a % x == 0) {
            cnt += 1;
            a /= x;
        }
        res = res * sum(x, cnt * b) % mod;
    }
    if(a == 0) res = 0;
    cout<<res<<endl;
    return 0;
}

AcWing 97. 约数之和【分治 数论】