$\pi(x)$ 的计算

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

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

理论基础: 参考潘承洞《数论基础》和 ams上的一篇论文另一篇论文

$\psi(x,s)$

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

其中 $m_s = p_1 \cdots p_s$ 为前$s$个素数的积。

另一方面,显然我们有

$\pi(x)$

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

$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$ 此时

其中

注意到上式中 $\frac{x}{p_k} < x^{\frac{2}{3}}$

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

取上面 $s = \pi(\sqrt[3]{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) $
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
#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 = 7;
const int PM = 2*3*5*7*11*13*17;
int phi[PM+1][M+1],sz[M+1];
void init(){
initprime();
pi[2] = 1;
for(int i=3;i<N;++i){
if(isp[i]) pi[i]=pi[i-1]+1;
else pi[i]=pi[i-1];
}
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<N && x/p[s]/p[s] <= p[s]){
int s2x = pi[(int)(sqrt(x+0.2))];
LL ans = pi[x]-LL(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 = pi[int(sqrt(x+0.2))];
int ps3x = pi[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;
}

int main(){
init();
LL n;
while(cin>>n){ // n=98765432109876 = 9.8*10^13 用时 38s, N 大点,耗时会小点
time_t now = time(0);
cout<<primepi(n)<<endl;
cout<<"Time used: "<<difftime(time(0),now)<<endl;
}
return 0;
}

lehmer计算公式

我自己写的代码没有上面的快,两种计算各有优势

令 $s = \pi(\sqrt[4]{x}), t= \pi(\sqrt[3]{x}) $。则,对任意$i>3, P_i(x,s) = 0$,

即:

注意到 $\frac{x}{p_k} < \sqrt{x}$ ,所以最后一个式子可以用下式求,最后计算复杂度在于$P_2(x,s)$

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
#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 = 7;
const int PM = 2*3*5*7*11*13*17;
int phi[PM+1][M+1],sz[M+1];
void init(){
initprime();
pi[2] = 1;
for(int i=3;i<N;++i){
if(isp[i]) pi[i] = pi[i-1]+1;
else pi[i]=pi[i-1];
}
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 lehmerpi(LL x);
LL P_2(LL x, int s);
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 lehmerpi(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 P_2(LL x,int s){
int ps2x = lehmerpi(sqrt(x+0.2));
LL ans = LL(s-ps2x)*(ps2x+s-1)/2;
for(int i=s+1;i<=ps2x;++i){
ans += lehmerpi(x/p[i]);
}
return ans;
}
LL lehmerpi(LL x){
if(x<N) return pi[x];
int ps3x = lehmerpi(int( cbrt(x+0.2) ));
int ps4x = lehmerpi(int( sqrt(sqrt(x+0.2)) ));
LL ans = primephi(x,ps4x) + ps4x-1 - P_2(x,ps4x);
for(int i =ps4x+1;i<=ps3x;++i){
ans -= P_2(x/p[i],i-1);
}
return ans;
}

int main(){
init();
LL n;
while(cin>>n){ // n=98765432109876 = 9.8*10^13 用时 42s, N 大点,耗时会小点
time_t now = time(0);
cout<<lehmerpi(n)<<endl;
cout<<"Time used: "<<difftime(time(0),now)<<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-1) = s-1$,最后一项可以写成 $dp(p_{s-1},s-1)$。虽然上面需要二维数组,但是实际上我们可以优化成一维数组的情况。因为

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

若$L[m]!=L[m-1]$ ,则说明 $m$ 不能被前 $s$ 个素数整除是第$s+1$个素数。

我们需要的是 $R[1]$

上述做法的时间复杂度为 $O(\frac{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
#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
const int N = 1e7+2;
LL L[N],R[N];
LL primepi(LL n){
LL rn = (LL)sqrt(n+0.2);
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];
}
}
return R[1];
}
int main(){
LL n;
while(cin>>n){ // n=98765432109876 = 9.8*10^13 用时 193s
cout<<primepi(n)<<endl;
}
return 0;
}

这个骚方法还目前还没有找到其它的应用0.0

主要还没法对它用树状数组

求第n个素数的方法:

  • 根据概率分布找到大致界限
  • 再牛顿梯度法(或者二分查找)使得 $\pi(m)= n $
  • 素数判断递减$m$ 直到 $m$ 为素数
  • 参考这里

$\psi(x,s)$ 计算太慢了,需要优化:

理论依据见此压缩包

我们知道,若 $\pi(\sqrt[3]{x}) \leq s \leq \pi(\sqrt{x})$ 则

其中

注意到上式中 $\frac{x}{p_k} < x^{\frac{2}{3}} $

我们之前的操作本质是递归让$x,s$ 变小,通过打表预处理来加速递归使得满足一定的效率需要。

现在我们来直接计算得到我们的答案。

符号约定:

取定整数 $\sqrt[3]{x} \leq y = x^{\frac{1}{3}+\epsilon }< \sqrt{x}, z = \frac{x}{y}, s = \pi(y)$ ,约定$p,q$ 是素数。 预处理 $y$ 以内的数组:isp[], p[], mu[], pi[] ,对 $ [1,z]$ 内的 $\psi(x,s)$ 用树状数组(可以在我的博客网站搜索:树状数组)维护(注意到这需要O(z)的内存,单次维护不现实,所以我们可以一段段的维护,保证每一段的长度为$2^{\lfloor \log_2(y) \rfloor +1}$ 来提高效率)

做完上面的预处理后,我们发现 $P_2(x,s)$ 可以直接计算了。

$\psi(x,s)$ 直接计算

在本博文的最开始有:

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

但是最右边本质上有很多项,所以直接算,其实复杂度特别高。

我们还有递归公式:

可以一直分解下去,如果一直分解下去就可以得到最上面的公式了。

所以我们约定:对于节点 $\mu(n) \psi(\frac{x}{n},b)$ ,如果满足

  • 原始节点: $b = c, n \leq y$
  • 特殊节点:$n>y$

就不再分解。这也等价于说 如果 $n < y$ 且 $b>c$ 就分解。因为一开始 $n=1,b=a>c$。而且 $n$ 会增大,$b$ 会减小,所以节点一定会有限步内,落入上述两种框架中的一种。并且 特殊节点的父节点 $-\mu(n) \psi(\frac{x}{\frac{n}{p_{d+1}}},b+1)$ 必然满足 $ \frac{n}{ p_{d+1} } \leq y c$。综上:

以前设置 $c = 7$,但是后来发现没必要,$c=0$ 就挺好。

$S_0$ 很好处理,计算$S$: 对 $p = \delta(n)$ 一起求:

注意到:$\frac{x}{mp} < z$ ,所以我们已经可以把 $\psi(x,s)$ 直接计算出来了。

但是我们可以避免很多计算来提高效率。于是我们有下列一系列的操作

  • $p \geq \sqrt{y}$,则 $m$为素数 ,且此时 $mp > p^2 \geq y $ (若$m$ 为合数,则 $m \geq \delta(m) ^2 >p^2 \geq y$ 矛盾)
  • $\frac{x}{mp} < p$ 时 $\psi(\frac{x}{mp},\pi(p)-1) = 1$
  • $\sqrt{\frac{x}{mp}} < p$ 时 $\psi(\frac{x}{mp},\pi(p)-1) = \pi(\frac{x}{mp}) - \pi(p) +2$
  • $\sqrt{y} < \sqrt{z} < x^{\frac{1}{4}} < x^{\frac{1}{3}} < y$

由此对$S$分段计算:

由限制关系式,我们化简 $S_1, S_2, S_3 $

$S_2$ 也可以用 $S_3$ 的式子求,只是效率不高。

$S_2$ 中 $\frac{x}{pq},p<y$,即可以直接求。

当然了还可以继续细化,但是我嫌麻烦就不想再细化了!

也就是现在的核心就是树状数组分段维护数据,然后每一段的总值要用数组存起来就好了。然后用这个数据结果计算 $S_2,S_3,S_4,P_2(x,s)$,然后就大功告成了0.0

用$\psi(x,s)$ 计算 $\pi(x)$,还是用 $\pi(x)$ 计算 $\psi(x,s)$ ,这是一个问题。

用树状数组维护的时候会有一个很大的问题就是:求和式中 每此动$p$ 整个维护就要从 $1 \to p$ 重新维护一次很麻烦。这个问题没解决所以我不想写代码。

想把 $\frac{x}{pq}$ 的所有可能的值单调递增排列但是又不现实。

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