Meissel–Lehmer algorithm

The Meissel–Lehmer algorithm (after Ernst Meissel and Derrick Henry Lehmer) is an algorithm that computes the prime-counting function.

Description

Key functions

Let $p_{1},p_{2},\ldots ,p_{n}$ be the first $n$ primes. For natural number a ≥ 1, define

$\varphi (x,a):=\left|\left\{n\leq x:p|n\implies p>p_{a}\right\}\right|,$

which counts natural numbers no greater than x with all prime factors greater $p_{a}$. Also define for natural number k,

$P_{k}(x,a):=\left|\left\{n\leq x:n=q_{1}q_{2}\cdots q_{k},~{\text{with}}~q_{1},\ldots ,q_{k}>p_{a}\right\}\right|,$

which counts natural numbers no greater than x with exactly k prime factors, all larger than $p_{a}$. With these, we have

$\varphi (x,a)=\sum _{k=0}^{\infty }P_{k}(x,a),$

where the sum only has finitely many nonzero terms, because $P_{k}(x,a)=0$ when $p_{a}^{k}>x$. Using the fact that $P_{1}(x,a)=\pi (x)-a$, we get

$\pi (x)=\varphi (x,a)+a-1-\sum _{k=2}^{\infty }P_{k}(x,a),$

which prove that one may compute $\pi (x)$ by computing $\varphi (x,a)$ and $P_{k}(x,a)$ for k ≥ 2. This is what the Meissel–Lehmer algorithm does.

Formula for P_k(x, a)

For $k=2$, we get the following formula for $P_{k}(x,a)$:

${\begin{aligned}P_{2}(x,a)&=\left|\left\{n:n\leq x,~n=p_{b}p_{c},~{\text{with}}~a<b\leq c\right\}\right|\\&=\sum _{b=a+1}^{\pi (x^{1/2})}\left|\left\{n:n\leq x,~n=p_{b}p_{c},~{\text{with}}~b\leq c\leq \pi \left({\frac {x}{p_{b}}}\right)\right\}\right|\\&=\sum _{b=a+1}^{\pi (x^{1/2})}\left(\pi \left({\frac {x}{p_{b}}}\right)-(b-1)\right)\\&={\binom {a}{2}}-{\binom {\pi (x^{1/2})}{2}}+\sum _{b=a+1}^{\pi (x^{1/2})}\pi \left({\frac {x}{p_{b}}}\right).\end{aligned}}$

For $k\geq 3$, the identities for $P_{3}(x,a)$ can be derived similarly.

Expanding 𝜑(x, a)

Using the recurrence

$\varphi (x,a)=\varphi (x,a-1)-\varphi \left({\frac {x}{p_{a}}},a-1\right),$

$\varphi (x,a)$ may be expanded. Each summand, in turn, may be expanded in the same way.

Combining the terms

The only thing that remains to be done is evaluating $\varphi (x,a)$ and $P_{k}(x,a)$ for k ≥ 2, for certain values of $x$ and $a$. This can be done by direct sieving and using the above formulas.

A modern variant

Jeffrey Lagarias, Victor Miller and Andrew Odlyzko published a realisation of this algorithm which computes $\pi (x)$ in time $O(x^{2/3+\varepsilon })$ and space $O(x^{1/3+\varepsilon })$ for any $\varepsilon >0$. Upon setting $a=\pi (x^{1/3})$, the tree of $\varphi (x,a)$ has $O(x^{2/3})$ leaf nodes.