自动搬运

查看原文

来自洛谷，原作者为

	ケロシ Blue Archive

搬运于2025-08-24 22:33:57，当前版本为作者最后更新于2024-07-10 09:24:19，作者可能在搬运后再次修改，您可在原文处查看最新版

自动搬运只会搬运当前题目点赞数最高的题解，您可前往洛谷题解查看更多

以下是正文

---

一道 Meissel-Lehmer 算法的模板题洛谷题解里没有 Meissel-Lehmer，大悲啊。

这里介绍 Meissel-Lehmer 算法。

先线性筛算出前大于等于 $n^\frac{1}{2}$ 个数内的素数和 $\pi(n)$。

然后对于求 $\pi(n)$，考虑对其产生贡献的素数分治。 考虑对于 $n^\frac{1}{3}$ 分治，先算出 $(n^\frac{1}{3},n]$ 中对于 $n^\frac{1}{3}$ 以内的素数筛不掉的数，暂定一个 dp，$dp_{i,j}$ 表示 $[1,i]$ 中不会被前 $j$ 个素数筛掉的个数。

接下来不难发现，现在筛剩下的 $(n^\frac{1}{3},n]$ 中的数，如果是合数，则其必定为两个大于 $n^\frac{1}{3}$ 的素数的乘积，其他即为素数。

计算这类由两个素数相乘的数的个数很简单，直接考虑枚举小的素数，去计算另外一个素数的方案数即可，假设当前选取了小素数 $p$，则其筛掉的数就是 $[p,\frac{n}{p}]$ 中的素数的个数，即为 $\pi(\frac{n}{p})-\pi(p)+1$。

结合一下上述，可得出 $\pi(n)$ 的递归公式：

$$\pi(n)=dp_{n,\pi(n^\frac{1}{3})}-1+\pi(n^\frac{1}{3}) -\sum_{p\in(n^\frac{1}{3},n^\frac{1}{2}]} \pi(\frac{n}{p})-\pi(p)+1 $$

公式前面减一是因为 $1$ 是筛不掉的，其中 $p$ 为素数。

然后考虑如何求 dp 数组，考虑转移。考虑计算最后一个素数的贡献，发现其贡献是其倍数中没有被前面素数筛掉的数的个数，发现可以再利用 dp 数组，所以转移为 $dp_{i,j}=dp_{i,j-1}-dp_{\frac{i}{p_j},j-1}$，其中 $p_j$ 表示第 $j$ 个素数。

考虑边界情况，显然的是 $dp_{i,0}=i$。但这样还不够优，考虑 $p_j^2\ge i$ 的情况，这时 $[1,i]$ 中的合数全被筛过，只剩下素数，所以值即为 $\max(0,\pi(i)-j+1)$，加一是为了留下 $1$，与前面同步。

这样一个基础的 Meissel-Lehmer 算法原理讲解就完成了。

接下来是实现，对于 dp 数组，可以预处理出一部分，例如 $i\le 1.8 \times 10^6,j\le 60$ 的情况。然后由于这道题空间有限，前面公式里的 $\pi(\frac{n}{p})$ 会超出预处理的范围，这是要继续递归。

由于不断递归了 $\pi(n)$，这份代码的时间复杂度我分析不出来，但是时间复杂度应该在 $O(n^\frac{2}{3})$ 级别，足以通过。

const int N = 8e6 + 10;
const int MI = 1.8e6;
const int MJ = 60;
ll n;
int f[MI + 5][MJ + 5];
int p[N / 10], ip[N], g[N], cnt;
void init() {
	FOR(i, 2, N - 1) {
		if(!ip[i]) {
			p[++cnt] = i;
		}
		for(int j = 1; j <= cnt && p[j] * i < N; j++) {
			ip[i * p[j]] = 1;
			if(i % p[j] == 0) {
				break;
			}
		}
	}
	// 警惕！这里的 g 数组表示的 pi(n) 含有 1，是 pi(n)+1
	FOR(i, 1, N - 1) g[i] = g[i - 1] + !ip[i];
	FOR(i, 1, MI) f[i][0] = i;
	FOR(i, 1, MI) FOR(j, 1, MJ)
		f[i][j] = f[i][j - 1] - f[i / p[j]][j - 1];
}
ll dp(ll i, ll j) {
	if(i <= MI && j <= MJ) 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);
}
ll pi(ll n) {
	if(n < N) return g[n] - 1;
	ll sn = pow(n, 1. / 3);
	ll m = g[sn] - 1;
	ll res = dp(n, m) + m - 1;
	for(m++; 1ll * p[m] * p[m] <= n; m++)
		res -= pi(n / p[m]) - pi(p[m]) + 1; // 继续递归
	return res;
}