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.^[1][2]

Description

The problem of counting the exact number of primes less than or equal to x, without actually listing them all, dates from Legendre. He observed from the Sieve of Eratosthenes that

${\displaystyle \pi (x)-\pi (x^{1/2})+1=\lfloor x\rfloor -\sum _{i}\lfloor x/p_{i}\rfloor +\sum _{i<j}\lfloor x/p_{i}p_{j}\rfloor -\ldots }$

where ${\displaystyle \lfloor x\rfloor }$ is the floor function, which denotes the greatest integer less than or equal to x and the ${\displaystyle p_{i}}$ run over all primes ${\displaystyle \leq x^{1/2}}$.^[1][2]

Since the evaluation of this sum formula becomes more and more complex and confusing for large x, Meissel tried to simplify the counting of the numbers in the Sieve of Eratosthenes. He and Lehmer therefore introduced certain sieve functions, which are detailed below.

Key functions

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

${\displaystyle \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 than ${\displaystyle p_{a}}$. Also define for a natural number k,

${\displaystyle 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 greater than ${\displaystyle p_{a}}$. With these, we have

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

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

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

which proves that one may compute ${\displaystyle \pi (x)}$ by computing ${\displaystyle \varphi (x,a)}$ and ${\displaystyle 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 ${\displaystyle P_{k}(x,a)}$:

${\displaystyle {\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 ${\displaystyle P_{k}(x,a)}$ can be derived similarly.^[1]

Expanding 𝜑(x, a)

With the starting condition

${\displaystyle \varphi (x,0)=\lfloor x\rfloor ,}$

and the recurrence

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

each value for ${\displaystyle \varphi (x,a)}$ can be calculated recursively.

Combining the terms

The only thing that remains to be done is evaluating ${\displaystyle \varphi (x,a)}$ and ${\displaystyle 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, ${\displaystyle P_{k}(x,a)=0}$ if ${\displaystyle a=\pi (x^{1/3})}$. He used the resulting equation for calculations of ${\displaystyle \pi (x)}$ for big values of ${\displaystyle x}$. ^[1]

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

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

Extended algorithm

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

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.^[3][4]

References