May 27 2012

## The Sieve of Eratosthenes

Published by

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:

1. Only check odd numbers for ‘primeness’. Hence we increment our loop counter by 2.
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:

1. The prime is already a known prime
2. 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.

• ## Time Marching On

 Hello world!