Computing $\pi (x)$: the Meissel-Lehmer method

Abstract:

E. D. F. Meissel, a German astronomer, found in the 1870’s a method for computing individual values of $\pi (x)$, the counting function for the number of primes $\leqslant x$. His method was based on recurrences for partial sieving functions, and he used it to compute $\pi ({10^9})$. D. H. Lehmer simplified and extended Meissel’s method. We present further refinements of the Meissel-Lehmer method which incorporate some new sieving techniques. We give an asymptotic running time analysis of the resulting algorithm, showing that for every $\varepsilon > 0$ it computes $\pi (x)$ using at most $O({x^{2/3 + \varepsilon }})$ arithmetic operations and using at most $O({x^{1/3 + \varepsilon }})$ storage locations on a Random Access Machine (RAM) using words of length $[{\log _2}x] + 1$ bits. The algorithm can be further speeded up using parallel processors. We show that there is an algorithm which, when given M RAM parallel processors, computes $\pi (x)$ in time at most $O({M^{ - 1}}{x^{2/3 + \varepsilon }})$ using at most $O({x^{1/3 + \varepsilon }})$ storage locations on each parallel processor, provided $M \leqslant {x^{1/3}}$. A variant of the algorithm was implemented and used to compute $\pi (4 \times {10^{16}})$.