P7884

简介

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

时间复杂度对比

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

前置知识

P3383 欧拉筛

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

定义

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

定理推导

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

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

证明

一：所有素数 $p_{i}$ 满足 ($p_{k}<p_{i}\le n$) 的数包含于 $dp_{n,k}$

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

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

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

三：最终证明

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

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

定理二：如何计算 $dp_{n,k}$

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

$$dp(x,k)=\begin{cases} x&(k=0)\\ 1&(p_k>x)\\ dp_{x,k}=\pi(x)-k+1&(p_k^2>x)\\ dp_{x,k}=dp_{x,k-1}-dp_{\lfloor x/p_k\rfloor,k-1} \end{cases}$$

证明

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

• $S=\{1≤m≤x|m\text{ 不被前 }k−1\text{ 个素数整除 }\}，\text{则 }|S|=dp_{x,k−1}$。
• $T=\{1≤m≤x|m\text{ 不被前 }k\text{ 个素数整除}\}，\text{则 }|T|=dp_{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{ 整除} \}|=dp_{\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$ 个素数整除的数为 $dp_{x,k}$，所以：$|\{m\in S|m \text{ 被 }p_k\text{ 整除} \}|=dp_{\lfloor x/p_k\rfloor,k-1}$
故得证。

---

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

定理三：最终公式

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

$$\pi(n)=\pi(a) + (dp_{n,a}-1) - Q(n,a)$$

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

• $\pi(a)$：$[1,a]$ 中的素数个数。
• $dp_{n,a}-1$：
  • 不被前 $a$ 个素数整除的数的个数减 1，包含大于 $a$ 的素数和 $Q(n,a)$
• $-Q(n,a)$：
  • 通过递归计算 $dp_{n,a}-1-(\pi(n)-a)$ 得到（$dp_{n,a}-1$ 是不被前 $a$ 个素数整除的数的个数，减去 $\pi(n)-a$ 个素数，剩余的为 $Q(n,a)$）。$Q(n,a) = \varSigma_{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) + (dp_{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=(dp_{n,a}-1)-Q(n,a)$。

---

现在来分析第二类素数的性质：
对于类 $2$ 中的素数 $q(q>a)$，考虑任意合数 $m=q\times r(r\ge q$，保证 $q$ 是 $m$ 的最小因子 $)$，则 $m\le n\implies 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$ 求和，所以第一种对于的素数个数为 $\varSigma_{i=a+1}^{b}(\pi(\lfloor n/p_i\rfloor)-i+1)$。

第二种

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

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

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

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

故得证。

实现思路

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

明确定义

变量

• $p_i$ 表示第 $i$ 个素数。
• $ip_i$ 表是第 $i$ 个数是否是素数。
• $g_i$ 表示 $[1,i]$ 中有多少个素数。
• $f_{i,j}$ 表示小于 $i$ 且不被前 $j$ 个素数整除的个数（也就是 $dp_{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];
}

$dp_{i,j}$ 的计算

直接储存 $dp_{i,j}$ 是存不下的，我们可以考虑预处理小部分。
对于 $dp_{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	$dp_{i,j-1}-dp_{i/p_j,j-1}$	不满足上述条件。	上文 $dp_{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=dp_{n,m}+m-1$，也就是 $dp_{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;
}

后记

最后，感谢观看！