Possible implementation of Prime Counting

Question

Can someone provide a computationally possible pseudo-code of any number counting function ? At first I tried to code the Hardy-Wright algorithm, but its factorials started to generate a miserable overflow, and many others seem to inevitably lead to similar problems. I searched Google for practical solutions, but at best I found very esoteric maths that I had never seen in regular programs.

Answer 1

The prime count function pi (x) calculates 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 the helper function phi (x, a), which counts numbers no more than x that are not affected by the sifting by the first prime numbers; for example, phi (50.3) = 14 for numbers 1, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, and 49. The phi function can be calculated 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 a prime (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 

Array p stores the a-th prime number for small a, calculated by sifting. Cache is important; without this, the lead time would be exponential. For this fi, the formula for calculating Legendre primes is of the form 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 78498, which, although it was incorrect, was surprisingly close to complicated manual calculation.

In the 1950s, Derrick Lehmer gave an improved algorithm for calculating prime numbers:

 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, their modern versions, and very fast computers, it remains extremely tedious to calculate large values ​​of pi; At the time of writing, the largest known value is pi (10 ^ 24) = 18435599767349200867866.

To use this algorithm to calculate the nth prime number, a consequence of the prime theorem restricts the nth prime number P (n) between n log n and n (log n + log log n) for n> 5; therefore, we calculate pi at the borders and use halving to determine the nth prime, switching to sieving when the borders are close.

I discuss prime numbers on several posts on my blog .

Answer 2

Wikipedia can also help. The article on counting primes contains several pointers. For starters, I would recommend the Meissel algorithm in the section "π (x) Estimation Algorithms", which is one of the simplest algorithms that does not generate all primes.

I also find Pomerance and Crandall's book , Prime Numbers in a Computational Perspective, useful. This book has a detailed and fairly accessible description of prime counting methods. But keep in mind that the topic is inherently too complex for most readers here.