Counting Primes

January 7, 2011

We celebrate our 200th exercise by implementing the prime-counting function, which has long been of interest to mathematicians. Generally denoted by the letter π, the prime-counting function π(n) returns the number of primes less than or equal to n; for instance, π(100) = 25. A related function, generally denoted as Pn, returns the nth prime number; thus, π(Pn) = n.

There are various analytic methods to compute the prime-counting function; the method invented by Ernst Meissel, improved by Derrick Henry Lehmer, and automated by Andrew Odlyzko is the best-known of them. But the obvious brute-force method of generating and counting the primes is hard to beat if you pre-compute a list of starting points, then use a segmented sieve to interpolate.

Your task is to write the prime-counting function and the nth-prime function. When you are finished, you are welcome to read or run a suggested solution, or to post your own solution or discuss the exercise in the comments below.

I don’t have time to do this properly today, so I’ll cheat with Haskell’s primes
package. It’s wicked fast, using “an efficient lazy wheel sieve for prime
generation.” I only tested up to 10,000, but the same code can work for 1e12
with patience and compilation. My compiler’s misbehaving today, though :-(

I wrote this one in a couple of forms a few months ago as a study in sieving. It’s not really pretty (nothing using pthreads is) but it chugs along pretty well. I can count all the primes up to 10^10 in about 12 seconds of elapsed time on my 8 processor machine, and it takes somewhere around 20 minutes to do all the way up to 10^12. There are actually ways to count primes that don’t involve finding all the primes, but they seemed impractical (hard to understand, and therefore hard to write).

I finally did this exercise to help me solve the Project Euler problems where having pi(x) comes in handy. The code is fairly big (156 lines of FORTH) so I won’t bore everyone with it, but here is a description of how I solved it:

First, for the table, a blocksize of 1e8 really makes little sense, since it takes 7.7 seconds to sieve a bucket that size (3.3 seconds in gforth-fast mode) and it takes less than one second to sieve a bucket of size 1e7 (in debug mode) However, the tables of pi for k*1e7 don’t go all the way to 1e12, so I wrote the code to generate the table myself and left my computer on overnight to generate it.

For the pi function, I use the 33K binary that contains the compressed primes < 1e6 from another exercise on this site. Writing pi() and nthprime() functions from the compressed table was probably the trickiest part of the code to get right. Lots of bit twiddling mistakes and off-by-one errors here.

The pi() function then considers three cases:
* For numbers 1e7, it uses the interpolation table for the base, then runs and counts a segmented sieve.

The next-prime function also has three cases:
* For numbers = pi(1e7), it does a binary search on the interpolation table, subtracts the pi value and runs a segmented sieve.

The segmented sieve code does not re-sieve if it is passed the same start and block size as what was previously computed, so there’s some caching for the nth-prime function.

Finally constantly shifting between 0 based indexing (the interpolation table and bit vectors) and 1 based indexing (the 0th entry in the pi table is really for 1e7, (i.e. k=1) and primes start counting at 1), was another part of the code that took a while to get right.