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

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

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

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

$\displaystyle 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;
}