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

问题描述

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

数据范围

1s，512M

对于 100%100\%100% 的数据，n≤1010n\leq 10^{10}n≤1010

解题思路

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

对于当前数 xxx，若未被判定为合数，即为素数，将其插入素数列表末尾。
将 xxx 与素数列表中的所有数从小到大相乘，并将乘积判定为合数。若 xxx 为当前素数的倍数，则不再与下一个素数相乘。

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

证明如下：

设 y=p1a1×p2a2×⋯×pnany=p_1^{a_1}\times p_2^{a_2}\times\cdots\times p_n^{a_n}y=p1a1×p2a2×⋯×pnan，其中 p1<p2<⋯<pnp_1 < p_2 < \cdots < p_np1<p2<⋯<pn，pip_ipi 为素数。

若 yyy 通过 pip_ipi 筛出，则 y=pi×(p1a1×p2a2×⋯×piai−1⋯×pnan)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})y=pi×(p1a1×p2a2×⋯×piai−1⋯×pnan)，而 p1a1×p2a2×⋯×piai−1⋯×pnanp_1^{a_1}\times p_2^{a_2}\times\cdots\times p_i^{a_i-1}\cdots \times p_n^{a_n}p1a1×p2a2×⋯×piai−1⋯×pnan 为 p1p_1p1 的倍数，因而在 p1p_1p1 筛完后操作终止，不可能继续到 pip_ipi。

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

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

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

那么问题所求的值为 S(n,⌊n⌋)mod998244353S(n, \lfloor\sqrt{n}\rfloor) \bmod 998244353S(n,⌊n⌋)mod998244353。

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

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

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

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

现在讨论最后一种情况，若 yyy 为素数，且 y2≤xy^2 \leq xy2≤x。

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

这些数可以表示为 y×ziy\times z_iy×zi，其中 zi∈[2,⌊xy⌋]z_i\in [2,\lfloor\frac{x}{y}\rfloor]zi∈[2,⌊yx⌋]，且 ziz_izi 的最小素因数不小于 yyy。

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

因而 zi∈T(⌊xy⌋,y−1)z_i\in T(\lfloor\frac{x}{y}\rfloor,y-1)zi∈T(⌊yx⌋,y−1)，且 {zi}+{prime≤y−1}=T(⌊xy⌋,y−1)\{z_i\} + \{prime \leq y-1\} = T(\lfloor\frac{x}{y}\rfloor,y-1){zi}+{prime≤y−1}=T(⌊yx⌋,y−1)，而 ∑prime≤y−1prime\sum\limits_{prime \leq y-1}{prime}prime≤y−1∑prime 可以表示为 S(y−1,y−1)S(y-1,y-1)S(y−1,y−1)。

得到 S(x,y)=S(x,y−1)−y×(S(⌊xy⌋,y−1)−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))S(x,y)=S(x,y−1)−y×(S(⌊yx⌋,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 = 0for p in range(2, r + 1):if S[p] != S[p - 1]:sp = S[p - 1]p2 = p * pfor v in V:if v < p2: breakS[v] -= p * (S[v // p] - sp)return S[n]if __name__ == '__main__':print(Solve(int(input())) % 998244353)

注意事项

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

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

复杂度分析

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

参考链接

Project Euler-Problem 10 Lucy_Hedgehog

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