May 27 2012

## The Sieve of Eratosthenes

This is an iterative solution to the Sieve of Eratosthenes problem.

This implementation cheats a little. First, primes are primes, so we can assume something about the sequence of prime numbers and take advantage of that knowledge. The first fact is that 2 is the only even prime number – all others are odd. Also, by the Prime Number Theorem we “assume” that all factors of a number *n* must be less than the square root of *n*. So, how do we exploit this knowledge to keep our processing down:

- Only check odd numbers for ‘primeness’. Hence we increment our loop counter by 2.
- We only test previously discovered primes < = sqrt(n) as factors for n

// return the nth prime number int getNthPrime(unsigned n) { assert(n>0); // if they ask for the first prime, cheat and just return it. if( n == 1 ) return 2; vector<unsigned > primes; // collection of primes (sans 2) // we start with 3, and go up to n-1. As 2 is the first prime, // and we've already accounted for it above, we only need to // find the nth-1 odd prime. for(unsigned cand(3); primes.size() < n; cand += 2) { unsigned upper_bound = static_cast<unsigned >( floor(sqrt(static_cast<double>( cand ))) ); bool isAPrime(true); for( auto itr = primes.begin(); itr != primes.end() && *itr < = upper_bound; itr++ ) if( (cand % *itr) == 0 ) { isAPrime = false; break; } if( isAPrime ) primes.push_back( cand ); } return primes.back() // return the last prime discovered }

This works fairly well, except that every time we request prime *n* it starts at 3 and works up to the *nth* prime. An improvement would be to make the vector static so that known primes are cached. Then for any requested prime we would check if:

- The prime is already a known prime
- If not, then start looking just after the last known prime

Here’s another version with these improvements:

unsigned getNthPrime(unsigned n) { // the list of known primes. The indexes for a given prime n will be // biased by 2. As we cheat with the first prime - 2 - then the index // for the second prime - 3 - will be 0, the third prime 1, and so on. static std::vector<unsigned> primes; assert(n>0); if( n == 1 ) return 2; // see if we already know about this prime. if( --n < primes.size() ) return primes[ n-1 ]; // all primes > 2 are odd. Since we know that 2 is a prime, // start at 3 and work up looking only at odds. unsigned cand = (primes.empty())? 3 : primes.back() + 2; for(; primes.size() < n; cand += 2) { unsigned upper_bound = static_cast<unsigned>( floor(sqrt(static_cast<double>( cand ))) ); bool isAPrime(true); for(auto itr = primes.begin(); itr != primes.end() && *itr <= upper_bound; itr++ ) if( (cand % *itr) == 0 ) { isAPrime=false; break; } if( isAPrime ) primes.push_back( cand ); } return primes[n-1]; }

The next step will be to see – just for fun – how this could be done with TMP.