算法疑难杂症

【分治】A的B次方的所有约数之和

字数统计: 1.8k字   |   阅读时长: 8分
, ,

摘要: 综合题：分治+分解质因数+模下快速幂

【对算法，数学，计算机感兴趣的同学，欢迎关注我哈，阅读更多原创文章】
我的网站：潮汐朝夕的生活实验室
我的公众号：算法题刷刷
我的知乎：潮汐朝夕
我的github：FennelDumplings
我的leetcode：FennelDumplings

在文章 分治 中，我们系统学习了分治算法，了解了分治的思想和流程，以及研究分治算法时间复杂度的主定理，然后用分治算法解决了几道 leetcode 的题目，最后提到了分治与减治和搜索的区别和联系。在文章 减治法 中，我们讨论了减治法相关的话题。

本文我们来看一个跟数论有关的问题，乍一看跟减治没关系，但是经过分析和拆解后，就容易想到减治了。

问题

求 $A^{B}$ 的所有约数之和模 9901 的值。

1
1 <= A, B <= 5e7

算法: 分治

(1) n 的所有约数之和

首先我们要解决给定整数 N，它的所有约数之和怎么求。这个问题我们以前解决过。

在文章 约数 中，我们提到了算术基本定理及其推论，其中一条推论就是关于 N 的约数之和的。

这里我们把算术基本定理及其推论再罗列一下。

算术基本定理说的是对于任何一个大于 1 的自然数 N，如果 N 不为质数，那么 N 可以唯一分解成有限个质数的乘积 $N = P_{1}^{c_{1}}P_{2}^{c_{2}}…P_{k}^{c_{k}}$，称该式为 N 的标准分解式。

其中 $c_{i} \in N^{*}$, $p_{i}$ 为质数且 $p_{1} < p_{2} < … < p_{k}$。

它有三条推论

  • N 的正约数集合：$\{p_{1}^{b_{1}}p_{2}^{b_{2}}…p_{k}^{b_{k}}\}$，其中 $0 \leq b_{i} \leq c_{i}$
  • N 的正约数个数：$\prod\limits_{i=1}\limits^{k}(c_{i} + 1)$
  • N 的正约数和：$\prod\limits_{i=1}\limits^{k}(\sum\limits_{j=0}\limits^{c_{i}}(p_{i})^{j})$

因此假设 $A = P_{1}^{c_{1}}P_{2}^{c_{2}}…P_{k}^{c_{k}}$，$A^{B}$ 可以写为

它的约数集合如下

其中 $0 \leq b_{i} \leq Bc_{i}$。

由乘法分配律，$A^{B}$ 的约数之和就是以下形式，这也是算术基本定理的推论

(2) 分解质因数 = 素数筛 + 试除法

有了以上基于算术基本定理推论的 $A^{B}$ 的约数之和的表达式，我们就可以计算了。首先需要求出 A 的所有质因数 $p_{i}$。

在文章 素数筛 中，我们提到用结合试除法素性测试+筛法筛素数的方法求 N 的质因数的方法。

结合试除法素性测试+筛法筛素数，可以扫描 $2\sim\lfloor\sqrt{N}\rfloor$ 的每个 d，若 d 能整除 N，则从 N 中除掉所有的因子 d 并累计除去的 d 的个数。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
void divide(int n, vector<int>& p, vector<int>& c)
{
    p.clear();
    c.clear();
    for(int i = 2; i <= sqrt(n); ++i)
    {
        if(n % i == 0)
        {
            // i 是质数
            int cnt = 0;
            while(n % i == 0)
            {
                n /= i;
                ++cnt;
            }
            p.push_back(i);
            c.push_back(cnt);
        }
    }
    if(n > 1)
    {
        // n 是质数
        p.push_back(n);
        c.push_back(1);
    }
}

一个合数，其因子一定在扫描到该合数之前就被筛掉了，即在试除的过程中，能整除 N 的一定是质数。时间复杂度 $O(\sqrt{N})$

(3) 等比数列求和

现在我们有了 $A$ 的标准分解式，也就是 p 数组和 c 数组。

然后回到约数和的公式：

可以发现，对于每一个 i，我们要求的是一个等比数列的和 sum(p[i], c[i])

1
sum(p[i], c[i]) = 1 + p[i] + p[i]^2 + ... + p[i]^(B * c[i])

其中首项为 1，公比为 p[i]，项数为 1 + Bc[i]。可以用等比数列求和公式。

下面的问题就是实现这个求和，同时需要将结果模 MOD。

1
2
3
4
int sum(int p, int n)
{
    // 返回首项为 1，项数为 n + 1，公比为 p 的等比数列的和，模 MOD
}

(4) 模下快速幂

在实现 sum(p, n)时，求和公式中有一步求 $p[i]^{(1 + Bc[i])}$，这需要模下快速幂，具体细节参考 快速幂 这篇文章。这里直接给出代码模板

1
2
3
4
5
6
7
8
9
10
11
12
using ll = long long;
int quickpower2(int a, int n, int MOD) // 处理不了n是负数
{
    int ans = 1;
    while(n)
    {
        if(n & 1) ans = ((ll)ans * a) % MOD;
        a = (ll)a * a % MOD;
        n >>= 1;
    }
    return ans % MOD;
}

(5) 减治

求和公式中的求幂解决了，但是还有除法的问题，有除法而答案需要对 9901 取模。MOD 只对加减乘有分配律，不能直接对分子分母取模后再做除法。

那就只能用一项一项加的方式实现 sum(p, n) 了。

通过观察发现，前半部分的和好像可以直接用到后半部分，于是就想到了用减治法进行等比数列求和。

如果 n 为奇数

如果 n 为偶数

1
2
3
4
5
6
7
8
9
int sum(int p, int n)
{
    // 首项为 1，项数为 n + 1，公比为 p 的等比数列的和，模 MOD
    if(n == 0)
        return 1;
    if((n & 1) == 0)
        return ((ll)quickpower2(p, n, MOD) + sum(p, n - 1)) % MOD;
    return ((1 + (ll)quickpower2(p, (n + 1) / 2, MOD)) * sum(p, (n - 1) / 2)) % MOD;
}

代码 (C++)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
#include <iostream>
#include <vector>
#include <cmath>

using namespace std;

void divide(int n, vector<int>& p, vector<int>& c)
{
    p.clear();
    c.clear();
    for(int i = 2; i <= sqrt(n); ++i)
    {
        if(n % i == 0)
        {
            // i 是质数
            int cnt = 0;
            while(n % i == 0)
            {
                n /= i;
                ++cnt;
            }
            p.push_back(i);
            c.push_back(cnt);
        }
    }
    if(n > 1)
    {
        // n 是质数
        p.push_back(n);
        c.push_back(1);
    }
}

const int MOD = 9901;
using ll = long long;

int quickpower2(int a, int n, int MOD) // 处理不了n是负数
{
    int ans = 1;
    while(n)
    {
        if(n & 1) ans = ((ll)ans * a) % MOD;
        a = (ll)a * a % MOD;
        n >>= 1;
    }
    return ans % MOD;
}

int sum(int p, int n)
{
    // 首项为 1，项数为 n + 1，公比为 p 的等比数列的和，模 MOD
    if(n == 0)
        return 1;
    if((n & 1) == 0)
        return ((ll)quickpower2(p, n, MOD) + sum(p, n - 1)) % MOD;
    return ((1 + (ll)quickpower2(p, (n + 1) / 2, MOD)) * sum(p, (n - 1) / 2)) % MOD;
}

int main()
{
    int A, B;
    cin >> A >> B;
    if(A == 0)
        cout << 0 << endl;
    else if(B == 0)
        cout << 1 << endl;
    else
    {
        vector<int> p, c;
        divide(A, p, c);
        int ans = 1;
        int k = p.size();
        for(int i = 0; i < k; ++i)
            ans = ((ll)ans * sum(p[i], (ll)B * c[i])) % MOD;
        cout << ans << endl;
    }
}
Share