约数之和

假设现在有两个自然数 $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代码如下：

 1 #include <bits/stdc++.h>
 2 using namespace std;
 3 
 4 const int mod = 9901;
 5 
 6 int qmi(int a, int k) {
 7     int ret = 1;
 8     while (k) {
 9         if (k & 1) ret = ret * a % mod;
10         a = a * a % mod;
11         k >>= 1;
12     }
13     return ret;
14 }
15 
16 int main() {
17     int a, b;
18     scanf("%d %d", &a, &b);
19     if (a == 0) {
20         printf("0");
21     }
22     else if (b == 0) {
23         printf("1");
24     }
25     else {
26         unordered_map<int, int> mp;
27         for (int i = 2; i <= a / i; i++) {
28             while (a % i == 0) {
29                 mp[i]++;
30                 a /= i;
31             }
32         }
33         if (a > 1) mp[a]++;
34         int ret = 1;
35         for (auto &it : mp) {
36             if ((it.first - 1) % mod == 0) ret = ret * (it.second * b % mod + 1) % mod;
37             else ret = ret * (qmi(it.first % mod, it.second * b + 1) - 1) % mod * qmi((it.first - 1) % mod, mod - 2) % mod;
38         }
39         printf("%d", (ret + mod) % mod);
40     }
41     
42     return 0;
43 }

　　给出一种分治求$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代码如下：

 1 #include <bits/stdc++.h>
 2 using namespace std;
 3 
 4 const int mod = 9901;
 5 
 6 int qmi(int a, int k) {
 7     int ret = 1;
 8     while (k) {
 9         if (k & 1) ret = ret * a % mod;
10         a = a * a % mod;
11         k >>= 1;
12     }
13     return ret;
14 }
15 
16 int dfs(int p, int k) {
17     if (k == 0) return 1;
18     if (k + 1 & 1) return (p * dfs(p, k - 1) + 1) % mod;
19     return (1 + qmi(p, k / 2 + 1)) * dfs(p, k / 2) % mod;
20 }
21 
22 int main() {
23     int a, b;
24     scanf("%d %d", &a, &b);
25     if (a == 0) {
26         printf("0");
27     }
28     else if (b == 0) {
29         printf("1");
30     }
31     else {
32         unordered_map<int, int> mp;
33         for (int i = 2; i <= a / i; i++) {
34             while (a % i == 0) {
35                 mp[i]++;
36                 a /= i;
37             }
38         }
39         int ret = 1;
40         if (a > 1) mp[a]++;
41         for (auto &it : mp) {
42             ret = ret * dfs(it.first % mod, it.second * b) % mod;
43         }
44         printf("%d", ret);
45     }
46     
47     return 0;
48 }

参考资料

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