不超过$x$的素数之和

我们可以在不求出不超过$x$的所有素数的情况下,求出最终结果。

$f_k(x,s)$: 最小素因子大于$p_s$ 且不超过$x$ 的数的$k$次方和。

其中, $\delta(m)$ 表示$m$的最小素因子(约定 $\delta(1) = + \infty$)。

$f_k(x,s)$ 的递推公式:

$f_k(x,0) = \sum_{i=1} ^ {\lfloor x \rfloor} i^k$

$f_k(x,s)$ 和 $f_k(x)$ 的关系:

若 $s >= \pi(\sqrt{x})$,则

只对$k=1$写代码(一般的$k$ 同理)

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
#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
const int N = 1e7+2;
int p[N],pi[N];
bool isp[N];
int initprime(){
int cnt=1;p[1]=2;isp[2] = true;
for(int i=3;i<N;i+=2) isp[i]=true;
for(int i=3;i<N;i+=2){
if(isp[i]) p[++cnt] = i;
for(int j = 2,t=(N-1)/i + 1;j <= cnt && p[j] < t; ++j){
isp[i * p[j]] = false;
if(i%p[j] == 0) break;
}
}
return cnt;
}
LL spi[N];
void init(){
initprime();
pi[2] = 1;spi[2]=2;
for(int i=3;i<N;++i){
if(isp[i]){
pi[i]=pi[i-1]+1;
spi[i] = spi[i-1]+i;
}else{
pi[i]=pi[i-1];
spi[i] = spi[i-1];
}
}
}
LL sump(LL x);
LL sumphi(LL x, int s){
if(s == 0) return (x+1)*x/2;
if(x/p[s] <= p[s]) return sump(x)+1-sump(p[s]);
return sumphi(x,s-1)-p[s]*sumphi(x/p[s],s-1);
}
LL sump(LL x){
if(x<N) return spi[x];
int s = pi[int(sqrt(x+0.2))];
return sumphi(x,s) -1+ sump(p[s]);
}

int main(){
init();
LL n;
while(cin>>n){ // n = 10^10 就超过int64了,此时用Int128 n < N*N
time_t now = time(0);
cout<<sump(n)<<endl;
cout<<"Time used: "<<difftime(time(0),now)<<endl;
}
return 0;
}

但是其实我们还可以,类似于 $\pi(x)$ 的计算一样, 取 $\pi(\sqrt[3]{x}) \leq s \leq \pi(\sqrt{x})$

只对$k=1$写代码(一般的$k$ 同理)

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
#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
const int N = 1e7+2;
int p[N],pi[N];
bool isp[N];
int initprime(){
int cnt=1;p[1]=2;isp[2] = true;
for(int i=3;i<N;i+=2) isp[i]=true;
for(int i=3;i<N;i+=2){
if(isp[i]) p[++cnt] = i;
for(int j = 2,t=(N-1)/i + 1;j <= cnt && p[j] < t; ++j){
isp[i * p[j]] = false;
if(i%p[j] == 0) break;
}
}
return cnt;
}
LL spi[N];
void init(){
initprime();
pi[2] = 1;spi[2]=2;
for(int i=3;i<N;++i){
if(isp[i]){
pi[i]=pi[i-1]+1;
spi[i] = spi[i-1]+i;
}else{
pi[i]=pi[i-1];
spi[i] = spi[i-1];
}
}
}
LL sump(LL x);
LL sumphi(LL x, int s){
if(s == 0) return (x+1)*x/2;
if(x/p[s] <= p[s]) return sump(x)+1-sump(p[s]);
if(x<N && x/p[s]/p[s] <= p[s]){
int s2x = pi[(int)(sqrt(x+0.2))];
LL ans = spi[x]+1-spi[p[s]];
for(int i=s+1;i<=s2x;++i){
ans+=p[i]*(spi[x/p[i]]-spi[p[i-1]]);
}
return ans;
}
return sumphi(x,s-1)-p[s]*sumphi(x/p[s],s-1);
}
LL sump(LL x){
if(x<N) return spi[x];
int ps2x = pi[int(sqrt(x+0.2))];
int ps3x = pi[int(cbrt(x+0.2))];
LL ans = sumphi(x,ps3x) -1 + spi[p[ps3x]];
for(int i =ps3x+1;i<=ps2x;++i){
ans -= p[i]*(sump(x/p[i])-sump(p[i-1]));
}
return ans;
}

int main(){
init();
LL n;
while(cin>>n){ // n = 10^10 就超过int64了,此时用Int128 n < N*N
time_t now = time(0);
cout<<sump(n)<<endl;
cout<<"Time used: "<<difftime(time(0),now)<<endl;
}
return 0;
}
欧尼酱,人家想喝可乐!