约数之和

约数之和

假设现在有两个自然数 ABSAB 的所有约数之和。

请你求出 Smod9901 的值是多少。

输入格式

在一行中输入用空格隔开的两个整数 AB

输出格式

输出一个整数,代表 Smod9901 的值。

数据范围

0A,B5×107

输入样例:

2 3

输出样例:

15 

注意: AB 不会同时为 0

 

解题思路

  把A分解成质因数乘积得到A=P1α1P2α2Pnαn,因此有AB=P1α1BP2α2BPnαnBAB的约数之和就是i=1nj=0αiBPij=(P10+P11++P1α1B)(P20+P21++P2α2B)(Pn0+Pn1++PnαnB)

  根据等比级数求和公式,有Pi0+Pi1++PiαiB=PiαiB+11Pi1

  其中分子直接通过快速幂求得。由于涉及到取模,因此需要通过费马小定理求分母Pi1的逆元(9901是质数)。问题是如果Pi19901不互质,即Pi19901的倍数,那么Pi1是不存在乘法逆元的。此时有Pi1mod9901=0,即有Pi1(mod9901),因此

     j=0αiBPijmod9901=j=0αiB1mod9901=(αiB+1)mod9901

  而如果Pi19901互质,那么直接求逆元然后相乘就可以了。

  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 }
复制代码

  给出一种分治求Pi0+Pi1++Pik的做法。

  如果k+1是偶数(有偶数项相加),那么把Pi0+Pi1++Pik看作两部分相加,即

      Pi0+Pi1++Pik=(Pi0+Pi1++Pik2)+(Pik2+1+Pik2+2++Pik)=(Pi0+Pi1++Pik2)+Pik2+1(Pi0+Pi1++Pik2)=(1+Pik2+1)(Pi0+Pi1++Pik2)

  其中记dfs(p,k)=p0+p1++pk,因此等式变成(1+Pik2+1)dfs(Pi,k2),问题规模就缩小了一半。

  如果k+1是奇数,那么我们可以提取出一个Pi,把项式变成偶数,即Pi0+Pi1++Pik=1+Pi×dfs(Pi,k1)

  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 }
复制代码

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

  定义Fi=[Sipi],其中Si=j=0ipj。为了得到Fi+1=[Si+1pi+1],构造矩阵A=[10pp],那么就会有[Sipi]×[10pp]=[Si+1pi+1]

  其中初始状态F0=[11],那么Fn=F0×An

  为了方便这里把向量扩展成2×2的矩阵即Fi=[Sipi00]

  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/

posted @   onlyblues  阅读(33)  评论(0编辑  收藏  举报
相关博文:
阅读排行:
· 单线程的Redis速度为什么快?
· 展开说说关于C#中ORM框架的用法!
· Pantheons:用 TypeScript 打造主流大模型对话的一站式集成库
· SQL Server 2025 AI相关能力初探
· 为什么 退出登录 或 修改密码 无法使 token 失效
历史上的今天:
2022-02-18 地宫取宝
Web Analytics
点击右上角即可分享
微信分享提示