Meissel–Lehmer algorithm

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

Description

The key identity

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

${\displaystyle \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 ${\displaystyle [1,x]\cap \mathbb {N} }$ (closed interval) when one has sieved out multiples of the first ${\displaystyle a}$ primes ${\displaystyle p_{1},\ldots ,p_{a}}$ (including these primes themselves); that is, ${\displaystyle 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

${\displaystyle 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|,}$

${\displaystyle k\in \mathbb {N} _{0};}$

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

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

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

The identities

${\displaystyle \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 ${\displaystyle \pi (x)}$ by computing ${\displaystyle \varphi (x,a)}$ and ${\displaystyle P_{k}(x,a)}$, ${\displaystyle k\geq 2}$. And this is what the Meissel–Lehmer algorithm does.

Formulae for the P_k(x, a)

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

${\displaystyle {\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 ${\displaystyle k=3}$, there is a similar expansion.[2]

Expanding 𝜑(x, a)

Using the formula

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

${\displaystyle \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 ${\displaystyle \varphi (x,a)}$ and ${\displaystyle P_{k}(x,a),k=1,2,\ldots }$ for certain values of ${\displaystyle x}$ and ${\displaystyle 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 ${\displaystyle \pi (x)}$ in time ${\displaystyle O(x^{2/3+\varepsilon })}$ and space ${\displaystyle O(x^{1/3+\varepsilon })}$ for any ${\displaystyle \varepsilon >0}$.[1] Upon setting ${\displaystyle a=\pi (x^{1/3})}$, the tree of ${\displaystyle \varphi (x,a)}$ has ${\displaystyle O(x^{2/3})}$ leaf nodes.[1]

References

1. Lagarias, Jeffrey; Miller, Victor; Odlyzko, Andrew (April 11, 1985). "Computing ${\displaystyle \pi (x)}$: The Meissel–Lehmer method" (PDF). Mathematics of Computation. 44 (170): 537–560. doi:10.1090/S0025-5718-1985-0777285-5. Retrieved September 13, 2016.
2. Lehmer, Derrick Henry (April 1, 1958). "ON THE EXACT NUMBER OF PRIMES LESS THAN A GIVEN LIMIT". Illinois J. Math. 3 (3): 381–388. Retrieved February 1, 2017.