求[1,n]中所有素数之和

问题描述

给定整数 $n$，求$[1,n]$中所有素数的和，答案对998244353取模。

数据范围

1s，512M

对于 $100\%$ 的数据，$n\le 10^{10}$

解题思路

回忆素数的欧拉筛法（Euler’s_sieve）。

1. 对于当前数 $x$，若未被判定为合数，即为素数，将其插入素数列表末尾。

2. 将 $x$ 与素数列表中的所有数从小到大相乘，并将乘积判定为合数。若 $x$ 为当前素数的倍数，则不再与下一个素数相乘。

因此当一个数被判定为合数时，必然是通过其最小素因数筛出。

证明如下：

设 $y=p_1^{a_1}\times p_2^{a_2}\times \cdots \times p_n^{a_n}$，其中 $p_1<p_2<\cdots <p_n$，$p_i$ 为素数。

若 $y$ 通过 $p_i$ 筛出，则 $y=p_i\times (p_1^{a_1}\times p_2^{a_2}\times \cdots \times p_i^{a_i-1}\cdots \times p_n^{a_n})$，而 $p_1^{a_1}\times p_2^{a_2}\times \cdots \times p_i^{a_i-1}\cdots \times p_n^{a_n}$ 为 $p_1$ 的倍数，因而在 $p_1$ 筛完后操作终止，不可能继续到 $p_i$。

定义 $S(x,y)$，表示在 $[2,x]$中，通过所有不大于 $y$ 的质数筛完后依然未被判定为合数的数之和。

定于 $T(x,y)$，表示满足 $S(x,y)$ 要求的所有数的集合。

根据上文分析，得到关于 $T(x,y)$ 的性质：属于 $T(x,y)$ 的数要么是素数，要么其最小素因子大于 $y$。

那么问题所求的值为 $S(n,\lfloor \sqrt{n}\rfloor )\mathrm{m}\mathrm{o}\mathrm{d}998244353$。

可以通过状态转移方程求解。

初始值 $S(x,1)=\sum _{i=2}^{x}i=\frac{(x+2)\times (x-1)}{2}$，此时未被任何素数筛选，所有数都符合要求。

若 $y$ 为合数，则 $S(x,y)=S(x,y-1)$。因为素数列表和筛选范围都没有发生变化。

若 $y$ 为素数，且 $y^2>x$，则 $S(x,y)=S(x,y-1)$。因为任何数 $p$ 的因数都不大于 $\sqrt{p}$，而 $y>\sqrt{x}$，因此 $y$ 不可能是 $[2,x]$ 中任何数的最小素因数。

现在讨论最后一种情况，若 $y$ 为素数，且 $y^2\le x$。

考虑 $S(x,y)$ 相比于 $S(x,y-1)$ 多筛去了哪些数，即 $[2,x]$ 中所有最小素因数为 $y$ 的数。

这些数可以表示为 $y\times z_i$，其中 $z_i\in [2,\lfloor \frac{x}{y}\rfloor ]$，且 $z_i$ 的最小素因数不小于 $y$。

而根据 $T(x,y)$ 的性质， $T(\lfloor \frac{x}{y}\rfloor ,y-1)$ 为所有最小素因数大于 $y-1$ 或其本身为素数的数的集合。

因而 $z_i\in T(\lfloor \frac{x}{y}\rfloor ,y-1)$，且 { z i } + { p r i m e ≤ y − 1 } = T ( ⌊ x y ⌋ , y − 1 ) \{z_i\} + \{prime \leq y-1\} = T(\lfloor\frac{x}{y}\rfloor,y-1) {zi​}+{prime≤y−1}=T(⌊yx​⌋,y−1)，而 $\sum _{prime\le y-1}prime$ 可以表示为 $S(y-1,y-1)$。

得到 $S(x,y)=S(x,y-1)-y\times (S(\lfloor \frac{x}{y}\rfloor ,y-1)-S(y-1,y-1))$。

参考代码

C++11:

#include <bits/stdc++.h>
#include <bits/extc++.h>

typedef long long ll;
const ll MOD = 998244353;

ll N;
std::vector<ll> A;
__gnu_pbds::cc_hash_table<ll, ll> S;

int main() {
	scanf("%d", &N);
	ll r = sqrt(N);
	for (ll i = 1; i <= r; ++i)
		A.push_back(N / i);
	for (ll i = A.back() - 1; i >= 1; --i)
		A.push_back(i);
	for (auto i : A) 
		S[i] = (i + 2) % MOD * ((i - 1) % MOD) / 2;
	for (ll p = 2; p <= r; ++p) {
		if (S[p] != S[p - 1]) {
			ll sp = S[p - 1];
			ll pp = p * p;
			for (auto i : A) {
				if (i < pp) break;
				S[i] = ((S[i] - p * (S[i / p] - sp)) % MOD + MOD) % MOD;
			}
		}
	}
	printf("%lld\n", S[N]);
	return 0;
}

python3:

def Solve(n):
    r = int(n ** 0.5)
    V = [n // i for i in range(1, r + 1)]
    V += list(range(V[-1] - 1, 0, -1))
    S = {i : i * (i + 1) // 2 - 1 for i in V}
    tms = 0
    for p in range(2, r + 1):
        if S[p] != S[p - 1]:
            sp = S[p - 1]
            p2 = p * p
            for v in V:
                if v < p2: break
                S[v] -= p * (S[v // p] - sp)
    return S[n]

if __name__ == '__main__':
    print(Solve(int(input())) % 998244353)

注意事项

注意到LLONG_MAX=9223372036854775807，即不超过 $10^{19}$，因而两个 $10^{10}$ 级别的数相乘可能溢出，需要先取模再相乘。

使用cc_hash_table 而不是map 的原因在于前者是哈希表，定位复杂度为 $O(1)$。而通过测试发现同为哈希表的 gp_hash_table 和 unordered_map 的表现都很不尽如人意，耗时至少5倍以上。

复杂度分析

实测大约是 $O(n^{0.75})$ 的量级，但是不太会计算具体值。

参考链接

Project Euler-Problem 10 Lucy_Hedgehog

知乎-求十亿内所有质数的和，怎么做最快? 菜鱼ftfish