## Counting Primes

### January 7, 2011

Our aspirations are modest: we provide versions of `prime-pi`

and `nth-prime`

that work only to 10^{12}. We note that π(1e12) = 37607912018 and P_{37607912018} = 999999999989. We pre-compute a list of π(*n*) for each *n* from 10^{8} to 10^{12} in steps of 10^{8}, which requires us to store ten thousand numbers in a vector:

`(define (make-pi-vector)`

(with-output-to-file "pi-vector.ss" (lambda ()

(display "#(") (newline) (display 0) (newline)

(let* ((ps (cdr (primes (+ (isqrt #e1e12) 1))))

(qs (map (lambda (p) (modulo (* -1/2 (+ #e1e8 1 p)) p)) ps))

(z (length (primes #e1e8))))

(do ((t #e1e8 (+ t #e1e5 #e1e5))

(qs qs (map (lambda (p q) (modulo (- q #e1e5) p)) ps qs)))

((> t #e1e12) (display ")") (newline))

(when (zero? (modulo t #e1e8)) (display z) (newline))

(let ((bv (make-vector #e1e5 #t)))

(do ((qs qs (cdr qs)) (ps ps (cdr ps))) ((null? qs))

(do ((j (car qs) (+ j (car ps)))) ((<= #e1e5 j))

(vector-set! bv j #f)))

(do ((j 0 (+ j 1))) ((= j #e1e5))

(if (vector-ref bv j) (set! z (+ z 1))))))))))

Our program is adapted from the segmented sieve of Eratosthenes that was the subject of a previous exercise, except that we count the primes instead of enrolling them in a list. That program takes a while to run; on my desktop computer at home it took a day and a half:

`(time (make-pi-vector))`

1661621 collections

124715985 ms elapsed cpu time, including 9523237 ms collecting

128828796 ms elapsed real time, including 9843783 ms collecting

7273859064840 bytes allocated, including 7273849368416 bytes reclaimed

If you are impatient, you can look at the next page, where the vector has been folded into 80-column lines, or steal the table from Tomás Oliveira e Silva.

Once we have `pi-vector`

, it is a simple matter to count primes within the appropriate range using a segmented sieve of Eratosthenes; for *n*<10^{8}, we use the regular sieve, not the segmented sieve, because the segmented sieve does the wrong thing when its left end is less than the square root of its right end:

`(define (prime-pi n)`

(if (< #e1e12 n) (error 'prime-pi "out of range")

(if (< n #e1e8) (length (primes n))

(let* ((n (if (even? n) (- n 1) n))

(b (quotient n #e1e8)) (l (* b #e1e8))

(ps (cdr (primes (+ (isqrt n) 1))))

(qs (map (lambda (p) (modulo (* -1/2 (+ l 1 p)) p)) ps))

(z (vector-ref pi-vector b)) (pi #f))

(do ((t l (+ t #e2e5))

(qs qs (map (lambda (p q) (modulo (- q #e1e5) p)) ps qs)))

(pi pi)

(let ((bv (make-vector #e1e5 #t)))

(do ((qs qs (cdr qs)) (ps ps (cdr ps))) ((null? qs))

(do ((j (car qs) (+ j (car ps)))) ((<= #e1e5 j))

(vector-set! bv j #f)))

(do ((j 0 (+ j 1))) ((= j #e1e5))

(when (vector-ref bv j) (set! z (+ z 1)))

(when (= (+ t j j 1) n) (set! pi z)))))))))

Computing the *n*th prime is yet another variant of the segmented sieve. We find the right bucket in `pi-vector`

by binary search, then count until we find the desired prime:

`(define (nth-prime n)`

(define (bsearch n)

(let loop ((lo 0) (hi 10000))

(if (< hi lo) hi

(let ((mid (quotient (+ lo hi) 2)))

(cond ((< (vector-ref pi-vector mid) n) (loop (+ mid 1) hi))

((< n (vector-ref pi-vector mid)) (loop lo (- mid 1)))

(else (- mid 1)))))))

(if (< 37607912018 n) (error 'nth-prime "out of range")

(if (< n #e1e8) (list-ref (primes #e1e8) (- n 1))

(let* ((b (bsearch n)) (l (* b #e1e8))

(ps (cdr (primes (+ (isqrt (+ l #e1e8)) 1))))

(qs (map (lambda (p) (modulo (* -1/2 (+ l 1 p)) p)) ps))

(z (vector-ref pi-vector b)) (nth #f))

(do ((t l (+ t #e2e5))

(qs qs (map (lambda (p q) (modulo (- q #e1e5) p)) ps qs)))

(nth nth)

(let ((bv (make-vector #e1e5 #t)))

(do ((qs qs (cdr qs)) (ps ps (cdr ps))) ((null? qs))

(do ((j (car qs) (+ j (car ps)))) ((<= #e1e5 j))

(vector-set! bv j #f)))

(do ((j 0 (+ j 1))) ((= j #e1e5))

(when (vector-ref bv j) (set! z (+ z 1)))

(when (and (not nth) (= z n)) (set! nth (+ t j j 1))))))))))

We test the program by selecting a thousand random integers 1 ≤ *n* ≤ 37607912018 and checking that `(prime-pi (nth-prime n))`

is equal to *n*, and also by checking the ten smallest *n*, the ten largest *n*, and ten *n* divisible by #e1e8. We also compared results to The Nth Prime Page of Andrew Booker. We used `isqrt`

from the Standard Prelude and `primes`

from a previous exercise. You can run the program at http://programmingpraxis.codepad.org/XmEwn8eE.

This is too easy in J, which has a _primes_ operator: p:

First, define ‘n’ as the Number of Primes function: _1 p: N

n =: _1&p:

n 100

25

d,:n d=:100+i.10

100 101 102 103 104 105 106 107 108 109

25 25 26 26 27 27 27 27 28 28

p: num is the primes function, which by default returns the first primes.

The first 25 primes:

p: i.25

2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97

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).

[…] week Programming Praxis proposed to write two […]

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.