题解：P7884 【模板】Meissel-Lehmer - 指南

P7884

简介

Meissel-Lehmer 是一种素数筛，适用于 $n$ 超大的范围。

时间复杂度对比

算法名称	预处理时间复杂度	单次查询时间复杂度	空间复杂度	适用场景
Meissel-Lehmer	$O(n^{2/3}/\log n)$	$O(1)$	$O(n^{2/3})$	超大范围素数计数（$n\ge 10^{12}$）
分段筛（Segment）	无	$O(n\log \log n)$	$O(\sqrt{n})$	中等范围素数生成（$n\le 10^8$）
欧拉筛（Euler）	$O(n)$	$O(1)$	$O(n)$	小范围素数预处理（$n\le 10^7$）

前置知识

P3383 欧拉筛

欧拉筛可以以 $O(N)$ 的方式求出区间 $[1,N]$ 内的素数。

定义

• $p_x$ 表示第 $x$ 个素数。
• $\pi (x)$ 表示 $[1,x]$ 内的素数个数。
• $d{p}_{x,k}$ 表示 $[1,x]$ 中不被前 $k$ 个素数整除的数的个数。（包含 $1$）

定理推导

定理一：Legendre 公式（基础拆分）

对于任意满足 $k=\pi (\sqrt{n})$ 的 $n$ 和 $k$，有：$\pi (n)=k+(d{p}_{n,k}-1)$

证明

一：所有素数 $p_i$ 满足 ($p_k<p_i\le n$) 的数包含于 $d{p}_{n,k}$

因为 $p_i$ 为素数，所以不会被前 $k$ 个素数整除，满足 $d{p}_{n,k}$ 的定义，得证。

二：$d{p}_{n,k}$ 中除 $1$ 外，其余数均满足 ($p_k<i\le n$)

考虑反证：设 $m$ 存在于符合 $d{p}_{n,k}$ 的定义且不满足 ($p_k<m\le n$)。
因为 $m$ 是合数，存在正整数 $a,b>1$ 满足 $m=a\times b$，因为 $k=\pi (\sqrt{n})$，所以有 $a\times b=m⟹\min (a,b)\le \sqrt{n}$。
又因为 $a,b>1$，所以 $\min (a,b)=p_j$（$j\le k$），与假设矛盾，故得证。

三：最终证明

由第二步得：$d{p}_{n,k}$ 中除 $1$ 外，其余数均满足 ($p_k<p_i\le n$) 。
设 $C$ 为大于 $p_k$ 的素数数量和。
显然有：

• $C=d{p}_{n,k}-1$
• $\pi (n)=k+C$ 等量代换得：$\pi (n)=k+d{p}_{n,k}-1$
  得证。

定理二：如何计算 $d{p}_{n,k}$

可以通过容斥原理与递归计算。
我们有：
不被前 $k$ 个素数整除的数 $=$ 不被前 $k-1$ 个素数整除的数 $-$ 被 $p_k$ 整除但不被前 $k-1$ 个素数整除的数。
而“被 $p_k$ 整除但不被前 $k-1$ 个素数整除的数” 等价于“$p_k$ 乘以不被前 $k-1$ 个素数整除的数（且 $\le x/p_k$）”，故数量为 $d{p}_{x/p_k,k-1}$。
故有：

$$dp(x,k)=\{\begin{array}{ll}x& (k=0)\\ 1& (p_k>x)\\ d{p}_{x,k}=\pi (x)-k+1& (p_k^2>x)\\ d{p}_{x,k}=d{p}_{x,k-1}-d{p}_{\lfloor x/p_k\rfloor ,k-1}& \end{array}$$

证明

对于任意 $x,k\ge 1$，有 $d{p}_{x,k}=d{p}_{x,k-1}-d{p}_{\lfloor x/p_k\rfloor ,k-1}$
设：

• $S=\{1\le m\le x|m\text{不被前}k-1\text{个素数整除}\}，\text{则}|S|=d{p}_{x,k-1}$。
• $T=\{1\le m\le x|m\text{不被前}k\text{个素数整除}\}，\text{则}|T|=d{p}_{x,k}$。

第一步：分析 $S$ 与 $T$ 的关系

前 $k$ 个素数是前 $k-1$ 个素数加 $p_k$，因此“不被前 $k$ 个素数整除” 等价于“不被前 $k-1$ 个素数整除，且不被 $p_k$ 整除”。
有：$|T|=|S|-|\{m\in S|m\text{被}{p}_k\text{整除}\}|$

第二步：证明 $|\{m\in S|m\text{被}{p}_k\text{整除}\}|=d{p}_{\lfloor x/p_k\rfloor ,k-1}$

若 $m\in S$ 且 $m$ 被 $p_k$ 整除，则存在整数 $t$ 使得 $m=p_k\times t$。
因为 $m\le x$ 所以 $t\le \lfloor x/p_k\rfloor$；又因为 $m\in S$，所以不被前 $k-1$ 素数整除，得 $t$ 不被前 $k-1$ 个素数整除。（因为 $p_k$ 与前 $k-1$ 个素数互质，假设 $t$ 被 $p_i(i<k)$ 整除，则 $m=p_k\times t$ 也被 $p_i$ 整除，与 $m\in S$ 矛盾，故假设不成立）。
因为 $t\le \lfloor x/p_k\rfloor$ 且 $t$ 不被前 $k-1$ 个素数整除，则 $m=p_k\times t\le x$，且 $m$ 不被前 $k-1$ 个素数整除，得出 $m\in S$ 且 $m$ 被 $p_k$ 整除。
因此 $S$ 中被 $p_k$ 整除的数与 $1\le t\le \lfloor x/p_k\rfloor$ 且 $t$ 不被前 $k-1$ 个素数整除的数是映射的，也就是两者的个数相等。
因为 $t$ 不被前 $k-1$ 个素数整除的数为 $d{p}_{x,k}$，所以：$|\{m\in S|m\text{被}{p}_k\text{整除}\}|=d{p}_{\lfloor x/p_k\rfloor ,k-1}$
故得证。

---

根据定义，有：
$d{p}_{x,k}=|S|-|\{m\in S|m\text{被}{p}_k\text{整除}\}|$
故得：$d{p}_{x,k}=d{p}_{x,k-1}-d{p}_{\lfloor x/p_k\rfloor ,k-1}$
公式得证。

定理三：最终公式

设 $Q(n,k)$ 表示为大于 $k$ 的两个素数的乘积小于 $n$ 的个数。
当 $n$ 很大时，直接计算 $d{p}_{x,k}$ 是很慢的，所以引入阈值 $a=\pi (n^{1/3})$，可以将 $\pi (n)$ 分为三部分：

$$\pi (n)=\pi (a)+(d{p}_{n,a}-1)-Q(n,a)$$

接下来介绍各部分的含义：

• $\pi (a)$：$[1,a]$ 中的素数个数。
• $d{p}_{n,a}-1$：
  • 不被前 $a$ 个素数整除的数的个数减 1，包含大于 $a$ 的素数和 $Q(n,a)$
• $-Q(n,a)$：
  • 通过递归计算 $d{p}_{n,a}-1-(\pi (n)-a)$ 得到（$d{p}_{n,a}-1$ 是不被前 $a$ 个素数整除的数的个数，减去 $\pi (n)-a$ 个素数，剩余的为 $Q(n,a)$）。$Q(n,a)={\Sigma }_{i=a+1}^{\pi (\sqrt{n})}(\pi (n/p_i)-i+1)$

证明

首先明确命题：
设 $a=\pi (n^{1/3})$，$b=\pi (\sqrt{n})$，则

$$\pi (n)=\pi (a)+(d{p}_{n,a}-1)-Q(n,a)$$

我们首先对 $\pi (n)$ 进行拆分：

• 第一类：对于 $p_i\le a$ 的素数，个数为 $\pi (a)$。
• 第二类：对于 $a<p_i\le n$ 的素数，设其个数为 $D$，则 $\pi (n)=\pi (a)+D$，故需要证明 $D=(d{p}_{n,a}-1)-Q(n,a)$。

---

现在来分析第二类素数的性质：
对于类 $2$ 中的素数 $q(q>a)$，考虑任意合数 $m=q\times r(r\ge q$，保证 $q$ 是 $m$ 的最小因子 $)$，则 $m\le n⟹r\le \lfloor n/q\rfloor$。
由于 $q>a$，得到 $q>\pi (n^{1/3})$，因 $p_{\pi (x)}>x$，故 $m=q\times r\ge q^2>n^{2/3}$，且 $r\le \lfloor n/q\rfloor <n^{2/3}$
我们继续细分第二类质数对应的合数贡献：

• 当 $r\le a$ 时，此时 $m=q\times r$ 的最小质因子为 $r$，对应的，$q$ 需要满足 $q\le \lfloor n/r\rfloor$ 且 $q>r$。
• 当 $r>a$ 时，此时 $m=q\times r$ 的最小质因子为 $q$，这类合数即为 $Q(n,a)$。

---

接下来分别计算两种情况的贡献

第一种

考虑 $a<p_i\le b$，有 $p_i\le n^{1/2}$（因为 $b=\pi (\sqrt{n})$）：

• 对于每个 $p_i(a<i\le b)$，对应的 $q$ 需满足 $q>p_i$ 且 $q\le \lfloor n/p_i\rfloor$
• 这类 $q$ 的个数是 $\pi (\lfloor n/p_i\rfloor )-i+1$。
  因为当 $q=p_i$ 时，$m=p_i^2$，但 $p_i>a=\pi (n^{1/3})>n^{1/3}$，故 $p_i^2>n^{2/3}$，但 $p_i\le n^{1/2}$ 等价于 $p_i^2\le n$，此时 $p_i^2\le n$，此时 $q=p_i$ 也需计入，所以实际个数要加 $1$。
  要对所以的 $i=a+1$ 到 $b$ 求和，所以第一种对于的素数个数为 ${\Sigma }_{i=a+1}^b(\pi (\lfloor n/p_i\rfloor )-i+1)$。

第二种

定义：$m=q\times r$ 是“双素数乘积”。
同时，第二类的素数 $q$ 满足“不是双素数乘积” 与“不被前 $a$ 个素数整除”，结合 $d{p}_{n,a}$ 的定义，有：

$$d{p}_{n,a}-1=Q(n,a)+\text{第二类的素数个数}$$

又因为（$d{p}_{n,a}-1$ 是剔除 $1$ 后的不被前 $a$ 个素数整除的数，包含第二类素数和双素数乘积），所以有第二类素数个数 $=(d{p}_{n,a}-1)-Q(n,a)$。
结合第一种情况，第二类素数个数 $D$ 等于第一种个数 $+$ 第二种对应的质数贡献，所以有：

$$D={\Sigma }_{i=a+1}^b(\pi (\lfloor n/p_i\rfloor )-i+1)-Q(n,a)$$

故得证。

实现思路

预处理小范围素数，再递归实现 $dp$ 与 $\pi (n)$。

明确定义

变量

• $p_i$ 表示第 $i$ 个素数。
• $i{p}_i$ 表是第 $i$ 个数是否是素数。
• $g_i$ 表示 $[1,i]$ 中有多少个素数。
• $f_{i,j}$ 表示小于 $i$ 且不被前 $j$ 个素数整除的个数（也就是 $d{p}_{i,j}$ 的预处理版本）。

定义常数部分

• $N$ 表示欧拉筛处理的范围上限。
• $MX$ 表示 $f_{i,j}$ 的第一维的大小。
• $MY$ 表示 $f_{i,j}$ 的第二维的大小。

一：欧拉筛

这个可以先套一下模板，但同时我们要维护 $g$ 与 $f$ 数组。
对于 $g$ 数组，采用类似前缀和的方式维护：

for(int i=1;i<N;i++)g[i]=g[i-1]+!ip[i];

对于 $f$ 数组，可以采用之前 $dp$ 的定义，有：

for(int i=1;i<=MX;i++)f[i][0]=i;
for(int i=1;i<=MX;i++)for(int j=1;j<=MY;j++)f[i][j]=f[i][j-1]-f[i/p[j]][j-1];

剩下的都是模板了，完整代码：

void shai(){
for(int i=2;i<N;i++){
if(!ip[i])p[++tot]=i;
for(int j=1;j<=tot&&p[j]*i<N;j++){
ip[i*p[j]]=1;
if(i%p[j]==0)break;
}
}
for(int i=1;i<N;i++)g[i]=g[i-1]+!ip[i];
for(int i=1;i<=MX;i++)f[i][0]=i;
for(int i=1;i<=MX;i++)for(int j=1;j<=MY;j++)f[i][j]=f[i][j-1]-f[i/p[j]][j-1];
}

$d{p}_{i,j}$ 的计算

直接储存 $d{p}_{i,j}$ 是存不下的，我们可以考虑预处理小部分。
对于 $d{p}_{i,j}$，有四种情况：

情况	返回值	判定条件	解释
1	$f_{i,j}$	$i\le MX\text{且}j\le MY$	在预处理范围内。
2	$i$	$!i\text{或}!j$	到达边界。
3	$\max \{0,g_i-j\}$	$i<N\text{且}1\times p_j\times p_j\ge i$	剩余数均为素数。
4	$d{p}_{i,j-1}-d{p}_{i/p_j,j-1}$	不满足上述条件。	上文 $d{p}_{i,j}$ 的计算公式。

对于特殊情况的解释：当 $p_j^2\ge i$ 时，小于 $i$ 且不被前 $j$ 个素数整除的数只能是素数（或 $1$），而 $g_i$ 是小于 $i$ 的素数总数，减去前 $j$ 个素数（已被排除），得到结果。
基于此，得出以下代码：

long long dp(long long i,long long j){
if(i<=MX&&j<=MY)return f[i][j];
if(!i||!j)return i;
if(i<N&&
1ll*p[j]*p[j]>=i)return max(0ll,g[i]-j);
return dp(i,j-1)-dp(i/p[j],j-1);
}

$P(N)$ 的计算 (其实就是计算 $\pi (n)$)

当然，如果在预处理范围内就直接返回了。
定义 $sum=d{p}_{n,m}+m-1$，也就是 $d{p}_{n,m}$ 是小于 $n$ 且不被前 $m$ 个素数整除的数，加 $m-1$ 是补上前 $m-1$ 个素数（因前 $m$ 个素数中排除了第 $m$ 个），再减去多计算的合数，代码：

long long P(long long n){
if(n<N)return g[n]-1;
long long m=g[(long long)(pow(n,1./3))]-1,sum=dp(n,m)+m-1;
for(m++;
1ll*p[m]*p[m]<=n;m++)sum-=P(n/p[m])-P(p[m])+1;
return sum;
}

完整代码

#include<bits/stdc++.h>
  using namespace std;
  const int N=8e6+10,MX=1.8e6,MY=60;
  long long n;
  int f[MX+5][MY+5],p[N/10],ip[N],g[N],tot;
  void shai(){
  for(int i=2;i<N;i++){
  if(!ip[i])p[++tot]=i;
  for(int j=1;j<=tot&&p[j]*i<N;j++){
  ip[i*p[j]]=1;
  if(i%p[j]==0)break;
  }
  }
  for(int i=1;i<N;i++)g[i]=g[i-1]+!ip[i];
  for(int i=1;i<=MX;i++)f[i][0]=i;
  for(int i=1;i<=MX;i++)for(int j=1;j<=MY;j++)f[i][j]=f[i][j-1]-f[i/p[j]][j-1];
  }
  long long dp(long long i,long long j){
  if(i<=MX&&j<=MY)return f[i][j];
  if(!i||!j)return i;
  if(i<N&&
  1ll*p[j]*p[j]>=i)return max(0ll,g[i]-j);
  return dp(i,j-1)-dp(i/p[j],j-1);
  }
  long long P(long long n){
  if(n<N)return g[n]-1;
  long long m=g[(long long)(pow(n,1./3))]-1,sum=dp(n,m)+m-1;
  for(m++;
  1ll*p[m]*p[m]<=n;m++)sum-=P(n/p[m])-P(p[m])+1;
  return sum;
  }
  int main(){
  shai();
  long long n;
  cin>>n;
  cout<<
  P(n);
  return 0;
  }

后记

学了好久啊！
锐评我的代码：先放个记录，可以看到，代码的内存刚好卡线（评测机波动分分钟给你 MLE 了），所以请不要直接交（真有那样的人吗），可以尝试改小 $N$ 和 $MX$ 和 $MY$ 的范围再交。（或者拼运气？）
这是本人第一次写黑题题解，难免会有错误，希望大家指正。
最后，感谢观看！