Programming Praxis

4 Responses to “Counting Primes”

1. Alan said

   January 7, 2011 at 10:07 AM

   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

   January 7, 2011 at 4:24 PM

   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. Mark VandeWettering said

   January 8, 2011 at 1:39 AM

   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. Counting primes in Factor | Alfredo Beaumont said

   January 12, 2011 at 11:08 AM

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