很早之前写过 $\pi(x)$ 的计算 ，在知乎上用它回答问题的时候，发现我怎么没有写 求第 $n$ 个素数。

做法依赖于 $\pi(x)$ 的计算，$\pi(x)$ 表示不超过 $x$ 的素数个数

素数定理( 这里就不证了)

从而我们知道：

其中，$p_n$ 为第 $n$ 个素数，显然 $p_n$ 是 $\pi(x) = n$ 最小的解。

$p_n$ 求解

  • 预处理小于 $N$ 的素数
  • 初始值 $n\ln n$
  • 牛顿梯度法
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
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
const int N = 1e8+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;
}
const int M = 8; // 不能再大了不然内存顶不住了
const int PM = 2*3*5*7*11*13*17*19;
int phi[PM+1][M+1],sz[M+1];
void init(){
    initprime();
    pi[2] = 1;
    for(int i=3;i<N;++i){
        pi[i]=pi[i-1];
        if(isp[i])    ++pi[i];
    }
    sz[0]=1;
    for(int i=0;i<=PM;++i)  phi[i][0]=i;
    for(int i=1;i<=M;++i){
        sz[i]=p[i]*sz[i-1];
        for(int j=1;j<=PM;++j){
            phi[j][i]=phi[j][i-1]-phi[j/p[i]][i-1];
        }
    }
}
LL primepi(LL x);
LL primephi(LL x, int s){
    if(s <= M)  return phi[x%sz[s]][s]+(x/sz[s])*phi[sz[s]][s];
    if(x/p[s] <= p[s])  return primepi(x)-s+1;
    if(x/p[s]/p[s] <= p[s] && x<N){
        int s2x = pi[(int)(sqrt(x+0.2))];
        LL ans = pi[x]-(s2x+s-2)*(s2x-s+1)/2;
        for(int i=s+1;i<=s2x;++i){
            ans+=pi[x/p[i]];
        }
        return ans;
    }
    return primephi(x,s-1)-primephi(x/p[s],s-1);
}
LL primepi(LL x){
    if(x<N) return pi[x];
    int ps2x = primepi(int(sqrt(x+0.2)));
    int ps3x = primepi(int(cbrt(x+0.2)));
    LL ans = primephi(x,ps3x) + LL(ps2x+ps3x-2)*(ps2x-ps3x+1)/2;
    for(int i =ps3x+1,ed = ps2x;i<=ed;++i){
        ans -= primepi(x/p[i]);
    }
    return ans;
}
bool isprime(LL n){  // 可以用概率判别来替换
    if(n<N) return isp[n];
    for(int i=1;p[i]<=n/p[i]&&i<N;++i){ 
        if(p[i]&&n%p[i]==0) return false;
    }
    return true;
}
LL primep(LL n){    // Newton 梯度法
    if(n <= pi[N-1])  return p[n];
    LL ans= n*log(n), err = log(n)/log(10);
    LL m = primepi(ans);
    while(m<n||m>n+err){
        ans += (n-m)/(log(m)-1)*log(m)*log(m);
        m = primepi(ans);
    }
    int step = m-n;
    for(int i=0;i<=step;++i){
        while(!isprime(ans))    --ans;
        --ans;
    }
    return ++ans;
}

int main(){
    init();
    LL n;   // n < N*N 
    while(cin>>n){ // n=987654321098 = 9.8*10^11 用时 128s 太慢了。
        cout<<primep(n)<<endl;
    }
    return 0;
}

第 $n$ 个素数和 $\pi(x)$ 的网站

世界纪录保持者的求法