Meissel–Lehmer algorithm

The Meissel–Lehmer algorithm (after Ernst Meissel and Derrick Henry Lehmer) is an algorithm that computes exact values of 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 ≥ 3, the identities for $P_{k}(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.

History

Meissel already found that for k ≥ 3, $P_{k}(x,a)=0$ if $a=\pi (x^{1/3})$. He used the resulting equation for calculations of $\pi (x)$ for big values of $x$.

Meissel calculated $\pi (x)$ for values of x up to $10^{9}$, but he narrowly missed the correct result for the biggest value of x.

Using this method and an IBM 701, Lehmer was able to compute the correct value of $\pi \left(10^{9}\right)$ and missed the correct value of $\pi \left(10^{10}\right)$ by 1.

Modern variants

Jeffrey Lagarias, Victor Miller and Andrew Odlyzko published a realisation of the 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.

This extended Meissel-Lehmer algorithm needs less computing time than the algorithm developed by Meissel and Lehmer, especially for big values of x.

Further improvements of the algorithm are given by M. Deleglise and J. Rivat in 1996.