$\pi(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 = () $ 因此问题最终转化成求 \(\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})\) 且常数特别小，代码十分简洁。

这个骚方法还目前还没有找到其它的应用 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) \]

注意到上式中 $ < x^{} $

我们之前的操作本质是递归让 \(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}\) 的所有可能的值单调递增排列但是又不现实。

第 n 个素数

做法依赖于 \(\pi(x)\) 的计算

素数定理( 这里就不证了)

\[ \lim _{x \to \infty} \frac{\pi(x)}{x/\ln x} = 1 \]

从而我们知道： \[ \lim _{x \to \infty} \frac{p_n}{n \ln n} = 1 \] 其中，\(p_n\) 为第 \(n\) 个素数，显然 \(p_n\) 是 \(\pi(x) = n\) 最小的解。

\(p_n\) 求解

• 预处理小于 \(N\) 的素数
• 初始值 \(n\ln n\)
• 牛顿梯度法

第 \(n\) 个素数和 \(\pi(x)\) 的网站

世界纪录保持者的求法

\(\sum_{p \leq n, p \text{ is prime}} p\)

我们可以在不求出 不超过 \(x\) 的所有素数 的情况下，求出最终结果。

\[ f_k(x) \doteq \sum_{p \leq x} p^k \]

• \(f_0(x) = \pi(x)\)
• \(f_1(x) = \sum_{p \leq x} p\) 这是我们关心的结果
• 对于一般的 \(k\) 借助自然数方幂和快速算法也可以求

\(f_k(x,s)\)：最小素因子大于 \(p_s\) 且不超过 \(x\) 的数的 \(k\) 次方和

\[ f_k(x,s) = \sum_{m \leq x, \delta(m)>p_s} m^k \] 其中，\(\delta(m)\) 表示 \(m\) 的最小素因子（约定 \(\delta(1) = + \infty\)）。

\(f_k(x,s)\) 的递推公式

\[ \begin{aligned} f_k(x,s) &= \sum_{m \leq x, \delta(m)>p_s} m^k \\ &= f(x,s-1) - \sum_{m \leq x, \delta(m) = p_s} m^k \\ &= f(x,s-1) - p_s ^k f(\frac{x}{p_s},s-1) \end{aligned} \]

\(\displaystyle f_k(x,0) = \sum_{i=1} ^ {\lfloor x \rfloor} i^k\)

\(f_k(x,s)\) 和 \(f_k(x)\) 的关系

若 \(s >= \pi(\sqrt{x})\)，则 \[ f_k(x) = f_k(x,s) - 1 + f_k(p_s) \]

其实我们还可以，类似于 \(\pi(x)\) 的计算一样, 取 \(\pi(\sqrt[3]{x}) \leq s \leq \pi(\sqrt{x})\) \[ f_k(x) = f_k(x,s) - 1 + f_k(p_s) - \sum_{i=s+1} ^ {\pi(\sqrt{x})} (f(\frac{x}{p_i})-f(p_{i-1})) p_i ^k \]