Project Euler 642

发表于 分类于 阅读次数： 本文字数： 3.9k 阅读时长 ≈ 4 分钟

Project Euler 642

题目

Sum of largest prime factors

Let \(f(n)\) be the largest prime factor of \(n\) and \(\displaystyle F(n) = \sum_{i=2}^{n}f(i)\).

For example \(F(10)=32\), \(F(100)=1915\) and \(F(10000)=10118280\).

Find \(F(201820182018)\). Give your answer modulus \(10^9\).

解决方案

注意，\(P^+(n)\)表示\(n\)的最大质因子，它和题目中的\(f(n)\)是等价的。

对任意 \(n\ge 2\)，令 \(p=f(n)\)。那么必存在整数 \(m\ge 1\) 使得\(n=p\cdot m.\)

并且有一个等价刻画：\(f(n)=p\Longleftrightarrow P^+(m)\le p.\)因为 \(p\) 是 \(n\) 的最大质因子，除了这一个因子 \(p\) 以外，剩下部分 \(m\) 的所有质因子都不能超过 \(p\)。

于是可以把总和写成

\[ F(N)=\sum_{\substack{p\le N\\p\text{ prime}}}p\cdot\#\left\{m\le \left\lfloor\frac{N}{p}\right\rfloor:P^+(m)\le p\right\}. \tag{1} \]

我们把上式里的先枚举 \(p\)改为先枚举 \(m\)。对于固定 \(m\ge 1\)，所有可行的 \(p\) 需要同时满足：\(p\) 是质数；\(p\le \lfloor N/m\rfloor\)（保证 \(pm\le N\)）；\(p\ge P^+(m)\)（保证 \(p\) 确实是最大质因子）。因此固定 \(m\) 的贡献为：\(\displaystyle{\sum_{\substack{p\text{ prime}\\P^+(m)\le p\le N/m}} p.}\)

而这就是一个质数区间和，我们令\(\displaystyle{\pi_1(n)=\sum_{\substack{p\le n\\p\text{ prime}}}p}\)，也就是质数和函数，那么上面\(m\)的贡献可用 \(\pi_1\) 写成闭式。

如果 \(A\) 是质数，则\(\displaystyle{\sum_{\substack{p\text{ prime}\\A\le p\le B}} p=\pi_1(B)-\pi_1(A)+A,}\)原因是 \(\pi_1(A)\) 包含了 \(A\) 本身，需要补回来才能变成闭区间 \([A,B]\)。

对 \(m>1\)，令 \(A=P^+(m)\ge 2\) 且为质数，于是\(m\)的贡献变为\(\pi_1\left(\left\lfloor\dfrac{N}{m}\right\rfloor\right)-\pi_1\left(P^+(m)\right)+P^+(m).\)

对 \(m=1\)，有 \(P^+(1)=1\)，贡献就是所有 \(\le N\) 的质数和：\(\displaystyle{\sum_{p\le N}p=\pi_1(N).}\)

将 \(m=1\) 单独提出，得到：

\[ F(N)=\pi_1(N)+\sum_{m\ge 2}\left(\pi_1\left(\left\lfloor\frac{N}{m}\right\rfloor\right)-\pi_1\left(P^+(m)\right)+P^+(m)\right), \tag{2} \]

我们接下来证明，只需要枚举到\(P^+(m)\le \sqrt N\)。要让固定 \(m\) 的区间 \([P^+(m),N/m]\) 非空，必须满足\(\left\lfloor\dfrac{N}{m}\right\rfloor \ge P^+(m)\Rightarrow\dfrac{N}{m}\ge P^+(m)\Rightarrow N\ge m\cdot P^+(m).\)

又因为 \(m\ge P^+(m)\)（最大质因子不可能大于整数本身），所以\(N\ge m\cdot P^+(m)\ge P^+(m)^2\Rightarrow P^+(m)\le \sqrt N.\)

也就是说，只有那些最大质因子不超过 \(\sqrt N\) 的 \(m\) 才会产生贡献。

因此式 \((2)\) 的求和只需要遍历满足\(m\ge 2, P^+(m)\le \sqrt N, m\cdot P^+(m)\le N\)的那些 \(m\)。

接下来我们尝试通过DFS枚举\(m\)。我们把一个 \(m\ge 2\) 写为质因数分解：\(\displaystyle{m=\prod_{i=1}^{t}p_i^{e_i},p_1>p_2>\cdots>p_t.}\)那么显然\(P^+(m)=p_1.\)

因此我们可以用先确定最大质因子 \(p_1\)，再往下乘更小的质数来无重无漏地枚举所有 \(m\)。

因为我们强制质数严格递减 \(p_1>p_2>\cdots\)，每个 \(m\) 的分解形式是唯一的，所以 DFS 路径也是唯一的：第一步选 \(p_1\)；之后只能选更小的质数；每个质数允许取任意正幂（对应循环乘 \(p,p^2,p^3,\dots\)）。因此不会出现重复计数。

当 DFS 构造出了某个 \(m\)，我们立刻累加贡献：\(\displaystyle{\pi_1\left(\left\lfloor\frac{N}{m}\right\rfloor\right)-\pi_1(P^+(m))+P^+(m).}\)

要保证贡献区间非空，必须满足：\(m\cdot P^+(m)\le N.\)令\(P=P^+(m)\)：当前 \(m\) 的最大质因子（当 \(m=1\) 时约定 \(P=1\)）。

当 DFS 正准备把一个新质数 \(p\) 乘进去时（并且我们始终让 \(p\le P\)，从而 \(P\) 始终表示最大质因子）：

  • 若当前 \(m=1\)（尚未选任何质数）：乘入 \(p\) 后得到新的 \(m=p^e\)，其最大质因子变为 \(P=p\)。为了保证以后至少能取到 \(p\) 作为区间端点（即区间不空），必须满足\((m\cdot p)\le N.\)
  • 若当前 \(m\ge2\)，其最大质因子为 \(P\)： 乘入更小的质数 \(p<P\) 后，新数变为 \(m\cdot p\)，但最大质因子仍为 \(P\)。为了保证新区间仍非空，需要\((m\cdot p)\cdot P \le N.\)

因此，DFS 在尝试扩展时，只要发现第二条分支不成立，就可以立刻停止这条分支，因为继续乘更小的质数只会让左边更大，不可能恢复可行。

接下来我们讲解高效快速计算\(\pi_1(x)\)。

令\(r=\left\lfloor\sqrt N\right\rfloor.\)在 DFS 中出现的自变量只会落在两类典型形式上：

  1. 小值区间：\(x\le r\)；
  2. 大值区间：\(x=\left\lfloor\dfrac{N}{i}\right\rfloor\)，其中 \(1\le i\le r\)。

因此我们只需要维护集合\(V=\left\{\left\lfloor\dfrac{N}{1}\right\rfloor,\left\lfloor\dfrac{N}{2}\right\rfloor,\dots,\left\lfloor\dfrac{N}{r}\right\rfloor\right\}\cup\{1,2,\dots,\left\lfloor N/r\right\rfloor\}.\)

该集合大小约为 \(2r\)，远小于 \(N\)。对每个 \(v\in V\)，我们维护一个值 \(S(v)\)，目标是在预处理结束后满足\(S(v)=\pi_1(v),(v\in V).\)

我们对任意 \(x\ge 1\)，先初始化为\(\displaystyle{S_0(x)=\sum_{k=2}^{x}k=\frac{x(x+1)}{2}-1.}\)随后通过筛法递推，把 \(S_0(x)\) 逐步修正成真正的质数和 \(\pi_1(x)\)。

当我们处理质数 \(p\)时，核心更新公式为：对所有 \(v\in V\) 且 \(v\ge p^2\)，\(S(v)\leftarrow S(v)-p\cdot\left(S(\lfloor v/p\rfloor)-\pi_1(p-1)\right).\)

这个式子的含义可以这样理解：

  • \(S(\lfloor v/p\rfloor)\) 表示在阈值尚未完全筛净时，\(\lfloor v/p\rfloor\) 范围内剩余对象的“总和估计”；
  • 其中小于 \(p\) 的质数部分 \(\pi_1(p-1)\) 不应被当作“以 \(p\) 为最小质因子产生的合数”重复扣除，因此需要减掉；
  • 整体乘上 \(p\)，对应把 \(\lfloor v/p\rfloor\) 中的对象乘上 \(p\) 映射回 \([1,v]\) 的那些数。

同时注意到：当 \(v<p^2\) 时，区间 \([1,v]\) 内已经不可能出现形如 \(p\cdot t\) 且 \(t\ge p\) 的合数，因此无需更新，递推自然截断。

当所有质数 \(p\le r\) 都被处理完之后，对于每个 \(v\in V\) 都有\(S(v)=\pi_1(v).\)

最终计算得到

\[ F(N)=\pi_1(N)+\sum_{\substack{m\ge 2,\\P^+(m)\le \sqrt N}}\left(\pi_1\left(\left\lfloor\frac{N}{m}\right\rfloor\right)-\pi_1(P^+(m))+P^+(m)\right) \]

代码

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
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
#include <bits/stdc++.h>
#include "tools.h"
using namespace std;

typedef long long ll;
const ll MOD = 1000000000LL;
const ll N = 201820182018LL;

// sum_{k=2..n} k = n(n+1)/2 - 1  (mod MOD)
ll sum_2_to_n_mod(ll n)
{
    if (n <= 1)
        return 0;
    ll a, b;
    if (n & 1)
    {
        a = (n + 1) / 2;
        b = n;
    }
    else
    {
        a = n / 2;
        b = n + 1;
    }
    ll res = (a % MOD) * (b % MOD) % MOD;
    res = (res - 1 + MOD) % MOD;
    return res;
}

// -------------------- Min_25 primesum precompute --------------------
ll g_n, g_r;
vector<ll> g_V;   // values: N//i and small numbers
vector<ll> g_S;   // S(x) = sum of primes <= x (mod MOD)
vector<int> g_pr; // primes up to sqrt(N)

// map x -> index in g_V
inline int idx_val(ll x)
{
    int m = (int)g_V.size();
    if (x <= g_r)
        return m - (int)x;
    return (int)(g_n / x - 1);
}

inline ll primesum(ll x)
{
    if (x < 2)
        return 0;
    return g_S[idx_val(x)];
}

void build_primesum_table(ll n)
{
    g_n = n;

    ll r = (ll)sqrt((long double)n);
    while ((r + 1) * (r + 1) <= n)
        ++r;
    while (r * r > n)
        --r;
    g_r = r;

    g_V.clear();
    for (ll i = 1; i <= r; i++)
        g_V.push_back(n / i);
    ll last = g_V.back();
    for (ll x = last - 1; x >= 1; x--)
        g_V.push_back(x);

    int m = (int)g_V.size();
    g_S.assign(m, 0);

    // init S(x) = sum_{k=2..x} k
    for (int i = 0; i < m; i++)
        g_S[i] = sum_2_to_n_mod(g_V[i]);

    g_pr = Factorizer(r).primes;

    // sieve:
    // S[v] -= p * ( S[v/p] - sum_{q<p} q )
    ll sum_prev = 0; // sum primes < p
    for (int p : g_pr)
    {
        ll p2 = 1LL * p * p;
        if (p2 > n)
            break;

        ll old_sum = sum_prev;
        sum_prev = (sum_prev + p) % MOD;

        for (int i = 0; i < m; i++)
        {
            ll v = g_V[i];
            if (v < p2)
                break;
            ll t = v / p;
            ll diff = (g_S[idx_val(t)] - old_sum + MOD) % MOD;
            g_S[i] = (g_S[i] - diff * p % MOD + MOD) % MOD;
        }
    }
}

// -------------------- Euler 642 DFS part --------------------
int g_lim;

ll dfs(int limit_idx, ll value, ll mp)
{
    ll res = 0;

    if (mp != 1)
    {
        ll a = primesum(g_n / value);
        ll b = primesum(mp);
        res = (a - b + mp + MOD) % MOD;
    }

    for (int i = 0; i < limit_idx; i++)
    {
        int p = g_pr[i];
        ll tmp = (value == 1 ? (ll)p : mp);

        // require value * p * tmp <= N
        if (value > g_n / ((ll)p * tmp))
            break;

        ll now = value * p;
        while (true)
        {
            res += dfs(i, now, tmp);
            if (res >= MOD)
                res %= MOD;

            if (now > g_n / ((ll)p * tmp))
                break;
            now *= p;
        }
    }
    return res % MOD;
}

int main()
{
    build_primesum_table(N);
    g_lim = int_square(N);
    ll ans = primesum(N);
    ans = (ans + dfs(g_lim, 1, 1)) % MOD;
    cout << ans << endl;
    return 0;
}