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

fk(x)pxpk

fk(x,s): 最小素因子大于 ps 且不超过 x 的数的 k 次方和

fk(x,s)=mx,δ(m)>psmk

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

fk(x,s) 的递推公式

fk(x,s)=mx,δ(m)>psmk=f(x,s1)mx,δ(m)=psmk=f(x,s1)pskf(xps,s1)

fk(x,0)=i=1xik

fk(x,s)fk(x) 的关系

s>=π(x),则

fk(x)=fk(x,s)1+fk(ps)

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

c++
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;
}

但是其实我们还可以,类似于 π(x) 的计算一样, 取 π(x3)sπ(x)

fk(x)=fk(x,s)1+fk(ps)i=s+1π(x)(f(xpi)f(pi1))pik

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

c++
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;
}