Project Euler 521

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

Project Euler 521

题目

Smallest prime factor

Let \(\text{smpf}(n)\) be the smallest prime factor of \(n\).

\(\text{smpf}(91)=7\) because \(91=7\times13\) and \(\text{smpf}(45)=3\) because \(45=3\times3\times5\).

Let \(S(n)\) be the sum of \(\text{smpf(i)}\) for \(2 \le i \le n\).

E.g. \(S(100)=1257\).

Find \(S(10^{12}) \bmod 10^9\).

解决方案

对“没有小于 \(p\) 的质因子”的数，定义两类函数（都只统计 \([2,x]\)，不含 \(1\)）：

\[ C_p(x)=\#\{m\in\mathbb Z:2\le m\le x,\operatorname{smpf}(m)\ge p\}, A_p(x)=\sum_{\substack{2\le m\le x\\\operatorname{smpf}(m)\ge p}} m. \]

并令\(\Omega\)是质数集合。

接下来先对合数部分进行计数，对固定质数 \(p\)，满足 \(\operatorname{smpf}(n)=p\) 且 \(n\) 为合数的数恰为 \(n=p\cdot m\)，其中

  • \(m\ge p\)（保证 \(n\ge p^2\) 合数），
  • \(\operatorname{smpf}(m)\ge p\)，
  • \(m\le \lfloor N/p\rfloor\)。

因此这类 \(m\) 的个数为\(C_p\left(\left\lfloor\dfrac{N}{p}\right\rfloor\right)-C_p(p-1).\) 它们对 \(S(N)\) 的贡献为 \[ p\cdot \left(C_p\left(\left\lfloor\frac{N}{p}\right\rfloor\right)-C_p(p-1)\right). \]

当我们把允许的最小质因子阈值”从 \(p\) 推进到下一个质数 \(p^+\)（也就是把质数 \(p\) 也排除掉），需要从集合里删掉所有“最小质因子为 \(p\) 的合数：这些数在 \([2,x]\) 中恰是 \(p\cdot t\)，其中 \(t\ge p\) 且 \(\operatorname{smpf}(t)\ge p\) 且 \(t\le \lfloor x/p\rfloor\)。因此：

\[ C_{p^+}(x)=C_p(x)-(C_p(\lfloor x/p\rfloor)-C_p(p-1)), A_{p^+}(x)=A_p(x)-p(A_p(\lfloor x/p\rfloor)-A_p(p-1)). \]

注意这里 刻意不删掉质数 \(p\) 本身：因为上式只删 \(p\cdot t\) 且 \(t\ge p\)，不会涉及 \(p\cdot 1\)。最终剩下来的恰好就是所有质数（它们的最小质因子就是它们自己），可以作为质数部分直接加回去。

设 \(v=\lfloor\sqrt N\rfloor\)。我们对上述\(p\)的处理，不会超过\(v\)，这是因为任意合数 \(n\le N\) 必有一个质因子 \(\le \sqrt n\le v\)，因此当我们把所有质数 \(p\le v\) 的“最小质因子为 \(p\) 的合数贡献”都加完后，剩下没覆盖的数只能是质数（以及 \(1\)，但我们从不统计 \(1\)）。

换言之：

  • 合数部分只需枚举 \(p\le v\)；
  • 最终剩下的 \(A_{v+1}(N)\) 就是 \(\displaystyle{\sum_{p\le N,p\in\Omega}p}\)（所有质数和）。

所以答案是 \[ S(N)= \left(\sum_{\substack{p\le v,p\in \Omega}}p\cdot \left(C_p(\lfloor N/p\rfloor)-C_p(p-1)\right)\right)+A_{v+1}(N). \]

在递推里会出现很多形如 \(\lfloor x/p\rfloor\) 的值。因此，对固定 \(N\)，所有可能出现的自变量只需要两类：

  • 小值：\(x\le v\)（直接数组下标存）。
  • 大值：\(x=\left\lfloor\dfrac{N}{i}\right\rfloor\)，其中 \(1\le i\le v\)（这类也只有 \(v\) 个）。

因此我们维护两套表：

  • s_cnt[x], s_sum[x] 表示当前阈值下的 \(C(x),A(x)\)，其中 \(x\in[0,v]\)；
  • l_cnt[i], l_sum[i] 表示当前阈值下的 \(C(\lfloor N/i\rfloor),A(\lfloor N/i\rfloor)\)，其中 \(i\in[1,v]\)。

初始化对应阈值 \(p=2\)（尚未筛掉任何质数）：\(\displaystyle{C_2(x)=x-1, A_2(x)=\sum_{k=2}^{x}k=\frac{x(x+1)}{2}-1.}\)

之后对每个质数 \(p\le v\) 用上一节的递推更新所有需要的 \(x\)。

代码

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
#include <bits/stdc++.h>
using namespace std;

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

// 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 - 1 + MOD) % MOD;
    return res;
}

ll S_min_prime_factor_mod(ll N)
{
    ll v = (ll)sqrt((long double)N);
    while ((v + 1) * (v + 1) <= N)
        ++v;
    while (v * v > N)
        --v;

    vector<ll> s_cnt(v + 1), s_sum(v + 1);
    for (ll i = 0; i <= v; i++)
    {
        s_cnt[i] = i - 1;             // exact
        s_sum[i] = sum_2_to_n_mod(i); // modded
    }

    vector<ll> l_cnt(v + 1), l_sum(v + 1);
    l_cnt[0] = 0;
    l_sum[0] = 0;
    for (ll i = 1; i <= v; i++)
    {
        ll x = N / i;
        l_cnt[i] = x - 1;             // exact
        l_sum[i] = sum_2_to_n_mod(x); // modded
    }

    ll ans = 0;

    for (ll p = 2; p <= v; p++)
    {
        // prime test: if s_cnt[p] != s_cnt[p-1], then p is prime
        if (s_cnt[p] == s_cnt[p - 1])
            continue;

        ll p_cnt = s_cnt[p - 1]; // exact
        ll p_sum = s_sum[p - 1]; // modded
        ll q = p * p;
        // contribution of composites with smallest prime factor = p:
        // p * ( C_p(N/p) - C_p(p-1) )
        ll diff_cnt = l_cnt[p] - p_cnt; // exact
        ans = (ans + p * diff_cnt % MOD + MOD) % MOD;

        ll end = min(v, N / q);

        // update large values: x = N//i, i = 1..end
        for (ll i = 1; i <= end; i++)
        {
            ll d = i * p;
            ll dc;
            ll ds;
            if (d <= v)
            {
                dc = l_cnt[d] - p_cnt;
                ds = l_sum[d] - p_sum;
            }
            else
            {
                ll t = N / d; // t <= v
                dc = s_cnt[t] - p_cnt;
                ds = s_sum[t] - p_sum;
            }
            l_cnt[i] -= dc;
            l_sum[i] = ((l_sum[i] - ds % MOD * p) % MOD + MOD) % MOD;
        }

        // update small values: x from v down to q
        for (ll x = v; x >= q; x--)
        {
            ll t = x / p;
            ll dc = s_cnt[t] - p_cnt;
            ll ds = (s_sum[t] - p_sum + MOD) % MOD;
            s_cnt[x] -= dc;
            s_sum[x] = ((s_sum[x] - p * ds) % MOD + MOD) % MOD;
        }
    }

    // remaining numbers are primes; their smpf equals themselves
    return (ans + l_sum[1]) % MOD;
}

int main()
{
    ios::sync_with_stdio(false);
    cin.tie(nullptr);
    cout << S_min_prime_factor_mod(N) << "\n";
    return 0;
}