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\mid 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=a+\varphi (x,a)-1-\sum _{k=2}^{\infty }P_{k}(x,a)\qquad {\text{(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 Pk(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

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 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^{2/3})$ , the tree of $\varphi (x,a)$  has $O(x^{2/3})$  leaf nodes.