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

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>

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

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

unsigned long long nextblock ;
unsigned long long primecount ;

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

for (;;) {
lo = nextblock ;
hi = min(lo+sqlimit, limit) ;
nextblock = hi  ;

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 */
/* printf("%lld %lld %lld %lld\n", lo, hi, primecount, cnt) ; */
primecount += cnt ;
}

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 ;