约数之和

约数之和

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

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

输入格式

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

输出格式

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

数据范围

$0 \leq A,B \leq 5 \times {10}^7$

输入样例：

2 3

输出样例：

15 

注意: $A$ 和 $B$ 不会同时为 $0$。

 

解题思路

　　把$A$分解成质因数乘积得到$A = P_1^{\alpha_1} \cdot P_2^{\alpha_2} \cdots P_n^{\alpha_n}$，因此有$A^B = P_1^{\alpha_1B} \cdot P_2^{\alpha_2B} \cdots P_n^{\alpha_nB}$，$A^B$的约数之和就是 $$\prod\limits_{i=1}^{n}{\sum\limits_{j=0}^{\alpha_iB}{P_i^j}} = \left( {P_1^0 + P_1^1 + \cdots  + P_1^{\alpha_1B}} \right) \cdot \left( {P_2^0 + P_2^1 + \cdots  + P_2^{\alpha_2B}} \right) \cdots \left( {P_n^0 + P_n^1 + \cdots  + P_n^{\alpha_nB}} \right)$$

　　根据等比级数求和公式，有$P_i^0 + P_i^1 + \cdots  + P_i^{\alpha_iB} = \dfrac{P_i^{\alpha_iB + 1} - 1}{P_i - 1}$。

　　其中分子直接通过快速幂求得。由于涉及到取模，因此需要通过费马小定理求分母$P_i - 1$的逆元（$9901$是质数）。问题是如果$P_i - 1$与$9901$不互质，即$P_i - 1$是$9901$的倍数，那么$P_i - 1$是不存在乘法逆元的。此时有$P_i - 1 \bmod 9901 = 0$，即有$P_i \equiv 1 \pmod{9901}$，因此

$$\begin{align*} & \ \ \ \ \ \sum\limits_{j=0}^{\alpha_iB}{P_i^j} \bmod 9901 \\ &= \sum\limits_{j = 0}^{\alpha_iB}{1} \bmod 9901 \\ &= (\alpha_iB + 1) \bmod 9901 \end{align*}$$

　　而如果$P_i - 1$与$9901$互质，那么直接求逆元然后相乘就可以了。

　　AC代码如下：

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

const int mod = 9901;

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

int main() {
    int a, b;
    scanf("%d %d", &a, &b);
    if (a == 0) {
        printf("0");
    }
    else if (b == 0) {
        printf("1");
    }
    else {
        unordered_map<int, int> mp;
        for (int i = 2; i <= a / i; i++) {
            while (a % i == 0) {
                mp[i]++;
                a /= i;
            }
        }
        if (a > 1) mp[a]++;
        int ret = 1;
        for (auto &it : mp) {
            if ((it.first - 1) % mod == 0) ret = ret * (it.second * b % mod + 1) % mod;
            else ret = ret * (qmi(it.first % mod, it.second * b + 1) - 1) % mod * qmi((it.first - 1) % mod, mod - 2) % mod;
        }
        printf("%d", (ret + mod) % mod);
    }
    
    return 0;
}

　　给出一种分治求$P_i^0 + P_i^1 + \cdots  + P_i^k$的做法。

　　如果$k + 1$是偶数（有偶数项相加），那么把$P_i^0 + P_i^1 + \cdots  + P_i^k$看作两部分相加，即

$$\begin{align*} & \ \ \ \ \ \ P_i^0 + P_i^1 + \cdots + P_i^k \\\\ &= \left( P_i^0 + P_i^1 + \cdots + P_i^{\left\lfloor \frac{k}{2} \right\rfloor} \right) + \left( {P_i^{\left\lfloor \frac{k}{2} \right\rfloor + 1} + P_i^{\left\lfloor \frac{k}{2} \right\rfloor + 2} + \cdots + P_i^k} \right) \\\\ &= \left( P_i^0 + P_i^1 + \cdots + P_i^{\left\lfloor \frac{k}{2} \right\rfloor} \right) + P_i^{\left\lfloor \frac{k}{2} \right\rfloor + 1} \cdot \left( P_i^0 + P_i^1 + \cdots + P_i^{\left\lfloor \frac{k}{2} \right\rfloor} \right) \\\\ &= \left( {1 + P_i^{\left\lfloor \frac{k}{2} \right\rfloor + 1}} \right) \cdot \left( P_i^0 + P_i^1 + \cdots + P_i^{\left\lfloor \frac{k}{2} \right\rfloor} \right) \end{align*}$$

　　其中记$\text{dfs}(p, k) = p^0 + p^1 + \cdots + p^k$，因此等式变成$\left( {1 + P_i^{\left\lfloor \frac{k}{2} \right\rfloor + 1}} \right) \cdot \text{dfs} \left(P_i, \left\lfloor \frac{k}{2} \right\rfloor \right)$，问题规模就缩小了一半。

　　如果$k + 1$是奇数，那么我们可以提取出一个$P_i$，把项式变成偶数，即$P_i^0 + P_i^1 + \cdots + P_i^k = 1 + P_i\times \text{dfs} \left( {P_i, k-1} \right)$。

　　AC代码如下：

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

const int mod = 9901;

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

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

int main() {
    int a, b;
    scanf("%d %d", &a, &b);
    if (a == 0) {
        printf("0");
    }
    else if (b == 0) {
        printf("1");
    }
    else {
        unordered_map<int, int> mp;
        for (int i = 2; i <= a / i; i++) {
            while (a % i == 0) {
                mp[i]++;
                a /= i;
            }
        }
        int ret = 1;
        if (a > 1) mp[a]++;
        for (auto &it : mp) {
            ret = ret * dfs(it.first % mod, it.second * b) % mod;
        }
        printf("%d", ret);
    }
    
    return 0;
}

　　再补充个矩阵乘法加快速幂的做法。

　　定义$F_i = \begin{bmatrix} S_i & p^i \end{bmatrix}$，其中$S_i = \sum\limits_{j=0}^{i}{p^j}$。为了得到$F_{i+1} = \begin{bmatrix} S_{i+1} & p^{i+1} \end{bmatrix}$，构造矩阵$A = \begin{bmatrix} 1 & 0 \\ p & p \end{bmatrix}$，那么就会有 $$\begin{bmatrix} S_i & p^i \end{bmatrix} \times \begin{bmatrix} 1 & 0 \\ p & p \end{bmatrix} = \begin{bmatrix} S_{i+1} & p^{i+1} \end{bmatrix}$$

　　其中初始状态$F_0 = \begin{bmatrix} 1 & 1 \end{bmatrix}$，那么$F_n = F_0 \times A^n$。

　　为了方便这里把向量扩展成$2 \times 2$的矩阵即$F_i = \begin{bmatrix} S_i & p^i \\ 0 & 0 \end{bmatrix}$。

　　AC代码如下：

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

const int mod = 9901;

void mul(int c[][2], int a[][2], int b[][2]) {
    int tmp[2][2] = {0};
    for (int i = 0; i < 2; i++) {
        for (int j = 0; j < 2; j++) {
            for (int k = 0; k < 2; k++) {
                tmp[i][j] = (tmp[i][j] + 1ll * a[i][k] * b[k][j]) % mod;
            }
        }
    }
    memcpy(c, tmp, sizeof(tmp));
}

int main() {
    int a, b;
    scanf("%d %d", &a, &b);
    if (a == 0) {
        printf("0");
    }
    else if (b == 0) {
        printf("1");
    }
    else {
        unordered_map<int, int> mp;
        for (int i = 2; i <= a / i; i++) {
            while (a % i == 0) {
                mp[i]++;
                a /= i;
            }
        }
        if (a > 1) mp[a]++;
        int ret = 1;
        for (auto it : mp) {
            int p = it.first, k = it.second * b;
            int f[2][2] = {{1, 1}, {0, 0}}, a[2][2] = {{1, 0}, {p, p}};
            while (k) {
                if (k & 1) mul(f, f, a);
                mul(a, a, a);
                k >>= 1;
            }
            ret = ret * f[0][0] % mod;
        }
        printf("%d", ret);
    }
    
    return 0;
}

 

参考资料

　　AcWing 97. 约数之和：https://www.acwing.com/video/116/

　　AcWing 97. 约数之和（两种解法）：https://www.acwing.com/solution/content/170582/