# 【算法学习笔记】Meissel-Lehmer 算法 （亚线性时间找出素数个数）

「Meissel-Lehmer 算法」是一种能在亚线性时间复杂度内求出 $1\sim n$ 内质数个数的一种算法。

## 代码实现

#include <bits/stdc++.h>
using namespace std;
using ll = long long;
//通过知道前面的 n^1/3 的质数可以推断后面n^2/3的质数所以可以适当减小
const int N  = 9e3;
const int M  = 2;         //为了减小内存可以不过是质数
const int PM = 2 * 3 * 5; //为了减小内存可以不过要按质数减小如去掉17
ll n;
bool np[N];
int prime[N], pi[N];
int phi[PM + 1][M + 1], sz[M + 1];

int getprime() {
int cnt = 0;
np[0] = np[1] = true;
pi[0] = pi[1] = 0;
for (int i = 2; i < N; ++i) {
if (!np[i]) prime[++cnt] = i;
pi[i] = cnt;
for (int j = 1; j <= cnt && i * prime[j] < N; ++j) {
np[i * prime[j]] = true;
if (i % prime[j] == 0) break;
}
}
return cnt;
}

void init() {
getprime();
sz[0] = 1;
for (int i = 0; i <= PM; ++i) phi[i][0] = i;
for (int i = 1; i <= M; ++i) {
sz[i] = prime[i] * sz[i - 1];
for (int j = 1; j <= PM; ++j) phi[j][i] = phi[j][i - 1] - phi[j / prime[i]][i - 1];
}
}

int sqrt2(ll x) {
ll r = (ll)sqrt(x - 0.1);
while (r * r <= x) ++r;
return int(r - 1);
}

int sqrt3(ll x) {
ll r = (ll)cbrt(x - 0.1);
while (r * r * r <= x) ++r;
return int(r - 1);
}

ll getphi(ll x, int s) {
if (s == 0) return x;
if (s <= M) return phi[x % sz[s]][s] + (x / sz[s]) * phi[sz[s]][s];
if (x <= prime[s] * prime[s]) return pi[x] - s + 1;
if (x <= prime[s] * prime[s] * prime[s] && x < N) {
int s2x = pi[sqrt2(x)];
ll ans  = pi[x] - (s2x + s - 2) * (s2x - s + 1) / 2;
for (int i = s + 1; i <= s2x; ++i) ans += pi[x / prime[i]];
return ans;
}
return getphi(x, s - 1) - getphi(x / prime[s], s - 1);
}

ll getpi(ll x) {
if (x < N) return pi[x];
ll ans = getphi(x, pi[sqrt3(x)]) + pi[sqrt3(x)] - 1;
for (int i = pi[sqrt3(x)] + 1, ed = pi[sqrt2(x)]; i <= ed; ++i) ans -= getpi(x / prime[i]) - i + 1;
return ans;
}

ll lehmer_pi(ll x) { //小于等于n的素数有多少个
if (x < N) return pi[x];
int a  = (int)lehmer_pi(sqrt2(sqrt2(x)));
int b  = (int)lehmer_pi(sqrt2(x));
int c  = (int)lehmer_pi(sqrt3(x));
ll sum = getphi(x, a) + (ll)(b + a - 2) * (b - a + 1) / 2;
for (int i = a + 1; i <= b; i++) {
ll w = x / prime[i];
sum -= lehmer_pi(w);
if (i > c) continue;
ll lim = lehmer_pi(sqrt2(w));
for (int j = i; j <= lim; j++) sum -= lehmer_pi(w / prime[j]) - (j - 1);
}
return sum;
}

int main() {
ios_base::sync_with_stdio(false), cin.tie(0);
init();
while (cin >> n && n) cout << lehmer_pi(n) << "\n";
return 0;
}

OJ 上的几道相关的题：Here

## 记号规定

$\left[x\right]$ 表示对 $x$ 下取整得到的结果。
$p_k$ 表示第 $k$ 个质数，$p_1=2$
$\pi\left(x\right)$ 表示 $1\sim x$ 范围内素数的个数。
$\mu\left(x\right)$ 表示莫比乌斯函数。

$\delta\left(x\right)$ 表示 $x$ 最小的质因子。
$P^+\left(n\right)$ 表示 $x$ 最大的质因子。

## Meissel-Lehmer 算法求 π(x)

$\phi\left(x,a\right)=\#\big\{n\le x\mid n\bmod p=0\Rightarrow p>p_a\big\}\tag{1}$

$P_k\left(x,a\right)=\#\big\{n\le x\mid n=q_1q_2\cdots q_k\Rightarrow\forall i,q_i>p_a\big\}\tag{2}$

$\phi\left(x,a\right)=P_0\left(x,a\right)+P_1\left(x,a\right)+\cdots+P_k\left(x,a\right)+\cdots$

$y$ 为满足 $x^{1/3}\le y\le x^{1/2}$ 的整数，再记 $a=\pi\left(y\right)$

$k\ge 3$ 时，有 $P_1\left(x,a\right)=\pi\left(x\right)-a$$P_k\left(x,a\right)=0$，由此我们可以推出：

$\pi\left(x\right)=\phi\left(x,a\right)+a-1-P_2\left(x,a\right)\tag{3}$

## 计算 P₂(x,a)

$P_2\left(x,a\right)=\sum\limits_{y<p\le \sqrt{x}}{\left(\pi\left(\dfrac{x}{p}\right)-\pi\left(p\right)+1\right)}\tag{4}$

$p\in \left[y+1,\sqrt{x}\right]$ 时，我们有 $\dfrac{x}{p}\in \left[1,\dfrac{x}{y}\right]$。因此，我们可以筛区间 $\left[1,\dfrac{x}{y}\right]$，然后对于所有的的质数 $p\in \left[y+1,\sqrt{x}\right]$ 计算 $\pi\left(\dfrac{x}{p}\right)-\pi\left(p\right)+1$。为了减少上述算法的空间复杂度，我们可以考虑分块，块长为 $L$。若块长 $L=y$，则我们可以在 $O\left(\dfrac{x}{y}\log{\log{x}}\right)$ 的时间复杂度，$O\left(y\right)$ 的空间复杂度内计算 $P_2\left(x,a\right)$

## 计算 ϕ(x,a)

1. 可以被 $p_b$ 整除的；
2. 不可以被 $p_b$ 整除的。

$\phi\left(u,0\right)=\left[u\right]\tag{5}$

$\phi\left(x,b\right)=\phi\left(x,b-1\right)-\phi\left(\dfrac{x}{p_b},b-1\right)\tag{6}$

$\phi\left(x,a\right)=\sum\limits_{1\le n\le x\\\\ P^+\left(n\right)\le y}{\mu\left(n\right)\left[x/n\right]}$

$\begin{matrix}&&\phi\left(x,a\right)&&\\ &\swarrow&&\searrow&\\ &\phi\left(x,a-1\right)&&-\phi\left(\frac{x}{p_a},a-1\right)&\\ \swarrow&\downarrow&&\downarrow&\searrow\\ \phi\left(x,a-2\right)&\phi\left(\frac{x}{p_{a-1}},a-2\right)&&-\phi\left(\frac{x}{p_a},a-2\right)&\phi\left(\frac{x}{p_ap_{a-1}},a-2\right)\end{matrix}\\ \vdots\\$

1. $b=0$$n\le y$
2. $n>y$

1. 如果叶子节点 $\mu\left(n\right)\phi\left(\dfrac xn,b\right)$ 满足 $n\le y$，则称这种叶子节点为 普通叶子
2. 如果叶子节点 $\mu\left(n\right)\phi\left(\dfrac xn,b\right)$ 满足 $n>y$$n=mp_b\left(m\le y\right)$，则称这种节点为 特殊叶子

$\phi\left(x,a\right)=S_0+S\tag{7}$

$S_0=\sum\limits_{n\le y}{\mu\left(n\right)\left[\dfrac xn\right]}\tag{8}$

$S$ 表示 特殊叶子 的贡献：

$S=\sum\limits_{n/\delta\left(n\right)\le y\le n}{\mu\left(n\right)\phi\left(\dfrac{x}{n},\pi\left(\delta\left(n\right)\right)-1 \right)}\tag{9}$

## 计算 S

$S=-\sum\limits_{p\le y}{\ \sum\limits_{\delta\left(m\right)>p\\\\ m\le y<mp}{\mu\left(m\right)\phi\left(\dfrac{x}{mp},\pi\left(p\right)-1\right)}}\tag{10}$

$S=S_1+S_2+S_3$

$S_1=-\sum\limits_{x^{1/3}<p\le y}{\ \sum\limits_{\delta\left(m\right)>p\\m\le y<mp}{\mu\left(m\right)\phi\left(\dfrac{x}{mp},\pi\left(p\right)-1\right)}}$

$S_2=-\sum\limits_{x^{1/4}<p\le x^{1/3}}{\ \sum\limits_{\delta\left(m\right)>p\\\\ m\le y<mp}{\mu\left(m\right)\phi\left(\dfrac{x}{mp},\pi\left(p\right)-1\right)}}$

$S_3=-\sum\limits_{p\le x^{1/4}}{\ \sum\limits_{\delta\left(m\right)>p\\\\ m\le y<mp}{\mu\left(m\right)\phi\left(\dfrac{x}{mp},\pi\left(p\right)-1\right)}}$

$S_1=\sum\limits_{x^{1/3}<p\le y}{\ \sum\limits_{p<q\le y}{\phi\left(\dfrac{x}{pq},\pi\left(p\right)-1\right)}}$

$S_2=\sum\limits_{x^{1/4}<p\le x^{1/3}}{\ \sum\limits_{p<q\le y}{\phi\left(\dfrac{x}{pq},\pi\left(p\right)-1\right)}}$

### 计算 S₁

$\dfrac{x}{pq}<x^{1/3}<p$

$\phi\left(\dfrac{x}{pq},\pi\left(p\right)-1\right)=1$

$S_1=\dfrac{\left(\pi\left(y\right)-\pi\left(x^{1/3}\right)\right)\left(\pi\left(y\right)-\pi\left(x^{1/3}\right)-1\right)}{2}$

### 计算 S₂

$S_2=\sum\limits_{x^{1/4}<p\le x^{1/3}}{\ \sum\limits_{p<q\le y}{\phi\left(\dfrac{x}{pq},\pi\left(p\right)-1\right)}}$

$S_2=U+V$

$U=\sum\limits_{x^{1/4}<p\le x^{1/3}}{\ \sum\limits_{p<q<y\\q>x/p^2}{\phi\left(\dfrac{x}{pq},\pi\left(p\right)-1 \right)}}$

$V=\sum\limits_{x^{1/4}<p\le x^{1/3}}{\ \sum\limits_{p<q<y\\q\le x/p^2}{\phi\left(\dfrac{x}{pq},\pi\left(p\right)-1 \right)}}$

### 计算 U

$q>\dfrac x{p^2}$ 可得 $p^2>\dfrac xq\le \dfrac xy,p>\sqrt{\dfrac xy}$，因此：

$U=\sum\limits_{\sqrt{x/y}<p\le x^{1/3}}{\ \sum\limits_{p<q\le y\\q>x/p^2}{\phi\left(\dfrac{x}{pq},\pi\left(p\right)-1 \right)}}$

$U=\sum\limits_{\sqrt{x/y}<p\le x^{1/3}}{\#\left\{q\mid \dfrac x{p^2}<q\le y \right\}}$

$U=\sum\limits_{\sqrt{x/y}<p\le x^{1/3}}{\left(\pi\left(y\right)-\pi\left(\dfrac{x}{p^2} \right) \right)}$

### 计算 V

$\phi\left(\dfrac{x}{pq},\pi\left(p\right)-1 \right)=1+\pi\left(\dfrac{x}{pq} \right)-\left(\pi\left(p\right)-1\right)=2-\pi\left(p\right)+\pi\left(\dfrac{x}{pq} \right)$

$V=V_1+V_2$

$V_1=\sum\limits_{x^{1/4}<p\le x^{1/3}}{\ \sum\limits_{p<q\le \min\left(x/p^2,y\right)}{\left(2-\pi\left(p\right)\right)}}$

$V_2=\sum\limits_{x^{1/4}<p\le x^{1/3}}{\ \sum\limits_{p<q\le \min\left(x/p^2,y\right)}{\pi\left(\dfrac{x}{pq} \right)}}$

$V_2=\sum\limits_{x^{1/4}<p\le \sqrt{x/y}}{\ \sum\limits_{p<q\le y}{\pi\left(\dfrac{x}{pq} \right)}}+\sum\limits_{\sqrt{x/y}<p\le x^{1/3}}{\ \sum\limits_{p<q\le x/p^2}{\pi\left(\dfrac{x}{pq} \right)}}$

$V_2=W_1+W_2+W_3+W_4+W_5$

$W_1=\sum\limits_{x^{1/4}<p\le x/y^2}{\ \sum\limits_{p<q\le y}{\pi\left(\dfrac{x}{pq} \right)}}$

$W_2=\sum\limits_{x/y^2<p\le \sqrt{x/y}}{\ \sum\limits_{p<q\le \sqrt{x/p}}{\pi\left(\dfrac{x}{pq} \right)}}$

$W_3=\sum\limits_{x/y^2<p\le \sqrt{x/y}}{\ \sum\limits_{\sqrt{x/p}<q\le y}{\pi\left(\dfrac{x}{pq} \right)}}$

$W_4=\sum\limits_{\sqrt{x/y}<p\le x^{1/3}}{\ \sum\limits_{p<q\le \sqrt{x/p}}{\pi\left(\dfrac{x}{pq} \right)}}$

$W_5=\sum\limits_{\sqrt{x/y}<p\le x^{1/3}}{\ \sum\limits_{\sqrt{x/p}<q\le x/p^2}{\pi\left(\dfrac{x}{pq} \right)}}$

## 算法的时空复杂度

1. 计算 $P_2\left(x,a\right)$
2. 计算 $W_1,W_2,W_3,W_4,W_5$
3. 计算 $S_3$

### 计算 W₁,W₂,W₃,W₄,W₅ 的复杂度

$\pi\left(\dfrac{x}{y^2} \right)\pi\left(y\right)=O\left(\dfrac{x}{y\log^2 x} \right)$

$O\left(\sum\limits_{x/y^2<p\le \sqrt{x/y}}{\pi\left(\sqrt{\dfrac xp}\right)} \right)=O\left(\dfrac{x^{3/4}}{y^{1/4}\log^2 x} \right)$

$O\left(\sum\limits_{x/y^2<p\le \sqrt{x/y}}{\pi\left(\sqrt{\dfrac xp}\right)} \right)=O\left(\dfrac{x^{3/4}}{y^{1/4}\log^2 x} \right)$

$O\left(\sum\limits_{\sqrt{x/y}<p\le x^{1/3}}{\pi\left(\sqrt{\dfrac xp}\right)} \right)=O\left(\dfrac{x^{2/3}}{\log^2 x} \right)$

$O\left(\sum\limits_{\sqrt{x/y}<p\le x^{1/3}}{\pi\left(\sqrt{\dfrac xp}\right)} \right)=O\left(\dfrac{x^{2/3}}{\log^2 x} \right)$

### 计算 S₃ 的复杂度

$O\left(\dfrac{x}{y}\log x\log\log x+yx^{1/4}\right)$

### 总复杂度

$O\left(\dfrac{x}{y}\log{\log x}+\dfrac{x}{y}\log x\log{\log x}+x^{1/4}y+\dfrac{x^{2/3}}{\log^2{x}} \right)$

## 一些改进

• 终止条件 $2$ 中，我们可以用一个 $z$ 来代替 $y$，其中 $z$ 满足 $z>y$。我们可以证明这样子计算 $S_3$ 的时间复杂度可以优化到：

$O\left(\dfrac{x}{z}\log x\log{\log x}+\dfrac{yx^{1/4}}{\log x}+z^{3/2} \right)$

这也为通过改变 $z$ 的值来检查计算提供了一个很好的方法。

• 为了清楚起见，我们在阐述算法的时候选择在 $x^{1/4}$ 处拆分来计算总和 $S$，但实际上我们只需要有 $p\le \dfrac{x}{pq}<p^2$ 就可以计算。我们可以利用这一点，渐近复杂性保持不变。

• 用前几个素数 $2,3,5$ 预处理计算可以节省更多的时间。

## References

posted @ 2021-04-20 21:09  RioTian  阅读(80)  评论(1编辑  收藏