10

What would be a computationally feasible pseudocode of any prime-counting function implementation?

I initially attempted coding the Hardy-Wright algorithm, but its factorials began generating miserable overflows, and many others appear bound to yield similar problems. I've scoured Google for practical solutions, but, at best, have found very esoteric mathematics which I haven't ever seen implemented in conventional programs.

Peter Mortensen
  • 30,738
  • 21
  • 105
  • 131
sevko
  • 1,402
  • 1
  • 18
  • 37
  • 1
    ins't it true that `floor(x/j) * j == x - (x%j)`. Then the formula you link to becomes `pi(x) = (-1) + SUM{j=3..n}( ((j-2)!) % j )` (?). Next use modular multiplication (i.e. `5! % 7 == (((((2*3)%7)*4)%7)*5)%7`). – Will Ness Sep 29 '13 at 13:44
  • @SeverynKozak The reason they say this is not a programming issue question is because your question contains no code. – brian beuning Sep 29 '13 at 16:05
  • As to providing much more than just pseudocode, [my answer in another SO thread](https://stackoverflow.com/a/70049343/549617) provides a working code snippet (runnable in the browser) of the basic Legarias-Miller-Odzlinko (LMO) 1985 algorithm/method in JavaScript. LMO forms the basis of all modern practical general prime counting functions, and a variation (Gourdon, 2001) has been used to compute the number of primes up to 1e28, although that takes months even using multi-threading on a modern computer with a very large number of cores and amount of memory. – GordonBGood Nov 22 '21 at 03:16

2 Answers2

22

The prime-counting function pi(x) computes the number of primes not exceeding x, and has fascinated mathematicians for centuries. At the beginning of the eighteenth century, Adrien-Marie Legendre gave a formula using an auxiliary function phi(x,a) that counts the numbers not greater than x that are not stricken by sieving with the first a primes; for instance, phi(50,3) = 14 for the numbers 1, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47 and 49. The phi function can be computed as phi(x,a) = phi(x,a-1) - phi(x/P(a),a-1), where phi(x,1) is the number of odd integers not exceeding x and P(a) is the a-th prime number (counting from P(1)=2).

function phi(x, a)
  if (phi(x, a) is in cache)
    return phi(x, a) from cache
  if (a == 1)
    return (x + 1) // 2
  t := phi(x, a-1) - phi(x // p[a], a-1)
  insert phi(x, a) = t in cache
  return t

An array p stores the a-th prime for small a, calculated by sieving. The cache is important; without it, run time would be exponential. Given phi, Legendre's prime-counting formula is pi(x) = phi(x,a) + a - 1, where a = pi(floor(sqrt(x))). Legendre used his formula to calculate pi(10^6), but he reported 78526 instead of the correct answer of 78498, which, even though wrong, was astonishingly close for an intricate manual calculation.

In the 1950s, Derrick H. Lehmer gave an improved algorithm for counting primes:

function pi(x)
  if (x < limit) return count(primes(x))
  a := pi(root(x, 4)) # fourth root of x
  b := pi(root(x, 2)) # square root of x
  c := pi(root(x, 3)) # cube root of x
  sum := phi(x,a) + (b+a-2) * (b-a+1) / 2
  for i from a+1 to b
    w := x / p[i]
    lim := pi(sqrt(w))
    sum := sum - pi(w)
    if (i <= c)
      for j from i to lim
        sum := sum - pi(w / p[j]) + j - 1
  return sum

For example, pi(10^12) = 37607912018. Even with these algorithms, and their modern variants, and very fast computers, it remains appallingly tedious to calculate large values of pi; at this writing, the largest known value is pi(10^24) = 18435599767349200867866.

To use this algorithm to calculate the n-th prime, a corollary to the Prime Number Theorem bounds the n-th prime P(n) between n log n and n(log n + log log n) for n > 5, so compute pi at the bounds and use bisection to determine the n-th prime, switching to sieving when the bounds are close.

I discuss prime numbers in several entries at my blog.

user448810
  • 17,381
  • 4
  • 34
  • 59
  • 2
    Recently, pi(10^25) and pi(10^26) have been calculated. See [page 40 here](http://dalspace.library.dal.ca/handle/10222/60524). – qwr Feb 10 '16 at 12:41
  • 1
    Your post was very helpful in my own implementation. Memoizing pi(x) as well is a large speedup. I had some rounding issues with calculating a, b, c so I used danaj's rounding: https://programmingpraxis.com/2011/07/22/counting-primes-using-legendres-formula/#comment-5958 – qwr Jun 08 '17 at 21:37
2

Wikipedia might help too. The article on prime counting contains a few pointers. For starters I'd recommend the algorithm by Meissel in the section "Algorithms for evaluating π(x)", which is one of simplest algorithm that does not generate all primes.

I also find the book by Pomerance and Crandall "Prime numbers a computational perspective" helpful. This book has a detailed and quite accessible description of prime counting methods. But keep in mind that the topic by its nature is a bit too advanced for most of the reader here.

code monkey
  • 29
  • 2