第$n$个素数

很早之前写过 $\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) $ 的网站

世界纪录保持者的求法

如有帮助,烦请资瓷(一块也是爱0.0)