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 [...]

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

Follow

Get every new post delivered to your Inbox.

Join 611 other followers

%d bloggers like this: