$\pi(x)$ 的计算

发表于 |

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

1. $\phi(x,s)$

设 $p_i$ 为第 $i$ 个素数，$\phi(x,s)$ 表示不超过 $x$ 且能不能被 $p_i,(1 \leq i \leq s)$ 整除的正整数个数。即
$$ \phi(x,s) = \sum_{n \leq x} \sum_{d|(n,p)} \mu(d) = \sum_{d|p} u(d)\lfloor \frac{x}{d} \rfloor $$
其中 $p = p_1 \cdots p_s$。

另一方面，显然我们有
$$ \phi(x,s) = \phi(x,s-1) - \phi(\frac{x}{p_s},s-1) $$

2. $\pi(x)$

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

3. $P_k(x,s)$

设 $P_k(x,s)$ 为不大于 $x$ 且正好有 $k$ 个大于 $p_s$ 的素因子(按重根计)的整数个数（方法属于 Lehmer）。进一步设 $P_0(x,s)=1$。则
$$ \phi(x,s) = \sum_{k=0} ^{\infty} P_k(x,s) $$
显然 $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$ 此时
$$ \phi(x,s) = 1 + \pi(x)-s + P_2(x,s) $$
其中
$$ P_2(x,s) = \sum_{k=s+1}^{\pi(\sqrt{x})} \left( \pi(\frac{x}{p_k}) - k + 1 \right) $$
即
$$ \phi(x,s) = \pi(x) + \sum_{k=s+1}^{\pi(\sqrt{x})} \pi(\frac{x}{p_k}) - \frac{(\pi(\sqrt{x})+s-2)(\pi(\sqrt{x})-s+1)}{2} $$

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

$$ \pi(x) = \phi(x,\pi(\sqrt[3]{x})) - \sum_{k=\pi(\sqrt[3]{x})+1}^{\pi(\sqrt{x})} \pi(\frac{x}{p_k}) + \frac{(\pi(\sqrt{x})+\pi(\sqrt[3]{x})-2)(\pi(\sqrt{x})-\pi(\sqrt[3]{x})+1)}{2} $$
因此问题最终转化成求 $\phi(x,\pi(\sqrt[3]{x}))$。它可以利用

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

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

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

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

令 $s = \pi(\sqrt[4]{x}),t=\pi(\sqrt[2]{x}),m=\pi(\sqrt[3]{x})$。则
$$ \begin{align}
\phi(x,s) &= 1 + \pi(x) - s + P_2(x,s) + P_3(x,s) \\
&= 1+ \pi(x) - s + \sum_{k=s+1}^{t} \left( \pi(\frac{x}{p_k}) - k +1 \right) + \sum_{k=s+1}^{m} P_2(\frac{x}{p_k},k-1) \\
&= \pi(x) - \frac{(t-s+1)(t+s-2)}{2} + \sum_{k=s+1}^{t} \pi(\frac{x}{p_k}) \sum_{k=s+1}^{m} P_2(\frac{x}{p_k},k-1)
\end{align}
$$
即
$$ \pi(x) = \phi(x,s) + \frac{(t-s+1)(t+s-2)}{2} -\sum_{k=s+1}^{t} \pi(\frac{x}{p_k}) - \sum_{k=s+1}^{m} P_2(\frac{x}{p_k},k-1)
$$

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) = \phi(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(x,s) = dp(x,s-1)-dp(\frac{x}{p_s},s-1)+s-1 $$
因为 $dp(p_{s-1},s) = s-1$,最后一项可以写成 $dp(p_{s-1},s)$。虽然上面需要二维数组，但是实际上我们可以优化成一维数组的情绪。因为
$$ dp(x,s) = dp(x,s-1)-dp(\frac{x}{p_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 getans(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<<getans(n)<<endl;
    }
    return 0;
}