COMPUTING p((x)): THE MEISSEL-LEHMER METHOD
J. C. Lagarias
AT&T Bell Laboratories
Murray Hill, New Jersey
V. S. Miller
IBM Watson Research Center
Yorktown Heights, New York
A. M. Odlyzko
AT&T Bell Laboratories
Murray Hill, New Jersey
ABSTRACT
E. D. F. Meissel, a German astronomer, found in the 1870’s a method for computing
individual values of p(x), the counting function for the number of primes £ x. His
method was based on recurrences for partial sieving functions, and he used it to compute
p( 109 ). 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 e > 0 it computes p(x) using at most O(x 2/3 + e ) arithmetic
operations and using at most O(x 1/3 + e ) 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 p(x) in time at most O(M- 1 x 2/3 + e ) using at most
O(x 1/3 + e ) storage locations on each parallel processor, provided M £ x 1/3. A variant of
the algorithm was implemented and used to compute p( 4´1016 ).
COMPUTING p((x)): THE MEISSEL-LEHMER METHOD
J. C. Lagarias
AT&T Bell Laboratories
Murray Hill, New Jersey
V. S. Miller
IBM Watson Research Center
Yorktown Heights, New York
A. M. Odlyzko
AT&T Bell Laboratories
Murray Hill, New Jersey
1. Introduction
The problem of computing p(x), the number of primes p £ x, has been studied for a
long time. The ancient Greeks had a method for locating all the prime numbers below a
given bound, the sieve of Eratosthenes. Legendre (see [4]) was the first to suggest a
method of calculating p(x) without locating all the primes less than x. He observed that
the principle of inclusion-exclusion implied that
p(x) - p(x 1/2 ) + 1 = [x] - S [
p i
__x _ ] + S [
p i p j
__ _x__ ] - S [
p i p j p k
__ __x___ ] + ... (1.1)
where [z] denotes the greatest integer £ z and the p i run over all primes £ x 1⁄2 , and
p i < p j in the second sum, p i < p j < p k in the third sum, and so on. This formula is
not of immediate utility in computing p(x) because it involves computing about
6p- 2 ( 1 - log 2 ) x nonzero terms. Actual calculations of p(x) by Gauss and others were
based on factor tables made using sieve methods. The first efficient algorithm for
computing p(x) which does not involve locating all the primes below x is due to the
astronomer E. D. F. Meissel ([11]-[14]). His method is economical of space, and can be
viewed as reducing the number of terms in ‘‘Legendre’s sum’’ (1.1). Using his methods
- 2 -
he calculated p( 108 ) = 5 , 761 , 455 in 1871, and then after an enormous calculation
announced in 1885 that p( 109 ) was 50,847,478. (His value of p( 109 ) was subsequently
shown to be too small by 56, see [9].) Modifications of Meissel’s method were
subsequently suggested by several authors (see [4], pp. 429-434) who did not, during the
era of hand computation, attempt to explicitly compute larger values of p(x). With the
advent of digital computers, the problem of computing p(x) was reconsidered by
D. H. Lehmer [9], who extended and simplified Meissel’s method. Lehmer used an IBM
701 to calculate p( 1010 ) = 455 , 052 , 512 (a value later shown [3] to be too large by 1).
Following this, other authors ([3],[10]) calculated p(x) for larger values of x using
variants of the Meissel-Lehmer method, the current record being the calculation of
p( 1013 ) by Bohman [3] in 1972. Bohman’s value for p( 1013 ) turns out to be too small
by 941, see Section 6.
In this paper we propose new algorithms of Meissel-Lehmer type for computing p(x)
and analyze their asymptotic computational complexity. Our interest in the asymptotic
complexity of computing p(x) was stimulated by a paper of H. S. Wilf [15], which cited
p(x) as an example of a function whose individual values are hard to compute. Indeed in
this connection it appears that the existing methods for computing p(x), which are based
on variants of the Meissel-Lehmer method, have asymptotic running times at least
c e x 1 - e for any e > 0 and some c e > 0. (It is hard to estimate the asymptotic
computational complexity of Meissel’s original method [11], since it is presented as a
collection of rules to be applied according to the human calculator’s judgment.) We
analyze asymptotic running times using as a model of computation a Random Access
Machine (RAM), which is a relatively realistic model of the addressible core storage area
- 3 -
of a digital computer. (The RAM model is described in detail in Aho, Hopcroft and
Ullman [1], Chapter 1.) An essential feature of the RAM model relevant to the
complexity analysis is that on a RAM it is possible to jump between any two storage
locations in a single operation. This feature is needed to quickly implement algorithms
which use sieve methods. We prove the following result.
Theorem A. The Extended Meissel-Lehmer algorithm computes p(x) on a Random
Access Machine using at most O(x 2/3 + e ) arithmetic operations and at most O(x 1/3 + e )
storage locations, for any fixed e > 0. All integers used in the course of the computation
have at most [ log 2 x] + 1 bits in their binary expansions.
By arithmetic operations we mean additions, subtractions, multiplications and
divisions. Since any such operation on k bit integers takes O(k 2 ) bit operations using the
usul algorithms, the difference between counting arithmetic operations and bit operations
is a multiplicative factor of O( ( log x)2 ). Consequently Theorem A implies that the
Extended Meissel-Lehmer method computes p(x) in O(x 2/3 + e ) bit operations using
O(x 1/3 + e ) bit locations of storage.
The constants implied by the O-symbols in Theorem A depend on e > 0. The O(x e )
terms in Theorem A can be replaced by lower order terms by a more detailed running
time analysis than we give; for example, it can be shown that the algorithm described in
Section 3 halts after at most O(x 2/3 ( log x)4 ) bit operations using at most
O(x 1/3 ( log x)2 ) bit locations of space. At the end of Section 3 we explain why we are
unable to improve the exponent of the running time bound below 2/3 for an algorithm of
Meissel-Lehmer type, even if more space locations are available. On the other hand, one
- 4 -
can find variants of the Extended Meissel-Lehmer algorithm which use less space,
provided a running time of more than O(x 2/3 ) is allowed.
In addition, one can further speed up the computation of p(x) using parallel
processing. We consider a model of parallel processing in which simultaneous transfers
of information are allowed. We assume that we have M processors, which are arranged
as the first M leaves of a complete balanced binary tree having 2k leaves, where
2k - 1 < M &pound; 2k, and that the nonleaf nodes are occupied by switches. In a single time
unit, a bit of information can be transmitted over any of the edges in the tree together
with its destination (which will consist of élog 2 Mù bits). This model allows many
messages to be transmitted simultaneously between different processors, provided the
message paths do not intersect. We propose an algorithm called the Extended Parallel
Meissel-Lehmer method, and prove in Section 4 the following result.
Theorem B. For any M with 1 &pound; M &pound; x 1/3 the Extended Parallel Meissel-Lehmer
algorithm computes p(x) using M Random Access Machine parallel processors using at
most O(M- 1 x 2/3 + e ) arithmetic operations, and at most O(x 1/3 + e ) storage locations for
each parallel processor, i.e. O(Mx1/3 + e ) storage locations in all, for any e > 0. All
integers used in the course of the computation have length at most [ log 2 x] + 1 bits in
their binary expansions.
The constants implied by the O-symbols in Theorem B depend on e > 0 and not on M in
the range 1 &pound; M &pound; x 1/3. The O(x e ) terms in Theorem B can be replaced by lower order
terms like O( ( log x)4 ) by a more detailed running time analysis of the algorithm
described in the text.
- 5 -
Theorem B logically includes Theorem A. We prove Theorem A separately because
a somewhat simplified algorithm and running time analysis is possible in this case.
The algorithms described in Sections 3 and 4 were designed to facilitate their
asymptotic running time analysis, and have inefficiences from the viewpoint of practical
computation. In Section 5 we describe modifications of the Extended Meissel-Lehmer
method which improve its performance in practice. The second author implemented a
version of the algorithm described in Section 5 on an IBM 3081, and used it to compute
various values of p(x), the largest being x = 4&acute;1016. These results are described in
Section 6.