Meissel–Lehmer algorithm

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

Description

The key identity

Given $a\in \mathbb {N}$, one may define the following functions: Firstly,

$\varphi (x,a):=\left|\left([1,x]\cap \mathbb {N} \right)\setminus \bigcup _{j=1}^{a}p_{j}\mathbb {Z} \right|$;

this function measures the set $[1,x]\cap \mathbb {N}$ (closed interval) when one has sieved out multiples of the first $a$ primes $p_{1},\ldots ,p_{a}$ (including these primes themselves); that is, $p_{1},p_{2},p_{3},\ldots =2,3,5,\ldots$ is the sequence of prime numbers in increasing order.

We may further split that up with the functions

$P_{k}(x,a):=\left|\left[\left([1,x]\cap \mathbb {N} \right)\setminus \bigcup _{j=1}^{a}p_{j}\mathbb {Z} \right]\cap \{n|n{\text{ has exactly }}k{\text{ prime factors}}\}\right|$, $k\in \mathbb {N} _{0}$;

these measure the set $\left([1,x]\cap \mathbb {N} \right)\setminus \bigcup _{j=1}^{a}p_{j}\mathbb {Z}$ but considers only the numbers that have exactly $k$ prime factors (this is well-defined by the fundamental theorem of arithmetic). With these, we have

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

the sum on the right is finite, since numbers $\leq x$ have, for instance, $\leq \lfloor x\rfloor$ prime factors.

The identities

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

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

Formulae for the 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\mid n\leq x,n=p_{b}p_{c},a<b\leq c\right\}\right|\\&=\sum _{b=a+1}^{\pi (x^{1/2})}\left|\left\{n~{\big |}~n=p_{b}p_{c},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)\\&=\sum _{b=a+1}^{\pi (x^{1/2})}\pi \left({\frac {x}{p_{b}}}\right)+{\binom {a}{2}}-{\binom {\pi (x^{1/2})}{2}}.\end{aligned}}$

For $k=3$, there is a similar expansion.

Expanding ϕ(x, a)

Using the formula

$\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, so that a very large alternating sum arises.

Combining the terms

The only thing that remains to be done is evaluating $\varphi (x,a)$ and $P_{k}(x,a),k=1,2,\ldots$ 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 crafted 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^{2/3})$, the tree of $\varphi (x,a)$ has $O(x^{2/3})$ leaf nodes.