Fastest Algorithm to generate primes in c/c++ – Sieve of Eratosthenes

Posted by Holmez on 26 Oct 2010. Filed under Optimization Techniques, Programming. You can follow any responses to this entry through the RSS 2.0. You can leave a response or trackback to this entry

In mathematics, the Sieve of Eratosthenes is a simple, ancient algorithm for finding all prime numbers up to a specified integer. It works efficiently for the smaller primes (below 10 million). It was created by Eratosthenes, an ancient Greek mathematician.
To see how efficiency this algorithm is we will take a challenging problem from projecteuler.net (problem 7) and see how fast we can implement it.

Problem:

By listing the first six prime numbers: 2, 3, 5, 7, 11, and 13, we can see that the 6th prime is 13.
What is the 10001st prime number?

Code:

#include <iostream>
#include "time.h"

using namespace std;

#define SET0(i) List[i>>5]&=~(1<<i)
#define SET1(i) List[i>>5]|=(1<<i)
#define GET(i)  ((List[i>>5]>>i)&1)

#define NUM 150000 //expected upper limit of the required prime num
#define SQRT 388 //ceeling of the squar root of the NUM

unsigned long List[NUM>>5] = {0};

void FillList(long num)
{
	long t;
	num = (num << 1) + 1;
	t = (num*num - 1) >> 1;
	while(t<(NUM>>1))
	{
		SET1(t);
		t += num;
	}
}
void MakeList()
{
	SET1(0);
	for(long i = 1; i<=((SQRT-1)>>1); i++)
	{
		if(!GET(i))
		{
			FillList(i);
		}
	}
}

int main()
{
	clock_t start, end;
	start = clock();

	int counter = 1; //number two counted

	MakeList();

	for(long i=1; i<(NUM>>1); i++)
	{
		if(!GET(i))
		{
			counter++;
			if(counter==10001)
			{
				cout<<"10001st prime: "<<(i<<1)+1<<endl;
				break;
			}
		}
	}
	end = clock();
	cout<<"Done! Total time: "<<end-start<<endl;
}

Output:

Analysis:

Using initially the Sieve of Eratosthenes method to generate primes, but rather than using array of bool for each number which takes 8 bits per element. I used an array of int (actually long to make sure it will be 4 bytes at all processors) and each int has a 32 bits, each one represent a 0 for prime and 1 for not prime.

Further more we will not looping through all numbers as we know there are no even primes (except 2) so our index of looping (i) will change following the equation

prime = i*2 +1

No prime will get out of this equation (except 2) and it will reduce the looping times to less than its half.

Now when we want to set number 13 to be prime we just get its (i) as [prime=i*2+1] -> [13=i*2+1] so call SET0(6) to set number 13 as a prime.

Similarly if GET(6) returns a zero then the number [6*2+1=13] is prime

The lines 6, 7 and 8 are made to make it easy setting or getting a value froma single bit of the array.

List[i>>5] is equivalent to List[i/32] and used to access the element in the array that has the bit number i.

The square root hint has been discussed in a previous post (Fastest algorithm to check if a given number is prime)

Initially we set all numbers to be primes (line 13), then set number 1 as not a prime (line 28). then looping through numbers and if there is a prime number we call function FillList() to fill all its multiples as a non prime numbers

Watch the following animation will help you imagine it