标签：psi frac 超过 个数 sqrt leq 素数 pi sum
\(\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\)。则
显然 \(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}))\)。它可以利用
- \(\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)\)
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) \]\[\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 \]注意到 \(\frac{x}{p_k} < \sqrt{x}\) ，所以最后一个式子可以用下式求，最后计算复杂度在于 \(P_2(x,s)\)
稳定简洁的 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\)。综上：
\[\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 \]以前设置 \(c = 7\)，但是后来发现没必要，\(c=0\) 就挺好。
\(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
本站声明： 1. iCode9 技术分享网（下文简称本站）提供的所有内容，仅供技术学习、探讨和分享； 2. 关于本站的所有留言、评论、转载及引用，纯属内容发起人的个人观点，与本站观点和立场无关； 3. 关于本站的所有言论和文字，纯属内容发起人的个人观点，与本站观点和立场无关； 4. 本站文章均是网友提供，不完全保证技术分享内容的完整性、准确性、时效性、风险性和版权归属；如您发现该文章侵犯了您的权益，可联系我们第一时间进行删除； 5. 本站为非盈利性的个人网站，所有内容不会用来进行牟利，也不会利用任何形式的广告来间接获益，纯粹是为了广大技术爱好者提供技术内容和技术思想的分享性交流网站。