题解：P7884 【模板】Meissel-Lehmer
ILearnedSomeCoding
·
2026-01-02 17:56:31
·
题解
我们打开这个算法的 OI wiki……啊，好复杂。看不懂。
这里公式是 AI 生成的。
Meissel-Lehmer 算法求 \pi(x)
定义 \phi(x,a) 为所有不超过 x 的正整数中，其所有质因子都大于 p_a 的数的个数，即
\phi(x,a)
=\#\{\,n\le x \mid n\bmod p=0 \Rightarrow p>p_a\,\}
再定义 P_k(x,a) 表示所有不超过 x 的正整数中，可重质因子恰好有 k 个，且所有质因子都大于 p_a 的数的个数，即
P_k(x,a)
=\#\{\,n\le x \mid n=q_1q_2\cdots q_k \Rightarrow \forall i,\ q_i>p_a\,\}
特殊地，定义
P_0(x,a)=1
于是有
\phi(x,a)=P_0(x,a)+P_1(x,a)+\cdots+P_k(x,a)+\cdots
该无穷和实际上是有限和，因为当 p_a^k>x 时，有 P_k(x,a)=0 。
设 y 为满足 x^{1/3}\le y\le x^{1/2} 的整数，并记
a=\pi(y)
当 k\ge3 时，有
P_1(x,a)=\pi(x)-a,\qquad P_k(x,a)=0
由此可得
\pi(x)=\phi(x,a)+a-1-P_2(x,a)
因此，计算 \pi(x) 可以转化为计算 \phi(x,a) 与 P_2(x,a) 。
后面其实就是计算 P_2 用的一堆公式，但是没必要优化，因为不优化也能过。
接下来只需要计算 P_2(x,a) 。
根据定义，P_2(x,a) 表示不超过 x 的正整数中，恰好由两个质数相乘，且这两个质数都大于 p_a 的数的个数。由于取 a=\pi(y) 且 x^{1/3}\le y\le x^{1/2} ，可以发现：
若 p\le x^{1/3} ，则 p\le p_a ，不符合条件；p>x^{1/2} ，则 pq>x ，不可能成立。
因此只需枚举较小的质因子 p ，其中
x^{1/3}<p\le\sqrt{x}.
对每个这样的质数 p ，满足 pq\le x 的质数 q 的个数为
\pi\!\left(\frac{x}{p}\right)-\pi(p)+1，
其中减去 \pi(p) 是为了保证 q\ge p ，避免重复计数。
于是有
P_2(x,a)
=
\sum_{x^{1/3}<p\le\sqrt{x}}
\left(
\pi\left(\frac{x}{p}\right)-\pi(p)+1
\right).
按定义，\phi(x,a) 表示 1..x 中不被前 a 个质数 2,3,5,\dots,p_a 整除的数的个数（等价于其所有质因子都大于 p_a ）。
计算 \phi 的经典递推为：
\phi(x,a)=\phi(x,a-1)-\phi\!\left(\left\lfloor\frac{x}{p_a}\right\rfloor,a-1\right).
这部分数形如 $p_a\cdot t$，且 $t\le\left\lfloor x/p_a\right\rfloor$，并且 $t$ 也不能含有前 $a-1$ 个质数作为因子，所以因此数量恰好是 $\phi(\lfloor x/p_a\rfloor,a-1)$。
为了加速，在较小范围内可以用 DP 预处理表格：
令 $s[n][j]=\phi(n,j)$，则
$$
s[n][0]=n,
$$
$$
s[n][j]=s[n][j-1]-s\!\left(\left\lfloor\frac{n}{p_j}\right\rfloor,j-1\right).
$$
算的方法是
- 先赋值 $s[i][0]=i$；
- 再按 $j=1\dots70$ 递推：
`s[i][j] = s[i][j-1] - s[i/prime[j]][j-1]`。
当需要 $\phi(x,a)$ 且 $x$ 在预处理范围内时可直接返回 `s[x][a]`，否则再用递归计算。
将其代入
$$
\pi(x)=\phi(x,a)+a-1-P_2(x,a)
$$
就可以算出来答案了。
时间复杂度：
前面预处理贼慢，后面是 $O(n^{\frac{2}{3}})$，后面递归算的其实影响不大，因为顶多递归两次。最后耗时 $3.5$ 秒钟运行~~只能说评测机太强了，自己测试要至少五秒钟~~。
代码：
```
// ILearnedSomeCoding 的万能头1.1
#include<bits/stdc++.h>
using namespace std;
#define ll long long
#define int64 long long
#define ull unsigned long long
#define uint64 unsigned long long
#define int128 __int128
#define uint128 __uint128_t
#define PQ priority_queue
#define fir(_S, _E, _T) for(int i = _S;i <= _E;i += _T)
#define fjr(_S, _E, _T) for(int j = _S;j <= _E;j += _T)
int s[1234567][78];
ll prime[12345678], ppi[12345678], cnt;
constexpr ll lim = 10000000;
constexpr ll lim2 = 1200000;
bool isprime[12345678];
void init(){
fir(2, lim, 1){
if(!isprime[i]){
prime[++cnt] = i;
}
for(int j = 1; j <= cnt && (ll)prime[j] * i <= lim; j++){
isprime[i * prime[j]] = 1;
if(i % prime[j] == 0) break;
}
}
isprime[1] = 1;
ppi[0] = 0;
fir(1, lim, 1) ppi[i] = ppi[i - 1] + !isprime[i];
fir(1, lim2, 1) s[i][0] = i;
fir(1, lim2, 1){
fjr(1, 70, 1){
s[i][j] = s[i][j - 1] - s[i / prime[j]][j - 1];
}
}
}
ll isqrt(ll n){
long double x = sqrtl((long double)n);
ll k = (ll)x;
if(k < 0) k = 0;
while((k + 1) * (k + 1) <= n) k++;
while(k * k > n) k--;
return k;
}
ll icbrt(ll n){
long double x = cbrt((long double)n);
ll k = (ll)x;
if(k < 0) k = 0;
while((k + 1) * (k + 1) * (k + 1) <= n) k++;
while(k * k * k > n) k--;
return k;
}
ll pi(ll n);
ll phi(ll a, ll b){
if(a < lim2 && b < 50) return s[a][b];
if(a == 0 || b == 0) return a;
if(a < lim && prime[b] * prime[b] >= a)
return max(0ll, ppi[a] - (ll)b + 1);
return phi(a, b - 1) - phi(a / prime[b], b - 1);
}
ll P(ll n){
ll k = icbrt(n);
ll m = ppi[k];
ll s2 = 0;
for(int i = m + 1; (ll)prime[i] * prime[i] <= n; i++)
s2 += pi(n / prime[i]) - pi(prime[i]) + 1;
return s2;
}
ll pi(ll n){
if(n < lim) return ppi[n];
ll k = icbrt(n);
ll a = ppi[k];
return phi(n, a) - 1 + a - P(n);
}
int main(){
init();
ll n;
cin >> n;
cout << pi(n);
return 0;
}
```