$\pi(x)$ 的计算

$\pi(x)$ 表示不超过 $x$ 的素数个数。容易看出可以在 $O(N)$时间复杂度,$O(N)$ 空间复杂度离线预处理求出小于 $N$ 的素数全体。但是如果 $N=10^{12}$ 或者更大,这种做法必然是不现实的。因此下面给出高效的求解方法…

再读一次还是觉得很巧妙,一步一步走向更优。

1. $\psi(x,s)$

$\psi(x,s)$ 表示不超过 $x$ 且能不能被前$s$个素数整除的正整数个数。即

其中 $p$ 为前$s$个素数的积。

另一方面,显然我们有

2. $\pi(x)$

我们知道一个数 $n>1$ 是素数当且仅当不存在素数 $p \leq \sqrt{n}$ 使得 $p \mid n$。因此当 $s \geq \pi(\sqrt{x})$ 时,

3. $P_k(x,s)$

设 $P_k(x,s)$ 为不超过 $x$ 且每个素因子都大于$p_s$ 且素因子(按重根计)个数为$k$ 的整数个数(方法属于 Lehmer)。进一步设 $P_0(x,s)=1$。则

显然 $P_1(x,s) = \pi(x)-s$。

若 $\pi(\sqrt[3]{x}) \leq s \leq \pi(\sqrt{x})$ 则 $P_k(x,s)=0,k \geq 3$ 此时

其中

4. $\pi(x)$ 的计算公式

因此问题最终转化成求 $\psi(x,\pi(\sqrt[3]{x}))$。它可以利用

  1. $\psi(x,0) = \lfloor x \rfloor$
  2. $\psi(x,s) = \psi(x,s-1) - \psi(\frac{x}{p_s},s-1)$

至此问题貌似就这么解决了。但是由于这个递归会使得程序效率大大降低,因此需要一些预处理操作。

  1. 若 $x<p_s$ 则 $\psi(x,s) = 1$
  2. 给定一个小整数M,预处理出 $\psi(x,s)$,其中 $x < q=p_1 \cdots p_s,\quad s<=M$
    则 $\psi(x,s) = \psi(x \mod q,s) + \lfloor \frac{x}{q} \rfloor \psi(q,s) $

5. 更为实用的lehmer计算公式

令 $s = \pi(\sqrt[4]{x}),t=\pi(\sqrt[2]{x}),m=\pi(\sqrt[3]{x})$。则

6. 例题:Codeforce 665F

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
93
94
95
const int N = 5e6+2;
bool np[N];
int p[N],pi[N];
int getprime(){
int cnt=0;
np[0]=np[1]=true;
pi[0]=pi[1]=0;
for(int i = 2; i < N; ++i){
if(!np[i]) p[++cnt] = i;
for(int j = 1;j <= cnt && i * p[j] < N; ++j) {
np[i * p[j]] = true;
if(i % p[j] == 0) break;
}
pi[i]=cnt;
}
return cnt;
}
const int M = 7;
const int PM = 2*3*5*7*11*13*17;
int phi[PM+1][M+1],sz[M+1];
void init(){
getprime();
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];
}
}
}
int sqrt2(LL x){
LL r = (LL)sqrt(x-0.1);
while(r*r<=x) ++r;
return int(r-1);
}
int sqrt3(LL x){
LL r = (LL)cbrt(x-0.1);
while(r*r*r<=x) ++r;
return int(r-1);
}
LL getphi(LL x,int s){
if(s == 0) return x;
if(s <= M) return phi[x%sz[s]][s]+(x/sz[s])*phi[sz[s]][s];
if(x <= p[s]*p[s]) return pi[x]-s+1;
if(x <= p[s]*p[s]*p[s] && x< N){
int s2x = pi[sqrt2(x)];
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 getphi(x,s-1)-getphi(x/p[s],s-1);
}
LL getpi(LL x){
if(x < N) return pi[x];
LL ans = getphi(x,pi[sqrt3(x)])+pi[sqrt3(x)]-1;
for(int i=pi[sqrt3(x)]+1,ed=pi[sqrt2(x)];i<=ed;++i){
ans -= getpi(x/p[i])-i+1;
}
return ans;
}
LL lehmer_pi(LL x){
if(x < N) return pi[x];
int a = (int)lehmer_pi(sqrt2(sqrt2(x)));
int b = (int)lehmer_pi(sqrt2(x));
int c = (int)lehmer_pi(sqrt3(x));
LL sum = getphi(x, a) + LL(b + a - 2) * (b - a + 1) / 2;
for (int i = a + 1; i <= b; i++) {
LL w = x / p[i];
sum -= lehmer_pi(w);
if (i > c) continue;
LL lim = lehmer_pi(sqrt2(w));
for (int j = i; j <= lim; j++) {
sum -= lehmer_pi(w / p[j]) - (j - 1);
}
}
return sum;
}
LL getans(LL x){ // x < 1e11
LL ans = pi[sqrt3(x)];
for(int i=1,ed=pi[sqrt2(x-1)];i<=ed;++i){
ans += lehmer_pi(x/p[i])-i;
}
return ans;
}
int main(){
init();
LL n;
while(cin>>n){
cout<<getans(n)<<endl;
}
return 0;
}

简洁的DP做法。

我们令 $dp(x,s) = \psi(x,s)-s+1$ 它的意义是,$2~x$ 中被前$s$ 个素数筛完后的伪素数个数。因此我们有 $dp(0,0)=0,dp(x,0)=x-1,x>1,dp(x,\pi(\sqrt{x})) = \pi(x)$ 且有状态转移

因为 $dp(p_{s-1},s) = s-1$,最后一项可以写成 $dp(p_{s-1},s)$。虽然上面需要二维数组,但是实际上我们可以优化成一维数组的情况。因为

另外我们也不可能开 $O(n)$ 的数组,但是可以利用一种黑科技开 $O(\sqrt{n})$ 的数组即可达到我们的目的。
即我们用 $L[x]$ 表示 $dp(x,s)$ 用 $R[x]$ 表示 $dp(\frac{n}{x},s)$。

但是上述做法的复杂度为 $\frac{O(n)}{\log 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
LL sqrt2(LL x){
LL r = (LL)sqrt(x-0.1);
while(r*r<=x) ++r;
return r-1;
}
LL sqrt3(LL x){
LL r = (LL)cbrt(x-0.1);
while(r*r*r<=x) ++r;
return r-1;
}
const int N = 1e6+2;
LL L[N],R[N];
LL pi(LL n){ // n < 1e12
LL rn = sqrt2(n);
for(LL i=1;i<=rn;++i) R[i]=n/i-1;
LL ln = n/(rn+1)+1;
for(LL i=1;i<=ln;++i) L[i]=i-1;
for(LL p=2;p<=rn;++p){
if(L[p]==L[p-1]) continue;
for(LL i=1,tn=min(n/(p*p),rn);i<=tn;++i){
R[i] -= (i*p<=rn?R[i*p]:L[n/(i*p)])-L[p-1];
}
for(LL i=ln;i>=p*p;--i){
L[i] -= L[i/p]-L[p-1];
}
}
LL ans = L[sqrt3(n)];
for(LL p=2;p<=rn;++p){
if(L[p] == L[p-1]) continue;
ans += R[p]-L[p];
}
return ans;
}
int main(){
LL n;
while(cin>>n){
cout<<pi(n)<<endl;
}
return 0;
}

欧尼酱,人家想喝可乐!