## Counting Primes

### January 7, 2011

Our aspirations are modest: we provide versions of `prime-pi` and `nth-prime` that work only to 1012. We note that π(1e12) = 37607912018 and P37607912018 = 999999999989. We pre-compute a list of π(n) for each n from 108 to 1012 in steps of 108, 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<108, 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 nth 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.

About these ads

Pages: 1 2 3

### 4 Responses to “Counting Primes”

1. Alan said

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 ```

2. Graham said

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 :-(

```import Data.Numbers.Primes

primesLessThan :: (Integral a) => a -> [a]
primesLessThan n = takeWhile (< n) primes

primePi :: (Integral a) => a -> Int
primePi = length . primesLessThan

nthPrime :: (Integral a) => Int -> a
nthPrime n = primes !! n

test :: Int -> Bool
test n = (primePi (nthPrime n)) == n

main :: IO ()
main = do
print \$ primesLessThan 100
print \$ all test [1..(round 1e4)]
```
3. 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).

```#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <pthread.h>

unsigned long long limit = 1000000000000LL ;
unsigned long long sqlimit, sqlimit2 ;
unsigned char *b0 ;

#define min(a, b)       ((a)<(b)?(a):(b))

int nthreads = 8 ;

pthread_mutex_t mutex ;
unsigned long long nextblock ;
unsigned long long primecount ;

void *
workerThread(void *parm)
{
unsigned long long lo, hi, m ;
unsigned long long idx, idx2, cnt = 0LL ;
unsigned char *b1 = (unsigned char *) parm ;
int i, rc ;

for (;;) {
pthread_mutex_lock(&mutex) ;
lo = nextblock ;
hi = min(lo+sqlimit, limit) ;
nextblock = hi  ;
pthread_mutex_unlock(&mutex) ;

if (lo >= limit)
break ;

/* now, sieve */
memset(b1, 0, sizeof(*b1)*sqlimit) ;
for (idx=2; idx<sqlimit; idx++) {
if (b0[idx] == 0) {
/* it's a prime, strike out multiples of it in b1 */
m = lo / idx ;
idx2 = m * idx ;
if (idx2 < lo)
idx2 += idx ;
for ( ; idx2 < hi; idx2 += idx)
b1[idx2-lo] = 1 ;
}
}
for (idx=lo, cnt=0LL; idx < hi; idx ++) {
if (b1[idx-lo] == 0) {
/* printf("%lld\n", idx) ; */
cnt ++ ;
}
}

/* add in the total */
pthread_mutex_lock(&mutex) ;
/* printf("%lld %lld %lld %lld\n", lo, hi, primecount, cnt) ; */
primecount += cnt ;
pthread_mutex_unlock(&mutex) ;
}

return NULL ;
}

main(int argc, char *argv[])
{
unsigned long long idx, idx2, cnt = 0LL ;
unsigned long long lo, hi, m ;
int i, rc ;

if (argc > 1)
limit = atoll(argv[1]) ;
sqlimit = (unsigned long long) ceil(sqrt((double) limit)) ;
sqlimit2 = (unsigned long long) ceil(sqrt((double) sqlimit)) ;

/* okay, allocate two buffers */
b0 = (unsigned char *) calloc(sqlimit, sizeof(unsigned char)) ;

/* sieve to find the first sqlimit primes */

for (idx=2; idx<=sqlimit2; idx++)
if (b0[idx] == 0) {
for (idx2 = idx*idx; idx2 < sqlimit; idx2 += idx)
b0[idx2] = 1 ;
}

for (idx=2LL; idx < sqlimit; idx++)
if (b0[idx] == 0) {
/* printf("%lld\n", idx) ; */
cnt ++ ;
}

fprintf(stderr, "%lld primes less than %lld\n", cnt, sqlimit) ;

/* now, fork a bunch of threads... */

nextblock = sqlimit ;
primecount = cnt ;
pthread_mutex_init(&mutex, NULL) ;

pthread_t *thread = (pthread_t *) calloc(nthreads, sizeof(pthread_t)) ;

for (i=0; i<nthreads; i++) {
rc = pthread_create(&thread[i], NULL, workerThread, (void *) calloc(sqlimit, sizeof(unsigned char))) ;
}

/* wait for everyone to finish */

for (i=0; i<nthreads; i++) {
pthread_join(thread[i], NULL) ;
}

printf("%lld primes less than %lld\n", primecount, limit) ;
}
```
4. [...] weekÂ Programming Praxis proposed to write two [...]