10

Can anyone provide computationally feasible pseudocode of any prime-counting function implementation? I initially attempted coding the Hardy-Wright algorithm, but its factorials began generating miserable overflow, 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.

Community
  • 1
  • 1
sevko
  • 1,321
  • 1
  • 17
  • 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

19

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
  • 16,961
  • 2
  • 32
  • 56
  • From the first code snippet it is implied that phi(x,a)=t, but here you stored phi(x, a-1)=t (If the latter was true, then there would be a much faster algorithm for this). I would have edited this myself, but we have the stupid restriction of the edit must be more than 6 characters. – Strategy Thinker Mar 09 '15 at 18:29
  • @StrategyThinker: Thanks for the correction. Fixed. – user448810 Mar 10 '15 at 03:51
  • 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
  • The line: w := n / p[i], seems to reference variable "n" but I don't see this variable anywhere. Did I miss something? – flcoder Jun 15 '19 at 10:54
  • It's x, not n. Fixed. Thank you. – user448810 Jun 15 '19 at 15:43
3

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
  • 31
  • 2