$\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}))$。它可以利用
- $\psi(x,0) = \lfloor x \rfloor$
- $\psi(x,s) = \psi(x,s-1) - \psi(\frac{x}{p_s},s-1)$
至此问题貌似就这么解决了。但是由于这个递归会使得程序效率大大降低，因此需要一些预处理操作。
- 若 $x<p_s$ 则 $\psi(x,s) ＝ 1$
- 给定一个小整数 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 |
|
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 |
|
稳定简洁的 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
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 = 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}$ 的所有可能的值单调递增排列但是又不现实。
