不超过 x 的素数个数

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

理论基础: 参考潘承洞《数论基础》以及论文包.zip

\(\psi(x,s)\)

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

\[\psi(x,s) = \sum_{n \leq x} \sum_{d|(n,m_s)} \mu(d) = \sum_{d|m_s} u(d)\lfloor \frac{x}{d} \rfloor \]

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

另一方面，显然我们有

\[\psi(x,s) = \psi(x,s-1) - \psi(\frac{x}{p_s},s-1) \]

\(\pi(x)\)

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

\[\psi(x,s) = \pi(x) - s + 1 \]

\(P_k(x,s)\)

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

\[\psi(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\) 此时

\[\psi(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) \]

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

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

\[\pi(x) = \psi(x,s) + \frac{(\pi(\sqrt{x})+s-2)(\pi(\sqrt{x})- s+1)}{2} - \sum_{k=s+1}^{\pi(\sqrt{x})} \pi(\frac{x}{p_k}) \]

取上面 $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)\)

lehmer 计算公式

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

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

\[\begin{array}{rl} \psi(x,s) &= 1 + \pi(x) - s + P_2(x,s) + P_3(x,s) \\ &= 1+ \pi(x) - s + P_2(x,s) + \sum_{k=s+1}^{t} P_2(\frac{x}{p_k},k-1) \\ \end{array} \]

即：

\[\pi(x) = \psi(x,s)-1+s-P_2(x,s) - \sum_{k=s+1}^{t} P_2(\frac{x}{p_k},k-1) \]

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

\[\sum_{k=s+1}^{t} P_2(\frac{x}{p_k},k-1) = \sum_{k=s+1}^{t} \sum_{j=k}^{\pi(\sqrt{\frac{x}{p_k}})} \pi(\frac{x}{p_k p_j}) - j+1 \]

稳定简洁的 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(x,s) = dp(x,s-1)-dp(\frac{x}{p_s},s-1)+s-1 \]

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

\[dp(x,s) = dp(x,s-1)-dp(\frac{x}{p_s},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})\) 且常数特别小，代码十分简洁。

#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})\) 则

\[\psi(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) = \sum_{k=s+1}^{\pi(\sqrt{x})} \left( \psi(\frac{x}{p_k},s)+ s - k \right) \]

注意到上式中 $\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)\) 直接计算

在本博文的最开始有：

\[\psi(x,s) = \sum_{n \leq x} \sum_{d|(n,p)} \mu(d) = \sum_{d|p} u(d)\lfloor \frac{x}{d} \rfloor \]

其中 \(p\) 为前 \(s\) 个素数的积。

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

我们还有递归公式：

\[\begin{array}{cl} \psi(x,s) &= \psi(x,s-1) - \psi(\frac{x}{p_s},s-1) \\ &= \psi(x,s-2) - \psi(\frac{x}{p_{s-1}},s-2) - \psi(\frac{x}{p_s},s-2) + \psi(\frac{x}{p_sp_{s-1}},s-2) \end{array} \]

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

所以我们约定：对于节点 \(\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 <n\) 且 \(b+1>c\)。综上：

以前设置 \(c = 7\)，但是后来发现没必要，\(c=0\) 就挺好。

\[\psi(x,s) = \sum_{n=1} ^y \mu(n) \lfloor \frac{x}{n} \rfloor + \sum_{\frac{n}{\delta(n)} \leq y < n} \mu(n) \psi(\frac{x}{n}, \pi(\delta(n))-1) = S_0 + S \]

\(S_0\) 很好处理，计算 \(S\)： 对 \(p = \delta(n)\) 一起求：

\[S = - \sum_{p \leq y} \sum_{\begin{array}{c} \delta(m) >p \\ m \leq y < mp \end{array}} \mu(m) \psi(\frac{x}{mp},\pi(p)-1) \]

注意到：\(\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 = - \sum_{\sqrt[3]{x} < p \leq y} \sum_{\begin{array}{c} \delta(m) >p \\ m \leq y < mp \end{array}} \mu(m) \psi(\frac{x}{mp},\pi(p)-1) \]

\[S_2 = - \sum_{\sqrt{z} < p \leq \sqrt[3]{x}} \sum_{\begin{array}{c} \delta(m) >p \\ m \leq y < mp \end{array}} \mu(m) \psi(\frac{x}{mp},\pi(p)-1) \]

\[S_3 = - \sum_{\sqrt{y} < p \leq \sqrt{z}} \sum_{\begin{array}{c} \delta(m) >p \\ m \leq y < mp \end{array}} \mu(m) \psi(\frac{x}{mp},\pi(p)-1) \]

\[S_4 = - \sum_{p \leq \sqrt{y}} \sum_{\begin{array}{c} \delta(m) >p \\ m \leq y < mp \end{array}} \mu(m) \psi(\frac{x}{mp},\pi(p)-1) \]

由限制关系式，我们化简 \(S_1, S_2, S_3\)

\[\begin{array}{cl} S_1 &= - \sum_{\sqrt[3]{x} < p \leq y} \sum_{p<q \leq y} \mu(q) \psi(\frac{x}{pq},\pi(p)-1) \\ &= \sum_{\sqrt[3]{x} < p < q \leq y} 1 \\ &= {\pi(y)-\pi(\sqrt[3]{x}) \choose 2} \end{array} \]

\[\begin{array}{cl} S_2 &= - \sum_{\sqrt{z} < p \leq \sqrt[3]{x}} \sum_{p<q \leq y} \mu(q) \psi(\frac{x}{pq},\pi(p)-1) \\ &= \sum_{\sqrt{z} < p \leq \sqrt[3]{x}} \sum_{p<q \leq y} \pi(\frac{x}{pq}) - \pi(p) +2 \\ &= \sum_{\sqrt{z} < p \leq \sqrt[3]{x}} \sum_{p<q \leq y} \pi(\frac{x}{pq}) + \sum_{\sqrt{z} < p \leq \sqrt[3]{x}} (\pi(p)-\pi(y))(\pi(p)-2) \end{array} \]

\[\begin{array}{cl} S_3 &= - \sum_{\sqrt{y} < p \leq \sqrt{z}} \sum_{p < q \leq y} \mu(q) \psi(\frac{x}{pq},\pi(p)-1) \\ &= \sum_{\sqrt{y} < p \leq \sqrt{z}} \sum_{p < q \leq y} \psi(\frac{x}{pq},\pi(p)-1) \end{array} \]

\(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}\) 的所有可能的值单调递增排列但是又不现实。

标签：psi,frac,超过,个数,sqrt,leq,素数,pi,sum
来源： https://www.cnblogs.com/izlyforever/p/computationOfPiX.html