约数之和
假设现在有两个自然数 $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 }
再补充个矩阵乘法加快速幂的做法。
定义$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代码如下：
1 #include <bits/stdc++.h> 2 using namespace std; 3 4 const int mod = 9901; 5 6 void mul(int c[][2], int a[][2], int b[][2]) { 7 int tmp[2][2] = {0}; 8 for (int i = 0; i < 2; i++) { 9 for (int j = 0; j < 2; j++) { 10 for (int k = 0; k < 2; k++) { 11 tmp[i][j] = (tmp[i][j] + 1ll * a[i][k] * b[k][j]) % mod; 12 } 13 } 14 } 15 memcpy(c, tmp, sizeof(tmp)); 16 } 17 18 int main() { 19 int a, b; 20 scanf("%d %d", &a, &b); 21 if (a == 0) { 22 printf("0"); 23 } 24 else if (b == 0) { 25 printf("1"); 26 } 27 else { 28 unordered_map<int, int> mp; 29 for (int i = 2; i <= a / i; i++) { 30 while (a % i == 0) { 31 mp[i]++; 32 a /= i; 33 } 34 } 35 if (a > 1) mp[a]++; 36 int ret = 1; 37 for (auto it : mp) { 38 int p = it.first, k = it.second * b; 39 int f[2][2] = {{1, 1}, {0, 0}}, a[2][2] = {{1, 0}, {p, p}}; 40 while (k) { 41 if (k & 1) mul(f, f, a); 42 mul(a, a, a); 43 k >>= 1; 44 } 45 ret = ret * f[0][0] % mod; 46 } 47 printf("%d", ret); 48 } 49 50 return 0; 51 }
参考资料
AcWing 97. 约数之和：https://www.acwing.com/video/116/
AcWing 97. 约数之和（两种解法）：https://www.acwing.com/solution/content/170582/
原文链接: https://www.cnblogs.com/onlyblues/p/17132273.html
欢迎关注
微信关注下方公众号，第一时间获取干货硬货；公众号内回复【pdf】免费获取数百本计算机经典书籍
原创文章受到原创版权保护。转载请注明出处：https://www.ccppcoding.com/archives/318967
非原创文章文中已经注明原地址，如有侵权，联系删除
关注公众号【高性能架构探索】，第一时间获取最新文章
转载文章受原作者版权保护。转载请注明原作者出处！
