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.

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 635 other followers

%d bloggers like this: