[Tapas Bose]

Sieve of Atkin is an optimized version of the ancient Sieve of Eratosthenes. Is this implies that execution time of Sieve of Atkin is faster than that of Sieve of Eratosthenes?
I wrote these two algorithm in C++.
Sieve of Atkin :

#include <new>
#include <cmath>
#include <iostream>
#include <exception>
#include <sys/time.h>

typedef unsigned long long ulong;

using namespace std;

class SieveOfAtkin
{
	private :
		ulong Limit;
		static bool *IsPrime;
		
	public :
		SieveOfAtkin() { }
		SieveOfAtkin(ulong Limit) : Limit(Limit)
		{			
			try
			{
				IsPrime = new bool[Limit];
			}
			catch(bad_alloc&)
			{
				cerr << "  Memory allocation unsuccessful." << endl;
			}
		}
		
		static long NoOfPrime;
		
		void PrimeGen();
		void Display();
		
};

bool *SieveOfAtkin :: IsPrime;
long SieveOfAtkin :: NoOfPrime = 0;

void SieveOfAtkin :: PrimeGen()
{
	ulong SqrtLimit = (ulong)ceil(sqrt((double)Limit));
	
	for (ulong x = 1; x <= SqrtLimit; x++)

        {
        	for (ulong y = 1; y <= SqrtLimit; y++)
        	{
        		ulong n = 4 * x * x + y * y;
        		
        		if (n <= Limit && (n % 12 == 1 || n % 12 == 5))
        		{
        			IsPrime[n] ^= true;
        		}
        		
        		n = 3 * x * x + y * y;



                	if (n <= Limit && n % 12 == 7)

                	{

                    		IsPrime[n] ^= true;

                	}



                	n = 3 * x * x - y * y;



                	if (x > y && n <= Limit && n % 12 == 11)

                	{

                    		IsPrime[n] ^= true;

                	}

            	}

        }
        
        for (ulong n = 5; n <= SqrtLimit; n++)

        {

		if (IsPrime[n])

            	{

                	ulong s = n * n;



                	for (ulong k = s; k <= Limit; k += s)

                	{

                    		IsPrime[k] = false;

                	}

            	}

        }
        
        IsPrime[2] = true;
        IsPrime[3] = true;
}

void SieveOfAtkin :: Display()
{	
	int Counter = 0;
	
	for (ulong i = 2; i <= Limit; i++)
	{
		if (IsPrime[i])
		{
			/*cout << i << "\t";
			
			Counter++;
		
			if (Counter == 19)
			{
				cout << endl;
				Counter = 0;
			}*/
			
			NoOfPrime++;
		}
	}
	
	cout << endl << endl;
}

unsigned GetTickCount()
{
        struct timeval tv;
        
        if (gettimeofday(&tv, NULL) != 0)
        {
                return 0;
        }

        return (tv.tv_sec * 1000) + (tv.tv_usec / 1000);
}

int main ()
{
	unsigned StartTime = GetTickCount();
	
	SieveOfAtkin Object = SieveOfAtkin(10000000LLU);
	
	Object.PrimeGen();
	
	Object.Display();

	cout << "  Number of Prime : " << Object.NoOfPrime << endl;
	cout << "  Execution time : " << (GetTickCount() - StartTime) << " ms." << endl << endl;
	
	return 0;
}

Output :
tapas@My-Child:~/Prime Number Program/Atkin$ g++ Atkin.cpp -o Atkin.o
tapas@My-Child:~/Prime Number Program/Atkin$ ./Atkin.o

Number of Prime : 664579
Execution time : 1370 ms.

And Sieve of Eratosthenes :

#include <new>
#include <list>
#include <cmath>
#include <iostream>
#include <exception>
#include <sys/time.h>

typedef unsigned long long ulong;

using namespace std;

class SieveOfEratosthenes
{
	private :
		ulong Limit;
		static bool *IsPrime;
	
	public :
		SieveOfEratosthenes() { }
		SieveOfEratosthenes(ulong Limit) : Limit(Limit)
		{			
			try
			{
				IsPrime = new bool[Limit + 1];
			}
			catch(bad_alloc&)
			{
				cerr << "  Memory allocation unsuccessful." << endl;
			}
		}
		
		static long NoOfPrime;
		
		void PrimeGen();
		void Display();
};

bool *SieveOfEratosthenes :: IsPrime;
long SieveOfEratosthenes :: NoOfPrime = 0;

void SieveOfEratosthenes :: PrimeGen()
{
	for (ulong i = 2; i*i <= Limit; i++) 
        {
        	if (!IsPrime[i]) 
            	{
                	for (ulong j = i; i*j <= Limit; j++)
                	{
                    		IsPrime[i*j] = true;
                	}
            	}
        }
}

void SieveOfEratosthenes :: Display()
{
	int Counter = 0;
	
	for (ulong i = 2; i <= Limit; i++)
	{
		if (!IsPrime[i])
		{
			/*cout << i << "\t";
			
			Counter++;
		
			if (Counter == 19)
			{
				cout << endl;
				Counter = 0;
			}*/
			
			NoOfPrime++;
		}
	}
	
	cout << endl << endl;
}

unsigned GetTickCount()
{
        struct timeval tv;
        
        if (gettimeofday(&tv, NULL) != 0)
        {
                return 0;
        }

        return (tv.tv_sec * 1000) + (tv.tv_usec / 1000);
}

int main ()
{
	unsigned StartTime = GetTickCount();
	
	SieveOfEratosthenes Object = SieveOfEratosthenes(10000000LLU);
	
	Object.PrimeGen();
	
	Object.Display();

	cout << "  Number of Prime : " << Object.NoOfPrime << endl;
	cout << "  Execution time : " << (GetTickCount() - StartTime) << " ms." << endl << endl;
	
	return 0;
}

Output :
tapas@My-Child:~/Prime Number Program/Eratosthenes$ g++ Eratosthenes.cpp -o Eratosthenes.o
tapas@My-Child:~/Prime Number Program/Eratosthenes$ ./Eratosthenes.o

Number of Prime : 664579
Execution time : 696 ms.

So which one is faster? I can not understand what is going on. Is there any problem in my code? Please help.

[PlasticineGuy]

Why are you using a class? In my opinion classes should not be used when you do not need operator overloading or more than one instace of something.

[Ferencn]

Apparently your implementation of Atkins is slower than your implementation of Eratosthenes...
In the wikipedia, description of the pseudocode, I found this line:

Quote

This pseudocode is written for clarity. Repeated and wasteful calculations mean that it would run slower than the sieve of Eratosthenes.

Perhaps that is also the cause of your findings.

EDIT: Optimiziation possibility: You could take out the 2nd
n = 3 * x * x + y * y;
out of PrimeGen(), because x nor y changed, and a recalculation is unecessary.

This post has been edited by Ferencn: 01 March 2010 - 03:49 AM

[Tapas Bose]

PlasticineGuy, on 01 March 2010 - 02:37 AM, said:

Why are you using a class? In my opinion classes should not be used when you do not need operator overloading or more than one instace of something.

There is a little history of using class. Actually I am pursuing Masters of computer application. In our college, to generate prime we set the upper limit as upperlimit/2 with huge time complexity. I am the first student who has great interest in prime number and its program. I told my teacher about these two prime sieve. He ask me to show the code and he like to use class in his every program. That is why I use class.

Ferencn, on 01 March 2010 - 02:44 AM, said:

Apparently your implementation of Atkins is slower than your implementation of Eratosthenes...
In the wikipedia, description of the pseudocode, I found this line:

Quote

This pseudocode is written for clarity. Repeated and wasteful calculations mean that it would run slower than the sieve of Eratosthenes.

Perhaps that is also the cause of your findings.

EDIT: Optimiziation possibility: You could take out the 2nd
n = 3 * x * x + y * y;
out of PrimeGen(), because x nor y changed, and a recalculation is unecessary.

If I omit n = 3 * x * x + y * y; from the code then the output is :
tapas@My-Child:~/Prime Number Program/Atkin$ g++ Atkin.cpp -o Atkin.o
tapas@My-Child:~/Prime Number Program/Atkin$ ./Atkin.o

Number of Prime : 517615
Execution time : 968 ms.

Which is wrong according to this site here.

This post has been edited by Tapas Bose: 01 March 2010 - 03:55 AM

[Ferencn]

EDIT:
I misread the 2nd! missed the minus sign...

                for (ulong y = 1; y <= SqrtLimit; y++)
                {

                        ulong xsquared = x * x;
                        ulong ysquared = y * y;
                        ulong n = (xsquared + ysquared)<<2 ;
                        
                        if (n <= Limit && (n % 12 == 1 || n % 12 == 5))
                        {
                                IsPrime[n] ^= true;
                        }
                        
                        n -= (xsqared+ysquared) ;

                        if (n <= Limit && n % 12 == 7)

                        {

                                IsPrime[n] ^= true;

                        }
                       n = ((xsquared-ysquared)<<1) +  (xsquared-ysquared);


                        if (x > y && n <= Limit && n % 12 == 11)

                        {

                                IsPrime[n] ^= true;

                        }

                }

        }

And see what it does?

This post has been edited by Ferencn: 01 March 2010 - 04:05 AM

[Tapas Bose]

@Ferencn:
Output :
tapas@My-Child:~/Prime Number Program/Atkin$ g++ Atkin.cpp -o Atkin.o
tapas@My-Child:~/Prime Number Program/Atkin$ ./Atkin.o

Number of Prime : 2
Execution time : 618 ms.

[Martyn.Rae]

Your SieveOfEratosthenes is incorrect. Here is the corrected code

class SieveOfEratosthenes 
{ 
    private : 
        ulong Limit; 
        static bool *IsPrime; 
         
    public : 
        SieveOfEratosthenes() { } 
        SieveOfEratosthenes(ulong Limit) : Limit(Limit) 
        {                        
            try 
            { 
                IsPrime = new bool[Limit + 1];
	        for (ulong i = 2; i*i <= Limit; i++) IsPrime[i] = true;
            } 
            catch(bad_alloc&) 
            { 
                cerr << "  Memory allocation unsuccessful." << endl; 
            } 
       } 
                 
       static long NoOfPrime; 
                
       void PrimeGen(); 
       void Display(); 
}; 
 
bool *SieveOfEratosthenes :: IsPrime; 
long SieveOfEratosthenes :: NoOfPrime = 0; 
 
void SieveOfEratosthenes :: PrimeGen() 
{ 
    for (ulong i = 2; i*i <= Limit; i++)  
    { 
        if (IsPrime[i])  
        { 
            for (ulong j = i+1; j <= Limit; j++) 
            {
	        if ( j%i == 0 ) IsPrime[j] = false; 
            } 
        } 
    } 
} 
 

Try running the correct code and you will have to wait a fair old time. As for the SieveOfAtkin class, this does not give the same number as the SieveOfEratosthenes. I have checked the SieveOfEratosthenes and it is absolutely correct, so check the SieveOfAtkin code as well.

[Ferencn]

Sorry again. I must be getting old :-)
i misread 4*x*x+y*y as 4 (x*x+y*y) and obviously that is wrong...

        for (ulong y = 1; y <= SqrtLimit; y++)
        {
            ulong xsquared = x*x;
            ulong ysquared = y*y;

            ulong n = (xsquared<<2)+ysquared; // n = 4*x*x + y*y

            if (n <= Limit && (n % 12 == 1 || n % 12 == 5))
            {
                IsPrime[n] ^= true;
            }

            n -= xsquared;  // n = 3*x*x + y*y
            if (n <= Limit && n % 12 == 7)
            {
                IsPrime[n] ^= true;
            }

            if (x>y)
            {
                // ONLY DO calculation WHEN X>Y!!
                n -= (ysquared<<1);  // n = 3*x*x - y*y
                if (n <= Limit && n % 12 == 11)
                {
                    IsPrime[n] ^= true;
                }
            }
        }
    }

[Ferencn]

Martyn.Rae, on 01 March 2010 - 11:14 AM, said:

Your SieveOfEratosthenes is incorrect. Here is the corrected code

void SieveOfEratosthenes :: PrimeGen() 
{ 
    for (ulong i = 2; i*i <= Limit; i++)  
    { 
        if (IsPrime[i])  
        { 
            for (ulong j = i+1; j <= Limit; j++) 
            {
	        if ( j%i == 0 ) IsPrime[j] = false; 
            } 
        } 
    } 
} 

Isn't this a bit faster?

void SieveOfEratosthenes :: PrimeGen() 
{ 
    for (ulong i = 2; i*i <= Limit; i++)  
    { 
        if (IsPrime[i])  
        { 
            for (ulong j = i*2; j <= Limit; j+=i) 
            {
	        IsPrime[j] = false; 
            } 
        } 
    } 
} 

This post has been edited by Ferencn: 01 March 2010 - 04:24 AM

[Martyn.Rae]

I forgot to include the modified Display routine as well.

void SieveOfEratosthenes :: Display() 
{ 
        int Counter = 0; 
         
        for (ulong i = 2; i <= Limit; i++) 
        { 
                if (IsPrime[i]) 
                { 
                        /*cout << i << "\t"; 
                         
                        Counter++; 
                 
                        if (Counter == 19) 
                        { 
                                cout << endl; 
                                Counter = 0; 
                        }*/ 
                         
                        NoOfPrime++; 
                } 
        } 
         
        cout << endl << endl; 
} 

@Ferencn that getting old thing is catching ... wait a minute though ... I am old!!!

This post has been edited by Martyn.Rae: 01 March 2010 - 04:25 AM

[Ferencn]

This code:

void SieveOfEratosthenes :: Display()
{
        int Counter = 0;
        
        for (ulong i = 2; i <= Limit; i++)
        {
                if (!IsPrime[i])
                {

seems to be counting non-primes!

Well ,maybe old, but your eyes were quicker than mine on the counting part ;-)

[Martyn.Rae]

Ferencn, on 01 March 2010 - 10:30 AM, said:

This code:

void SieveOfEratosthenes :: Display()
{
        int Counter = 0;
        
        for (ulong i = 2; i <= Limit; i++)
        {
                if (!IsPrime[i])
                {

seems to be counting non-primes!

See my previous post!

[Tapas Bose]

@ Martyn.Rae. Sir I think the initialization of IsPrime array with 'true' is unnecessary because I declare it as static since static declaration initialize the array by default by 'false' and it increase the time complexity. That is why I use if(!IsPrime[i]) instead of if(!IsPrime[i]) and my method also checked the multiples of i. Am I wrong sir?

This post has been edited by Tapas Bose: 01 March 2010 - 04:36 AM

[Tapas Bose]

I just change the PrimeGen function as :

void SieveOfEratosthenes :: PrimeGen()
{
	for (ulong i = 2; i*i <= Limit; i++) 
        {
        	if (!IsPrime[i]) 
            	{
                	for (ulong j = i*2; j <= Limit; j+=i) 
            		{
                		IsPrime[j] = true; 
            		} 
            	}
        }
}

The output is:
tapas@My-Child:~/Prime Number Program/Eratosthenes$ g++ Eratosthenes.cpp -o Eratosthenes.o
tapas@My-Child:~/Prime Number Program/Eratosthenes$ ./Eratosthenes.o

Number of Prime : 664579
Execution time : 946 ms.

[Martyn.Rae]

You are correct in saying that static variables are initialized to zero. Unfortunately, you are using a static pointer to an array. This will be set to zero, but then you allocate the array using the new operator. This provides you with a chunk of memory that contains well, garbage!

If you use the same initialization code for both prime number classes, it will not make any difference to the overall timings.