[baavgai]

Ferencn, on 01 March 2010 - 09:59 AM, said:

@ baavgai: Atkin was beaten by Eratosthenes from the beginning!
Martyn and I already found the improvements that you propose for Eratosthenes...

Granted. I don't really care about the math; I want the code to be done well. Presenting both lack of static and the idea of using a base class or interface as a test harness for this sort of thing.

I tried to squeeze a little more out of Atkin, but couldn't do it. I don't much care for it. Eratosthenes is intuitive and elegant. Atkin is a mess of pretentious quadratic diddling that fails to deliver.

[Ferencn]

baavgai, on 01 March 2010 - 05:05 PM, said:

Eratosthenes is intuitive and elegant. Atkin is a mess of pretentious quadratic diddling that fails to deliver.

I agree. 100%.

[Tapas Bose]

Here is a link.

[baavgai]

Tapas Bose, on 01 March 2010 - 11:47 AM, said:

Here is a link.

Thanks. Saw the site on wiki, but the paper was locked up. For completeness, looked at the source code ( http://cr.yp.to/primegen.html ). It did run very quickly:

Finished L=16580609: 11029697.000000
50847534 primes up to 1000000000.
Timings are in ticks. Nanoseconds per tick: approximately 0.655067.
Overall seconds: approximately 1.633908.

However, it's not a simple algorithm. It's cached and precalculated; a lot. You'll find little gems like this in the code:

static const uint32 two[32] = {
  0x00000001 , 0x00000002 , 0x00000004 , 0x00000008
, 0x00000010 , 0x00000020 , 0x00000040 , 0x00000080
, 0x00000100 , 0x00000200 , 0x00000400 , 0x00000800
, 0x00001000 , 0x00002000 , 0x00004000 , 0x00008000
, 0x00010000 , 0x00020000 , 0x00040000 , 0x00080000
, 0x00100000 , 0x00200000 , 0x00400000 , 0x00800000
, 0x01000000 , 0x02000000 , 0x04000000 , 0x08000000
, 0x10000000 , 0x20000000 , 0x40000000 , 0x80000000
} ;

/* squares of primes >= 7, < 240 */
uint32 qqtab[49] = {
  49,121,169,289,361,529,841,961,1369,1681
 ,1849,2209,2809,3481,3721,4489,5041,5329,6241,6889
 ,7921,9409,10201,10609,11449,11881,12769,16129,17161,18769
 ,19321,22201,22801,24649,26569,27889,29929,32041,32761,36481
 ,37249,38809,39601,44521,49729,51529,52441,54289,57121
} ;

There are also a lot of "register" declarations in that old C code.

Overall, still unimpressed. I'm reasonably certain your implementation was correct. I don't believe their fundamental algorithm, as stated, is optimized in any meaningful way. Or truly reflective of the code they used as a proof.

[Tapas Bose]

baavgai, on 01 March 2010 - 10:19 AM, said:

Tapas Bose, on 01 March 2010 - 11:47 AM, said:

Here is a link.

Thanks. Saw the site on wiki, but the paper was locked up. For completeness, looked at the source code ( http://cr.yp.to/primegen.html ). It did run very quickly:

Finished L=16580609: 11029697.000000
50847534 primes up to 1000000000.
Timings are in ticks. Nanoseconds per tick: approximately 0.655067.
Overall seconds: approximately 1.633908.

However, it's not a simple algorithm. It's cached and precalculated; a lot. You'll find little gems like this in the code:

static const uint32 two[32] = {
  0x00000001 , 0x00000002 , 0x00000004 , 0x00000008
, 0x00000010 , 0x00000020 , 0x00000040 , 0x00000080
, 0x00000100 , 0x00000200 , 0x00000400 , 0x00000800
, 0x00001000 , 0x00002000 , 0x00004000 , 0x00008000
, 0x00010000 , 0x00020000 , 0x00040000 , 0x00080000
, 0x00100000 , 0x00200000 , 0x00400000 , 0x00800000
, 0x01000000 , 0x02000000 , 0x04000000 , 0x08000000
, 0x10000000 , 0x20000000 , 0x40000000 , 0x80000000
} ;

/* squares of primes >= 7, < 240 */
uint32 qqtab[49] = {
  49,121,169,289,361,529,841,961,1369,1681
 ,1849,2209,2809,3481,3721,4489,5041,5329,6241,6889
 ,7921,9409,10201,10609,11449,11881,12769,16129,17161,18769
 ,19321,22201,22801,24649,26569,27889,29929,32041,32761,36481
 ,37249,38809,39601,44521,49729,51529,52441,54289,57121
} ;

There are also a lot of "register" declarations in that old C code.

Overall, still unimpressed. I'm reasonably certain your implementation was correct. I don't believe their fundamental algorithm, as stated, is optimized in any meaningful way. Or truly reflective of the code they used as a proof.

I have this primegen program but I can not understand the code.

[Tapas Bose]

Here is another Sieve, this is known as Sieve of Sundaram. Here is the code :

void SieveOfSundaram :: PrimeGen()
{
	ulong iLim = Limit / 2;
	
	for (ulong i = 1; i < iLim; i++)
	{
		for (ulong j = 1; j <= i; j++)
		{
			ulong n = i + j + 2 * i * j;
			
			if (n >= Limit)
			{
				break;
			}
			
			IsPrime[n] = false;
		}
	}
}

Output :
Total number of prime : 664578 (except 2)
Execution time : 913 ms.

Is there any optimization possible?

This post has been edited by Tapas Bose: 01 March 2010 - 01:02 PM

[Tellus]

Sorry to BUMP this post but:
I made a test to GENERATE (not COUNT) the prime numbers up to 70 000 000:

Sieve ot Atkin:

#include<stdio.h>
#include<time.h>
#include<iostream>
#define MAXN 70000000
#include<math.h>
using namespace std;
unsigned long sieve[MAXN+1];


void PrimeFind(unsigned long Limit)
{           double sqrtLimit=sqrt(Limit)+1;
            unsigned long xsqd = 1;
            unsigned long xadd = 3;
            for (unsigned long x = 1; x<=sqrtLimit; x++)

            {       
                unsigned long ysqd = 1;
                unsigned long yadd = 3;
                for (unsigned long y = 1; y<=sqrtLimit; y++)
                {
                        

                        unsigned long n = (xsqd << 2) + ysqd; // n = 4*x*x + y*y

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

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

                        if (x > y)
                        {
                                n -= (ysqd << 1);  // n = 3*x*x - y*y
                                
                                if (n <= Limit && n % 12 == 11)
                                {
                                        sieve[n] ^= true;
                                }
                        }
                        ysqd+=yadd;
                        yadd+=2;
                }
                 xsqd+=xadd;
                 xadd+=2;
        }                 
}
int main()
{
 clock_t start = clock();
 PrimeFind(MAXN);
 sieve[2]=1;sieve[3]=1;
 unsigned long sqrtLimit=(unsigned long)sqrt(MAXN)+1;
 for(unsigned long i=5 ; i<=sqrtLimit ; i++){
              if(sieve[i]){
                           unsigned long j=i*i;
                           for (unsigned long k = j; k <= MAXN; k += j)
                                {
                                sieve[k] = false;
                                }
                           }
              }   
 printf("%10lu", 2);//printing the only even prime
 for(unsigned long i=1 ; i<=MAXN ; i+=2)// without checking for odds
 if(sieve[i])printf("%10lu", i);
 printf("\nTime elapsed: %f\n", ((double)clock() - start) / CLOCKS_PER_SEC);
 system("pause");
 return 0;
}

Sieve of Eratosthenes:

#include<stdio.h>
#define MAXN 70000000
#include <iostream>
#include <time.h>
const unsigned long n=MAXN;
char sieve[MAXN];

void eratosten (unsigned long n)
{
  unsigned long j, i=2;
  while (i <= n ){
       if(sieve[i] == 0){
            printf("%10lu", i);
            j=i*i;
            while (j <= n) {
              sieve[j] = 1;
              j += i;
            }
       }
       i++;
  }
}
int main(void) {
    clock_t start = clock();
    eratosten( n );
    printf("\nTime elapsed: %f\n", ((double)clock() - start) / CLOCKS_PER_SEC);
    system("pause");
    return 0;
}

I am using pure C and the only difference with baavgai's solution on the Sieve of Eratosthenes is that I start the actual sieving from i*i not from i+i (which speeds it up a bit).

THE RESULTS:
Sieve of Eratosthenes: 197.968 sec
Sieve ot Atkin: 197.828 sec

This post has been edited by Tellus: 17 April 2010 - 05:52 AM

[baavgai]

Sorry, your implementation is NOT pure C ( #include<iostream> ). It could be C99, but even then you'd need stdbool.h.

Also, printing while you figure things out will skew your results. And why two programs?

I was curious, so I gave what you had a spin. I rewrote it to be consistent.

#include <stdlib.h>
#include <stdio.h>
#include <time.h>
#include <math.h>

#define BOOL int
#define TRUE 1
#define FALSE 0

// #define SHOW_OUTPUT

typedef void(*PrimeFunc)(BOOL [], unsigned long);

void fillSieve(BOOL sieve[], unsigned long size, BOOL value) {
	if (size==0) { return; }
	do { sieve[--size] = value; } while(size!=0);
}

void printSieve(BOOL sieve[], unsigned long size) {
#ifdef SHOW_OUTPUT
	unsigned long i = 2;
	BOOL value = sieve[2];
	printf("%10lu", i++);
	for(; i<=size; i+=2) {
		if(sieve[i]==value) {
			printf("%10lu", i);
		}
	}
#endif
}


void PrimeFind(BOOL sieve[], unsigned long Limit) {
	double sqrtLimit=sqrt(Limit)+1;
	unsigned long x, xsqd = 1, xadd = 3;
	for (x = 1; x<=sqrtLimit; x++) {       
		unsigned long y, ysqd = 1, yadd = 3;
		for (y = 1; y<=sqrtLimit; y++) {
			unsigned long n = (xsqd << 2) + ysqd; // n = 4*x*x + y*y
			if (n <= Limit && (n % 12 == 1 || n % 12 == 5)) { sieve[n] ^= TRUE; }
			n -= xsqd;  // n = 3*x*x + y*y
			if (n <= Limit && n % 12 == 7) { sieve[n] ^= TRUE; }
			if (x > y) {
				n -= (ysqd << 1);  // n = 3*x*x - y*y
				if (n <= Limit && n % 12 == 11) { sieve[n] ^= TRUE; }
			}
			ysqd+=yadd;
			yadd+=2;
		}
		xsqd+=xadd;
		xadd+=2;
	}
}

void sieveOfAtkin(BOOL sieve[], unsigned long MAXN) { 
	unsigned long i, sqrtLimit=(unsigned long)sqrt(MAXN)+1;
	
	fillSieve(sieve, MAXN, FALSE);
	PrimeFind(sieve, MAXN);
	sieve[2]=TRUE;
	sieve[3]=TRUE;
	for(i=5 ; i<=sqrtLimit ; i++){
		if(sieve[i]){
			unsigned long k, j=i*i;
			for (k = j; k <= MAXN; k += j) { sieve[k] = FALSE; }
		}
	}
	printSieve(sieve, MAXN);
}



void sieveOfEratosthenes(BOOL sieve[], unsigned long n) {
	unsigned long j, i=2;
	
	fillSieve(sieve, n, FALSE);
	while (i <= n ) {
		if(!sieve[i]) {
			j=i*i;
			while (j <= n) {
				sieve[j] = TRUE;
				j += i;
			}
		}
	i++;
	}
	printSieve(sieve, n);
}


void timeAndShow(const char *name, PrimeFunc func, BOOL sieve[], unsigned long n) {
	clock_t start;
	
	printf("Testing: %s\n", name);
	start = clock();
	func(sieve, n);
	printf("\nTime elapsed: %f\n", ((double)clock() - start) / CLOCKS_PER_SEC);
}

void process(unsigned long size) {
	BOOL *sieve = malloc(sizeof(BOOL)*size);
	
	printf("Sample Size: %lu\n", size);
	
	timeAndShow("Sieve of Eratosthenes", sieveOfEratosthenes, sieve, size);
	timeAndShow("Sieve of Atkin", sieveOfAtkin, sieve, size);

	free(sieve);
}


int main(void) {
	process(70000000L);
}

Results:

Sample Size: 1000000
Sieve of Eratosthenes: 0.060000
Sieve of Atkin: 0.080000


Sample Size: 70000000
Sieve of Eratosthenes: 5.660000
Sieve of Atkin:: 5.900000

Conclusion: Atkin still sucks.

[Tellus]

Ok, I made some tests of my own again. First, I want to apologise - it seems that starting from i*i instead of i*2 is a mistake (Erathostenes omits primes that way). I used <iostream> only for the system("pause") and it can be replaced with <stdlib.h> ( or <dos.h> in some cases ). So I am starting the sieving from i*2 this time :

#include<stdio.h>
#define MAXN  220000000
#include <stdlib.h>
#include <time.h>
#include <math.h>
const unsigned long n=MAXN;
char sieveE[MAXN];
char sieveA[MAXN];

void EratosthenesSieve (unsigned long n);
void AtkinSieve(unsigned long Limit);

int main(void) {
   
   
    AtkinSieve( n ) ;
    printf("\n");
    EratosthenesSieve( n );
    system("pause");    
    return 0;
}


void EratosthenesSieve (unsigned long n)
{ clock_t start = clock();
  unsigned long j, i=2, br=0;
  while (i <= n ){
       if(sieveE[i] == 0){
            br++;
            for(j=i*2; j<= n ; j+=i)
              sieveE[j] = 1;
              
            
       }
       i++;
  }
printf("Eratosthenes counted: %10lu", br);
printf("\nTime elapsed: %f\n", ((double)clock() - start) / CLOCKS_PER_SEC);
}

void AtkinSieve(unsigned long Limit)
{           clock_t start = clock();
            double sqrtLimit=sqrt(Limit)+1;
            unsigned long xsqd = 1;
            unsigned long xadd = 3;
            for (unsigned long x = 1; x<=sqrtLimit; x++)

            {       
                unsigned long ysqd = 1;
                unsigned long yadd = 3;
                for (unsigned long y = 1; y<=sqrtLimit; y++)
                {
                        

                        unsigned long n = (xsqd << 2) + ysqd; // n = 4*x*x + y*y

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

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

                        if (x > y)
                        {
                                n -= (ysqd << 1);  // n = 3*x*x - y*y
                                
                                if (n <= Limit && n % 12 == 11)
                                {
                                        sieveA[n] ^= true;
                                }
                        }
                        ysqd+=yadd;
                        yadd+=2;
                }
                 xsqd+=xadd;
                 xadd+=2;
        }                 
        unsigned long br=0;             
        sieveA[2]=1;sieveA[3]=1;
        for(unsigned long i=5 ; i<=sqrtLimit ; i++){
              if(sieveA[i]){
                           unsigned long j=i*i;
                           for (unsigned long k = j; k <= MAXN; k += j)
                                {
                                sieveA[k] = false;
                                }
                           }
        }           
        for(unsigned long i=1 ; i<=MAXN ; i+=2)// without checking for odds
        if(sieveA[i])br++;
        printf("Atkin counted: %10lu", br+1);
        printf("\nTime elapsed: %f\n", ((double)clock() - start) / CLOCKS_PER_SEC);
        
}

And the results:

Up to 1000000:
Eratosthenes counted:      78498
Time elapsed: 0.014000

Atkin counted:      78498
Time elapsed: 0.018000

Up to 10000000:
Eratostenes counted:     664579
Time elapsed: 0.308000

Atkin counted:     664579
Time elapsed: 0.242000

Up to 220000000:
Eratosthenes counted:   12122540
Time elapsed: 11.174000

Atkin counted:   12122540
Time elapsed: 8.594000

It seems that Atkin gets faster when testing with larger numbers (: .

This post has been edited by Tellus: 18 April 2010 - 01:57 AM

[baavgai]

Tellus, on 18 April 2010 - 02:33 AM, said:

It seems that Atkin gets faster when testing with larger numbers (: .

Interesting. You understand, of course, that I had to do my own test.

Test program looks like so:

#include <stdlib.h>
#include <stdio.h>
#include <time.h>
#include <math.h>
#include <stdbool.h>

// gcc sieve4.c -lm -std=c99
typedef unsigned long ulong;

typedef void(*PrimeFunc)(bool [], ulong);

typedef struct {
	ulong found;
	double elapsedTime;
} TestResults;

typedef struct {
	char *name;
	PrimeFunc func;
	bool fillValue;
} TestItem;

void fillSieve(bool sieve[], ulong size, bool value) {
	if (size==0) { return; }
	do { sieve[--size] = value; } while(size!=0);
}

ulong countSieve(bool sieve[], ulong size) {
	ulong i, count = 1;
	bool value = sieve[2];
	for(i=3; i<=size; i+=2) {
		if(sieve[i]==value) { count++; }
	}
	return count;
}

void sieveOfAtkin(bool sieveA[], ulong Limit) {
	double sqrtLimit=sqrt(Limit)+1;
	ulong x, xsqd = 1, xadd = 3;
	for (x = 1; x<=sqrtLimit; x++) {       
		 ulong y, ysqd = 1, yadd = 3;
		 for (y = 1; y<=sqrtLimit; y++) {
			ulong n = (xsqd << 2) + ysqd; // n = 4*x*x + y*y
			if (n <= Limit && (n % 12 == 1 || n % 12 == 5)) { sieveA[n] ^= true; }
			n -= xsqd;  // n = 3*x*x + y*y
			if (n <= Limit && n % 12 == 7) { sieveA[n] ^= true; }
			if (x > y) { 
				n -= (ysqd << 1);  // n = 3*x*x - y*y
				if (n <= Limit && n % 12 == 11){ sieveA[n] ^= true; }
			}
			ysqd+=yadd;
			yadd+=2;
		}
		xsqd+=xadd;
		xadd+=2;
	}    
	
	sieveA[2]=1; sieveA[3]=1;
	for(ulong i=5 ; i<=sqrtLimit ; i++){
		if(sieveA[i]){
			ulong k, j=i*i;
			for (ulong k = j; k < Limit; k += j) { sieveA[k] = false; }
		}
	}
}




void sieveOfEratosthenes(bool sieve[], ulong n) {
	ulong j, i=2;
	
	while (i <= n ) {
		if(!sieve[i]) {
			j=i*i;
			while (j <= n) {
				sieve[j] = true;
				j += i;
			}
		}
	i++;
	}
}

void sieveOfEratosthenes2(bool sieve[], ulong size) {
	for (ulong i = 2; i*i <= size; i++) {
		if (sieve[i]) { 
			for (ulong j = i+i; j <=size; j += i) { sieve[j] = false; } 
		} 
	} 
}

void sieveOfEratosthenes3(bool sieve[], ulong size) {
	ulong limit = sqrt(size)+1;
	for (ulong i = 2; i <= limit; i++) {
		if (sieve[i]) { 
			for (ulong j = i+i; j <=size; j += i) { sieve[j] = false; } 
		} 
	} 
}



double timeAndShow(PrimeFunc func, bool sieve[], ulong size, bool fillValue) {
	clock_t start;
	
	fillSieve(sieve, size, fillValue);
	start = clock();
	func(sieve, size);
	return ((double)clock() - start) / CLOCKS_PER_SEC;
}


void processForSample(TestItem list[], int listSize, ulong size, TestResults *resultList) {
	bool *sieve = malloc(sizeof(bool)*(size+1));
	int i;
	
	printf("Sample Size: %lu\n", size);
	for(i=0; i<listSize; i++) {
		TestResults results;
		TestItem *item = &list[i];
		
		printf("     Testing: %s\n", item->name);
		results.elapsedTime = timeAndShow(item->func, sieve, size, item->fillValue);
		results.found = countSieve(sieve, size);
		printf("Time elapsed: %f\n", results.elapsedTime);
		printf("       Found: %lu\n", results.found);
		resultList[i] = results;
	}
	printf("\n");
	
	free(sieve);
}


void process() {
	TestItem list[] = {
		{"Era", sieveOfEratosthenes, false},
		{"Era 2", sieveOfEratosthenes2, true},
		{"Era 3", sieveOfEratosthenes3, true},
		{"Atkin", sieveOfAtkin, false},
	};
	ulong samples[] = {100000, 1000000, 10000000, 70000000, 220000000};
	TestResults **resultList;
	int listSize = sizeof(list)/sizeof(TestItem);
	int samplesSize = sizeof(samples)/sizeof(int);
	
	resultList = malloc(sizeof(TestResults *) * samplesSize);
	
	for(int i=0; i<samplesSize; i++) {
		resultList[i] = malloc(sizeof(TestResults) * listSize);
		processForSample(list, listSize, samples[i], resultList[i]);
	}

	printf("%12s %12s ", " ", " ");
	for(int i=0; i<listSize; i++) { printf("%12s ", list[i].name); }
	printf("\n");
	
	for(int i=0; i<samplesSize; i++) {
		TestResults *results = resultList[i];
		
		printf("%12lu %12s ", samples[i], "found");
		for(int j=0; j<listSize; j++) { printf("%12lu ", results[j].found); }
		printf("\n");
		
		printf("%12s %12s ", " ", "time (secs)");
		for(int j=0; j<listSize; j++) { printf("%12.5f ", results[j].elapsedTime); }
		printf("\n");
		
		free(results);
	}
	
	free(resultList);
}


int main(void) {
	process();
}

Results:

                                   Era        Era 2        Era 3        Atkin 
      100000        found         9592         9592         9592         9592 
              time (secs)      0.00000      0.00000      0.00000      0.01000 
     1000000        found        78491        78498        78498        78498 
              time (secs)      0.02000      0.03000      0.03000      0.05000 
    10000000        found       664116       664579       664579       664579 
              time (secs)      0.46000      0.46000      0.45000      0.71000 
    70000000        found      4096478      4118064      4118064      4118064 
              time (secs)      3.71000      3.66000      3.65000      5.27000 
   220000000        found     11912035     12122540     12122540     12122540 
              time (secs)     12.55000     12.19000     12.20000     16.79000 

That first Era is the one you already spotted the bug for. Era 2 and 3 are how I'd write it, using the sqrt or not.

Hmm, I want a better display. I'm going to prefer Era 2, the one without the sqrt, just because I'd prefer not to implement it. It does appear that Atkins begins to catch up on larger sets. How large can I go?

Next test:

void process2() {
	TestItem list[] = {
		{"Era 2", sieveOfEratosthenes2, true},
		{"Atkin", sieveOfAtkin, false},
	};
	ulong samples[] = {10000000, 70000000, 220000000, 1000000000};
	
	TestResults **resultList;
	int listSize = sizeof(list)/sizeof(TestItem);
	int samplesSize = sizeof(samples)/sizeof(int);
	
	resultList = malloc(sizeof(TestResults *) * samplesSize);
	
	for(int i=0; i<samplesSize; i++) {
		resultList[i] = malloc(sizeof(TestResults) * listSize);
		processForSample(list, listSize, samples[i], resultList[i]);
	}

	printf("%12s %12s ", " ", " ");
	for(int i=0; i<listSize; i++) { printf("%12s ", list[i].name); }
	printf("\n");
	
	for(int i=0; i<samplesSize; i++) {
		TestResults *results = resultList[i];
		
		printf("%12lu %12s ", samples[i], "time (secs)");
		for(int j=0; j<listSize; j++) { 
			printf("%12.5f ", results[j].elapsedTime);
		}
		printf("\n");
		
		if (results[0].elapsedTime!=0) {
			printf("%12s %12s ", " ", "pct");
			for(int j=0; j<listSize; j++) { 
				printf("%11.5f%% ", (results[j].elapsedTime/results[0].elapsedTime) * 100.0);
				
			}
			printf("\n");
		}
		
		free(results);
	}
	
	free(resultList);
}

Results:

                                 Era 2        Atkin 
    10000000  time (secs)      0.45000      0.71000 
                      pct   100.00000%   157.77778% 
    70000000  time (secs)      3.66000      5.25000 
                      pct   100.00000%   143.44262% 
   220000000  time (secs)     12.19000     16.81000 
                      pct   100.00000%   137.89992% 
  1000000000  time (secs)     59.49000     78.88000 
                      pct   100.00000%   132.59371% 

Now, 1,000,000,000 is dangeriously close to the end of my ulong. Also, my gain seemed to have slowed down. Scrolling back, I find this:

Sample Size: 220000000
     Testing: Era 2
Time elapsed: 12.190000
       Found: 12122540
     Testing: Atkin
Time elapsed: 16.810000
       Found: 12122540

Sample Size: 1000000000
     Testing: Era 2
Time elapsed: 59.490000
       Found: 50847534
     Testing: Atkin
Time elapsed: 78.880000
       Found: 53327624

Oops, our found count doesn't match at 1000000000. I don't know who's right, but the internet (http://primes.utm.edu/howmany.shtml) tells me good old Eratosthenes is! Looks like the complexity to Atkin rolled the buffer.

So, to go bigger, you'd need a bigint implementation. Perhaps Atkin could catch up, but it would have to be through large number algorithms where every operation is more expensive, so I believe not.

[Tellus]

Quote

Oops, our found count doesn't match at 1000000000

Hm...it is logical that the Atkin will overflow because it's more memory consuming of course, but from math point of view it should be faster and it's odd that it is not presented in your tests. I noticed something in your Era2 tho:

void sieveOfEratosthenes2(bool sieve[], ulong size) {
	for (ulong i = 2; i*i <= size; i++) {
		if (sieve[i]) { 
			for (ulong j = i+i; j <=size; j += i) { sieve[j] = false; } 
		} 
	} 
}

I don't understand why i*i<=size. If primes are being counted that will lead to a significant lower number at least from my point of view. I haven't programmed for a year or so and now that I decided to take it serious, finding prime numbers was one of the simple but significant algorithms to start with. I still don't understand why in your tests Atkin is slower (:

#include<stdio.h>
#define MAXN  900000000
#include <stdlib.h>
#include <time.h>
#include <math.h>
const unsigned long n=MAXN;
char sieveE[MAXN];
char sieveA[MAXN];

void EratosthenesSieve (unsigned long n);
void AtkinSieve(unsigned long Limit);

int main(void) {
   
   
    AtkinSieve( n ) ;
    printf("\n");
    EratosthenesSieve( n );
    system("pause");    
    return 0;
}


void EratosthenesSieve (unsigned long n)
{ clock_t start = clock();
  unsigned long j, i, br=0;
   
  for(i=2; i<=n ; i++)
       if(!sieveE[i]){br++;
            for(j=i+i; j<= n ; j+=i)
              {sieveE[j] = 1;}
       }
                

printf("Eratosthenes counted: %10lu", br);
printf("\nTime elapsed: %f\n", ((double)clock() - start) / CLOCKS_PER_SEC);
}

void AtkinSieve(unsigned long Limit)
{           clock_t start = clock();
            double sqrtLimit=sqrt(Limit)+1;
            unsigned long xsqd = 1;
            unsigned long xadd = 3;
            for (unsigned long x = 1; x<=sqrtLimit; x++)

            {       
                unsigned long ysqd = 1;
                unsigned long yadd = 3;
                for (unsigned long y = 1; y<=sqrtLimit; y++)
                {
                        

                        unsigned long n = (xsqd << 2) + ysqd; // n = 4*x*x + y*y

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

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

                        if (x > y)
                        {
                                n -= (ysqd << 1);  // n = 3*x*x - y*y
                                
                                if (n <= Limit && n % 12 == 11)
                                {
                                        sieveA[n] ^= true;
                                }
                        }
                        ysqd+=yadd;
                        yadd+=2;
                }
                 xsqd+=xadd;
                 xadd+=2;
        }                 
        unsigned long br=0;             
        sieveA[2]=1;sieveA[3]=1;
        for(unsigned long i=5 ; i<=sqrtLimit ; i++){
              if(sieveA[i]){
                           unsigned long j=i*i;
                           for (unsigned long k = j; k <= MAXN; k += j)
                                {
                                sieveA[k] = false;
                                }
                           }
        }           
        for(unsigned long i=1 ; i<=MAXN ; i+=2)// without checking for odds
        if(sieveA[i])br++;
        printf("Atkin counted: %10lu", br+1);
        printf("\nTime elapsed: %f\n", ((double)clock() - start) / CLOCKS_PER_SEC);
        
}

Using for over while this time (which should not lead to any difference but still):

Up to 900000000
Atkin counted: 46229727
Time elapsed: 38.570000

Eratosthenes counted: 46009215
Time elapsed: 51.675000

It seems that Atkin fails here (because of the overflow I suppose) but it is still significantly faster. Weird...

Edit: Also your source gives me errors on 'malloc' and since I don't understand it I can't fix it for now. But oh well I will google it and try my best. (It might be because I use Dev-C++... still not sure):
In function `void processForSample(TestItem*, int, ulong, TestResults*)':
bool *sieve = malloc(sizeof(bool)*(size+1));
invalid conversion from `void*' to `bool*'

In function `void process()':
resultList = malloc(sizeof(TestResults *) * samplesSize);
invalid conversion from `void*' to `TestResults**'

In function `void process()':
resultList[i] = malloc(sizeof(TestResults) * listSize);
invalid conversion from `void*' to `TestResults**'

For the first error if I type it this way:
bool *sieve = (bool *) malloc(sizeof(bool)*(size+1));
there is no error but I don't know if this is an actual fix. About the other errors I have no idea what to do (:

This post has been edited by Tellus: 18 April 2010 - 11:36 AM

[baavgai]

Tellus, on 18 April 2010 - 11:59 AM, said:

I don't understand why i*i<=size.

These are roughly equal: i*i<=size and i<=sqrt(size)+1. You already use that in your code, right? After the square of the size in Eratosthenes it's done. Akins does the same thing.

Tellus, on 18 April 2010 - 11:59 AM, said:

I still don't understand why in your tests Atkin is slower (:

Well, your Eratosthenes imposes has issues. Also, you never seem to init your arrays...

for(i=2; i<=n ; i++) // spin (n-sqrt(n)) longer than required
	if(!sieveE[i]){br++; // incrementing on every hit, marginal my not required
}

Implementing your version like this:

void sieveOfEratosthenes4(bool sieve[], ulong size) {
	ulong br = 0;
	for (ulong i = 2; i <= size; i++) {
		if (sieve[i]) { 
			br++;
			for (ulong j = i+i; j <=size; j += i) { sieve[j] = false; } 
		} 
	} 
}

I get these results.

                                 Era 2        Atkin        Era 4 
   900000000  time (secs)     53.86000     72.86000     89.28000 
                              46009215     46229727     46009215 
                      pct   100.00000%   135.27664%   165.76309% 

Hell of a difference an extra 899970000 iterations makes.

Tellus, on 18 April 2010 - 11:59 AM, said:

Edit: Also your source gives me errors on 'malloc' and since I don't understand it I can't fix it for now. But oh well I will google it and try my best. (It might be because I use Dev-C++... still not sure):

You said pure C. However, you're still using a C++ compiler. They're different bloody languages. C++ is hypocritically pedantic where it comes to data type. There are no void pointers and you have to cast you malloc.

If you want to use C++? Fine. Use it. No malloc, no #define, explict typing, implicit bool, const, OO, classes, all of it. It leads to a different style and approach.

Here's the prior code I did in C++.

#include <iostream>
#include <ctime>
#include <cmath>

typedef unsigned long ulong;

struct PrimeFinder {
	bool *sieve;
	const int size;
	double elapsedTime;
	ulong found;
	
	PrimeFinder(int n) : size(n) { sieve = new bool[size+1]; }
	~PrimeFinder() { delete [] sieve; }
	
	virtual const char *getName() const = 0;
	
	void findPrimes() {
		clock_t start;
		
		fill();
		start = clock();
		generate();
		elapsedTime = ((double)clock() - start) / CLOCKS_PER_SEC;
		found = count();
	}
	
	void print(std::ostream &out) const {
		out << "     Testing: " << getName() << std::endl;
		out << "Time elapsed: " << elapsedTime << std::endl;
		out << "       Found: " << found << std::endl;
		out << std::endl;
	}
	
protected:	
	virtual bool getInitValue() const = 0;
	virtual void generate() = 0;

	void fill() {
		bool value = getInitValue();
		for(ulong i=0; i<=size; i++) { sieve[i]=value; }
	}

	ulong count() const {
		ulong count = 1;
		bool value = sieve[2];
		for(ulong i=3; i<=size; i+=2) {
			if(sieve[i]==value) { count++; }
		}
		return count;
	}
};


struct Atkin : public PrimeFinder {
	const char *getName() const { return "Atkin"; }
	Atkin(int n) : PrimeFinder(n) {}
protected:
	bool getInitValue() const { return false; }
	void generate() {
		bool *sieveA = this->sieve;
		ulong Limit = this->size;
		double sqrtLimit=sqrt(Limit)+1;
		ulong x, xsqd = 1, xadd = 3;
		for (x = 1; x<=sqrtLimit; x++) {       
			 ulong y, ysqd = 1, yadd = 3;
			 for (y = 1; y<=sqrtLimit; y++) {
				ulong n = (xsqd << 2) + ysqd; // n = 4*x*x + y*y
				if (n <= Limit && (n % 12 == 1 || n % 12 == 5)) { sieveA[n] ^= true; }
				n -= xsqd;  // n = 3*x*x + y*y
				if (n <= Limit && n % 12 == 7) { sieveA[n] ^= true; }
				if (x > y) { 
					n -= (ysqd << 1);  // n = 3*x*x - y*y
					if (n <= Limit && n % 12 == 11){ sieveA[n] ^= true; }
				}
				ysqd+=yadd;
				yadd+=2;
			}
			xsqd+=xadd;
			xadd+=2;
		}    
		
		sieveA[2]=1; sieveA[3]=1;
		for(ulong i=5 ; i<=sqrtLimit ; i++){
			if(sieveA[i]){
				ulong k, j=i*i;
				for (ulong k = j; k < Limit; k += j) { sieveA[k] = false; }
			}
		}
	}
};


struct Eratosthenes : public PrimeFinder {
	const char *getName() const { return "Eratosthenes"; }
	Eratosthenes(int n) : PrimeFinder(n) {}
protected:
	bool getInitValue() const { return true; }
	void generate() {
		ulong limit = sqrt(size)+1;
		for (ulong i = 2; i <= limit; i++) {
			if (sieve[i]) { 
				for (ulong j = i+i; j <=size; j += i) { sieve[j] = false; } 
			} 
		} 
	}
};


void run(PrimeFinder &pf) {
	pf.findPrimes();
	pf.print(std::cout);
}

void process(ulong sampleSize) {
	Eratosthenes s1(sampleSize);
	Atkin s2(sampleSize);

	std::cout << " Sample Size: " << sampleSize << std::endl;
	
	run(s1);
	run(s2);
	
	std::cout << std::endl;
}


int main(void) {
	process(900000000);
	return 0;
}

Results:

 Sample Size: 900000000
     Testing: Eratosthenes
Time elapsed: 55.4
       Found: 46009215

     Testing: Atkin
Time elapsed: 70.57
       Found: 46229727

Try the C++ code in the C++ compiler and see how it runs.

[Tellus]

Hm.. I am not that good at C++ syntax - classes and stuff (Willing to take up C# tho) and that's the reason I use C style. (Also Dev-C++ is supposed to be C/C++ compiler). Anyway.. I finally understood when the difference in speed is coming from.. When I asked for explanation on i*i<=size it was because in my version:

void EratosthenesSieve (unsigned long n)
{ clock_t start = clock();
  unsigned long j, i=2, br=0;
  while (i <= n ){
       if(sieveE[i] == 0){
            br++;
            for(j=i+i; j<= n ; j+=i)
              sieveE[j] = 1;
              
            
       }
       i++;
  }
printf("Eratosthenes counted: %10lu", br);
printf("\nTime elapsed: %f\n", ((double)clock() - start) / CLOCKS_PER_SEC);
}

I am counting the primes in the process of finding and I need to go to i<=n. (I got that just before I went to bed but was too lazy to edit (: ). As it seems I was willing too much Atkin to win over Era that I was thinking of ways to optimize Atkin only. So if I rewrite Era this way:

void EratosthenesSieve (unsigned long n)
{ clock_t start = clock();
  unsigned long j, i, br=0;
  sieveE[1]=1;
  for(i=2; i*i<=n ; i++)
       if(!sieveE[i]){
            for(j=i+i; j<= n ; j+=i)
              {sieveE[j] = 1;}
       }
                
  for(i=1; i<=n ; i+=2) //without counting evens
  if(!sieveE[i])br++;
printf("Eratosthenes counted: %10lu", br+1);//counting the only even prime (+1)
printf("\nTime elapsed: %f\n", ((double)clock() - start) / CLOCKS_PER_SEC);
}

There are significantly less checks.

#include<stdio.h>
#define MAXN  900000000
#include <stdlib.h>
#include <time.h>
#include <math.h>
const unsigned long n=MAXN;
char sieveE[MAXN];
char sieveA[MAXN];

void EratosthenesSieve (unsigned long n);
void AtkinSieve(unsigned long Limit);

int main(void) {
   
   
    AtkinSieve( n ) ;
    printf("\n");
    EratosthenesSieve( n );
    system("pause");    
    return 0;
}


void EratosthenesSieve (unsigned long n)
{ clock_t start = clock();
  unsigned long j, i, br=0;
  sieveE[1]=1;
  for(i=2; i*i<=n ; i++)
       if(!sieveE[i]){
            for(j=i+i; j<= n ; j+=i)
              {sieveE[j] = 1;}
       }
                
  for(i=1; i<=n ; i+=2)
  if(!sieveE[i])br++;
printf("Eratosthenes counted: %10lu", br+1);
printf("\nTime elapsed: %f\n", ((double)clock() - start) / CLOCKS_PER_SEC);
}

void AtkinSieve(unsigned long Limit)
{           clock_t start = clock();
            double sqrtLimit=sqrt(Limit)+1;
            unsigned long xsqd = 1;
            unsigned long xadd = 3;
            for (unsigned long x = 1; x<=sqrtLimit; x++)

            {       
                unsigned long ysqd = 1;
                unsigned long yadd = 3;
                for (unsigned long y = 1; y<=sqrtLimit; y++)
                {
                        

                        unsigned long n = (xsqd << 2) + ysqd; // n = 4*x*x + y*y

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

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

                        if (x > y)
                        {
                                n -= (ysqd << 1);  // n = 3*x*x - y*y
                                
                                if (n <= Limit && n % 12 == 11)
                                {
                                        sieveA[n] ^= true;
                                }
                        }
                        ysqd+=yadd;
                        yadd+=2;
                }
                 xsqd+=xadd;
                 xadd+=2;
        }                 
        unsigned long br=0;             
        sieveA[2]=1;sieveA[3]=1;
        for(unsigned long i=5 ; i<=sqrtLimit ; i++){
              if(sieveA[i]){
                           unsigned long j=i*i;
                           for (unsigned long k = j; k <= MAXN; k += j)
                                {
                                sieveA[k] = false;
                                }
                           }
        }           
        for(unsigned long i=1 ; i<=MAXN ; i+=2)// without checking for evens
        if(sieveA[i])br++;
        printf("Atkin counted: %10lu", br+1);
        printf("\nTime elapsed: %f\n", ((double)clock() - start) / CLOCKS_PER_SEC);
        
}

The result speak for itself:

Up to 900000000
Atkin counted: 46229727
Time elapsed: 37.361000

Eratosthenes counted: 46009215
Time elapsed: 30.799000

Also executed your code:

#include <iostream>
#include <ctime>
#include <cmath>

typedef unsigned long ulong;

struct PrimeFinder {
        bool *sieve;
        const int size;
        double elapsedTime;
        ulong found;
        
        PrimeFinder(int n) : size(n) { sieve = new bool[size+1]; }
        ~PrimeFinder() { delete [] sieve; }
        
        virtual const char *getName() const = 0;
        
        void findPrimes() {
                clock_t start;
                
                fill();
                start = clock();
                generate();
                elapsedTime = ((double)clock() - start) / CLOCKS_PER_SEC;
                found = count();
        }
        
        void print(std::ostream &out) const {
                out << "     Testing: " << getName() << std::endl;
                out << "Time elapsed: " << elapsedTime << std::endl;
                out << "       Found: " << found << std::endl;
                out << std::endl;
        }
        
protected:      
        virtual bool getInitValue() const = 0;
        virtual void generate() = 0;

        void fill() {
                bool value = getInitValue();
                for(ulong i=0; i<=size; i++) { sieve[i]=value; }
        }

        ulong count() const {
                ulong count = 1;
                bool value = sieve[2];
                for(ulong i=3; i<=size; i+=2) {
                        if(sieve[i]==value) { count++; }
                }
                return count;
        }
};


struct Atkin : public PrimeFinder {
        const char *getName() const { return "Atkin"; }
        Atkin(int n) : PrimeFinder(n) {}
protected:
        bool getInitValue() const { return false; }
        void generate() {
                bool *sieveA = this->sieve;
                ulong Limit = this->size;
                double sqrtLimit=sqrt(Limit)+1;
                ulong x, xsqd = 1, xadd = 3;
                for (x = 1; x<=sqrtLimit; x++) {       
                         ulong y, ysqd = 1, yadd = 3;
                         for (y = 1; y<=sqrtLimit; y++) {
                                ulong n = (xsqd << 2) + ysqd; // n = 4*x*x + y*y
                                if (n <= Limit && (n % 12 == 1 || n % 12 == 5)) { sieveA[n] ^= true; }
                                n -= xsqd;  // n = 3*x*x + y*y
                                if (n <= Limit && n % 12 == 7) { sieveA[n] ^= true; }
                                if (x > y) { 
                                        n -= (ysqd << 1);  // n = 3*x*x - y*y
                                        if (n <= Limit && n % 12 == 11){ sieveA[n] ^= true; }
                                }
                                ysqd+=yadd;
                                yadd+=2;
                        }
                        xsqd+=xadd;
                        xadd+=2;
                }    
                
                sieveA[2]=1; sieveA[3]=1;
                for(ulong i=5 ; i<=sqrtLimit ; i++){
                        if(sieveA[i]){
                                ulong k, j=i*i;
                                for (ulong k = j; k < Limit; k += j) { sieveA[k] = false; }
                        }
                }
        }
};


struct Eratosthenes : public PrimeFinder {
        const char *getName() const { return "Eratosthenes"; }
        Eratosthenes(int n) : PrimeFinder(n) {}
protected:
        bool getInitValue() const { return true; }
        void generate() {
                ulong limit = sqrt(size)+1;
                for (ulong i = 2; i <= limit; i++) {
                        if (sieve[i]) { 
                                for (ulong j = i+i; j <=size; j += i) { sieve[j] = false; } 
                        } 
                } 
        }
};


void run(PrimeFinder &pf) {
        pf.findPrimes();
        pf.print(std::cout);
}

void process(ulong sampleSize) {
        Eratosthenes s1(sampleSize);
        Atkin s2(sampleSize);

        std::cout << " Sample Size: " << sampleSize << std::endl;
        
        run(s1);
        run(s2);
        
        std::cout << std::endl;
}


int main(void) {
        process(900000000);
        system("pause");
        return 0;
}

The results from it:

Sample Size: 900000000
Testing: Eratosthenes
Time elapsed: 31.721
Found: 46009215

Testing: Atkin
Time elapsed: 36.816
Found: 46229727


Thanks.. mystery solved. It seems that Era is not only shorter (sourcewise), much easier to understand and remember, less memory consuming but also faster than its 'newer optimized version' - Atkin (which is supposed to win a competition for speed from mathematical point of view but oh well it is coding after all and the fact is a fact)

This post has been edited by Tellus: 19 April 2010 - 08:04 AM

[Tellus]

On the second thought I optimized Atkin a bit more so it is efficient now even for larger numbers (It is still slower than Era but it is getting closer):

#include<stdio.h>
#define MAXN  990009999
#include <stdlib.h>
#include <time.h>
#include <math.h>
const unsigned long n=MAXN;
char sieveE[MAXN];
char sieveA[MAXN];

void EratosthenesSieve (unsigned long n);
void AtkinSieve(unsigned long Limit);

int main(void) {
   
   
    AtkinSieve( n ) ;
    printf("\n");
    EratosthenesSieve( n );
    system("pause");    
    return 0;
}


void EratosthenesSieve (unsigned long n)
{ clock_t start = clock();
  unsigned long j, i, br=0;
  sieveE[1]=1;
  for(i=2; i*i<=n ; i++)
       if(!sieveE[i]){
            for(j=i+i; j<= n ; j+=i)
              {sieveE[j] = 1;}
       }
                
  for(i=1; i<=n ; i+=2)
  if(!sieveE[i])br++;
printf("\nEratosthenes counted: %10lu", br+1);
printf("\nTime elapsed: %f\n", ((double)clock() - start) / CLOCKS_PER_SEC);
}

void AtkinSieve(unsigned long Limit)
{           clock_t start = clock();
            unsigned long sqrtLimit=sqrt(Limit);
            unsigned long xsqd = 1;
            unsigned long xadd = 3;
            unsigned long LimitFour = Limit/4;            
            unsigned long LimitThree = Limit/3;
            unsigned long n;
            for (unsigned long x = 1; x<=sqrtLimit; x++)

            {       
                unsigned long ysqd = 1;
                unsigned long yadd = 3;
               
               
                for (unsigned long y = 1; y<=sqrtLimit; y++)
                {   
                              
                        
                         if(xsqd <= LimitFour-ysqd/4) //Since we check for Limit >= 4*x*x + y*y there is no point in
                                                      //going further than Sqrt(Limit/4)-ysqd/4 for the xsqd
                        {
                                             n = (xsqd << 2) + ysqd;  // n = 4*x*x + y*y
                                                                                         
                                              if (n % 12 == 5 || n % 12 == 1)
                                                {
                                                sieveA[n] ^= true;
                                                }
                                             
                        
                        }
                        
                        if(xsqd <= LimitThree-ysqd/3) //Since we check for Limit >= 3*x*x + y*y there is no point in 
                                                      //going further than Sqrt(Limit/3)-ysqd/3 for the xsqd
                        {
                                             n = (xsqd << 2) - xsqd + ysqd;  // n = 3*x*x + y*y
                                             
                                             
                                              if (n % 12 == 7)
                                                {
                                                sieveA[n] ^= true;
                                                }
                                          
                        }
                        
                        if ( (x > y) && (xsqd>ysqd/3) && (xsqd <= LimitThree+ysqd/3) ) 
                        //Here we have the system Limit >= 3*x*x-y*y >0
                        {
                                n = (xsqd << 2) - xsqd - ysqd;  // n = 3*x*x - y*y
                                
                                
                                 if (n % 12 == 11)
                                   {
                                        sieveA[n] ^= true;
                                   }
                              
                        }
                        ysqd+=yadd;
                        yadd+=2;
                }
                 xsqd+=xadd;
                 xadd+=2;
        }                 
        unsigned long br=0;             
        sieveA[2]=1;sieveA[3]=1;
        for(unsigned long i=5 ; i<=sqrtLimit ; i++){
              if(sieveA[i]){
                           unsigned long j=i*i;
                           for (unsigned long k = j; k <= MAXN; k += j)
                                {
                                sieveA[k] = false;
                                }
                           }
        }           
        for(unsigned long i=1 ; i<=MAXN ; i+=2) // without checking for evens
        if(sieveA[i])br++;
        printf("\nAtkin counted: %10lu", br+1);
        printf("\nTime elapsed: %f\n", ((double)clock() - start) / CLOCKS_PER_SEC);
        
}

Results:

Up to 900000000
Atkin counted:   46009215
Time elapsed: 32.491000


Eratosthenes counted:   46009215
Time elapsed: 30.896000
Press any key to continue . . .

Up to 990009999
Atkin counted:   50365183
Time elapsed: 35.143000


Eratosthenes counted:   50365183
Time elapsed: 34.188000

At least Atkin is counting proper this time and it is getting closer to Era as the numbers go larger.

EDIT: Tweaked Atkin a bit more:

#include<stdio.h>
#define MAXN  900000000
#include <stdlib.h>
#include <time.h>
#include <math.h>
const unsigned long n=MAXN;
char sieveE[MAXN];
char sieveA[MAXN];

void EratosthenesSieve (unsigned long n);
void AtkinSieve(unsigned long Limit);

int main(void) {
   
   
    AtkinSieve( n ) ;
    printf("\n");
    EratosthenesSieve( n );
    system("pause");    
    return 0;
}


void EratosthenesSieve (unsigned long n)
{ clock_t start = clock();
  unsigned long j, i, br=0;
  sieveE[1]=1;
  for(i=2; i*i<=n ; i++)
       if(!sieveE[i]){
            for(j=i+i; j<= n ; j+=i)
              {sieveE[j] = 1;}
       }
                
  for(i=1; i<=n ; i+=2)
  if(!sieveE[i])br++;
printf("\nEratosthenes counted: %10lu", br+1);
printf("\nTime elapsed: %f\n", ((double)clock() - start) / CLOCKS_PER_SEC);
}

void AtkinSieve(unsigned long Limit)
{           clock_t start = clock();
            unsigned long sqrtLimit=sqrt(Limit);
            unsigned long xsqd = 1;
            unsigned long xadd = 3;
            unsigned long LimitFour = Limit/4;            
            unsigned long LimitThree = Limit/3;
            unsigned long n;
            unsigned long fl;
            unsigned long fll=0;
            unsigned long ysqd = 1;
            unsigned long yadd = 3; 
            for (unsigned long x = 1; x<=sqrtLimit; x++)

            {       
                
                ysqd=1;
                yadd =3;
                fl=0;
                fll=0;
                for (unsigned long y = 1; y<=sqrtLimit; y++)
                {     if(fl==2 && fll==1)break; //no point checking further
                      else if(fl==1)goto one;
                      else if(fl==2)goto two;       
                        
                         if(xsqd <= LimitFour-ysqd/4-1) //Since we check for Limit >= 4*x*x + y*y there is no point in
                                                        //going further than Sqrt(Limit/4)-ysqd/4 for the xsqd
                        {
                                             n = (xsqd << 2) + ysqd;  // n = 4*x*x + y*y
                                                                                         
                                              if (n % 12 == 5 || n % 12 == 1)
                                                {
                                                sieveA[n] ^= true;
                                                }
                                             
                        
                        }
                        else fl++; // once this is exceeded no point checking again (goto one:)
                         
                         one:
                         if(xsqd <= LimitThree-ysqd/3) //Since we check for Limit >= 3*x*x + y*y there is no point in 
                                                      //going further than Sqrt(Limit/3)-ysqd/3 for the xsqd
                        {
                                             n = (xsqd << 2) - xsqd + ysqd;  // n = 3*x*x + y*y
                                             
                                             
                                              if (n % 12 == 7)
                                                {
                                                sieveA[n] ^= true;
                                                }
                                          
                        }
                        else fl++; // once this is exceeded no point checking again (goto two:)
                         
                         
                        
                        two:
                        
                        if(xsqd<=ysqd/3){fll=1;goto three;}//once n goes lower than 0 no point checking here (goto three:)
                        if ( (x > y) && (xsqd>ysqd/3) && (xsqd <= LimitThree+ysqd/3) ) 
                        //Here we have the system Limit >= 3*x*x-y*y >0
                        {
                                n = (xsqd << 2) - xsqd - ysqd;  // n = 3*x*x - y*y
                                
                                
                                 if (n % 12 == 11)
                                   {
                                        sieveA[n] ^= true;
                                   }
                              
                        }
                        
                        
                        three:
                        ysqd+=yadd;
                        yadd+=2;
                        
                        
                        
                }
                 xsqd+=xadd;
                 xadd+=2;
        }                 
        unsigned long br=0;             
        sieveA[2]=1;sieveA[3]=1;
        for(unsigned long i=5 ; i<=sqrtLimit ; i++){
              if(sieveA[i]){
                           unsigned long j=i*i;
                           for (unsigned long k = j; k <= MAXN; k += j)
                                {
                                sieveA[k] = false;
                                }
                           }
        }           
        for(unsigned long i=1 ; i<=MAXN ; i+=2) // without checking for evens
        if(sieveA[i])br++;
        printf("\nAtkin counted: %10lu", br+1);
        printf("\nTime elapsed: %f\n", ((double)clock() - start) / CLOCKS_PER_SEC);
        
}

Up to 900000000

Atkin counted:   46009215
Time elapsed: 31.512000


Eratosthenes counted:   46009215
Time elapsed: 30.695000

EDIT2: (lol)
It seems that goto is slow and it can be done without it:

#include<stdio.h>
#define MAXN  990009999
#include <stdlib.h>
#include <time.h>
#include <math.h>
const unsigned long n=MAXN;
char sieveE[MAXN];
char sieveA[MAXN];

void EratosthenesSieve (unsigned long n);
void AtkinSieve(unsigned long Limit);

int main(void) {
   
   
    AtkinSieve( n ) ;
    printf("\n");
    EratosthenesSieve( n );
    system("pause");    
    return 0;
}


void EratosthenesSieve (unsigned long n)
{ clock_t start = clock();
  unsigned long j, i, br=0;
  sieveE[1]=1;
  for(i=2; i*i<=n ; i++)
       if(!sieveE[i]){
            for(j=i+i; j<= n ; j+=i)
              {sieveE[j] = 1;}
       }
                
  for(i=1; i<=n ; i+=2)
  if(!sieveE[i])br++;
printf("\nEratosthenes counted: %10lu", br+1);
printf("\nTime elapsed: %f\n", ((double)clock() - start) / CLOCKS_PER_SEC);
}

void AtkinSieve(unsigned long Limit)
{           clock_t start = clock();
            unsigned long sqrtLimit=sqrt(Limit);
            unsigned long xsqd = 1;
            unsigned long xadd = 3;
            unsigned long LimitFour = Limit/4;            
            unsigned long LimitThree = Limit/3;
            unsigned long n;
            unsigned long fl; //flag
            unsigned long fll=0; //flag
            unsigned long ysqd = 1;
            unsigned long yadd = 3; 
            for (unsigned long x = 1; x<=sqrtLimit; x++)

            {       
                
                ysqd=1;
                yadd =3;
                fl=0;
                fll=0;
                for (unsigned long y = 1; y<=sqrtLimit; y++)
                {     if(fl==2 && fll==1)break; //no point checking further
                     
                            
                         if(fl<1)
                         {
                                 if(xsqd <= LimitFour-ysqd/4) //Since we check for Limit >= 4*x*x + y*y there is no point in
                                                        //going further than Sqrt(Limit/4)-ysqd/4 for the xsqd
                                 {
                                              n = (xsqd << 2) + ysqd;  // n = 4*x*x + y*y
                                                                                         
                                              if (n % 12 == 5 || n % 12 == 1)
                                                {
                                                sieveA[n] ^= true;
                                                }
                                             
                        
                                }
                                else fl++; // once this is exceeded no point checking again 
                        } 
                         
                         if(fl<2)
                         {
                                 if(xsqd <= LimitThree-ysqd/3) //Since we check for Limit >= 3*x*x + y*y there is no point in 
                                                      //going further than Sqrt(Limit/3)-ysqd/3 for the xsqd
                                 {
                                             n = (xsqd << 2) - xsqd + ysqd;  // n = 3*x*x + y*y
                                             
                                             
                                             if (n % 12 == 7)
                                                {
                                                sieveA[n] ^= true;
                                                }
                                          
                                 }
                                 else fl++; // once this is exceeded no point checking again 
                        } 
                         
                        
                        
                        if(!fll)
                        {
                                if(xsqd*3>ysqd){
                                //once n goes lower than 0 no point checking here 
                                       if ( (x > y) && (xsqd <= LimitThree+ysqd/3) ) 
                                       //Here we have the system Limit >= 3*x*x-y*y >0
                                              {
                                              n = (xsqd << 2) - xsqd - ysqd;  // n = 3*x*x - y*y
                                
                                
                                              if (n % 12 == 11)
                                                 {
                                                 sieveA[n] ^= true;
                                                 }
                              
                                              }
                              }
                              else fll=1;
                        }
                        
                       
                        ysqd+=yadd;
                        yadd+=2;
                        
                        
                        
                }
                 xsqd+=xadd;
                 xadd+=2;
        }                 
        unsigned long br=0;             
        sieveA[2]=1;sieveA[3]=1;
        for(unsigned long i=5 ; i<=sqrtLimit ; i++){
              if(sieveA[i]){
                           unsigned long j=i*i;
                           for (unsigned long k = j; k <= MAXN; k += j)
                                {
                                sieveA[k] = false;
                                }
                           }
        }           
        for(unsigned long i=1 ; i<=MAXN ; i+=2) // without checking for evens
        if(sieveA[i])br++;
        printf("\nAtkin counted: %10lu", br+1);
        printf("\nTime elapsed: %f\n", ((double)clock() - start) / CLOCKS_PER_SEC);
        
}

Final Results:

Up to 990009999
Atkin counted:   50365183
Time elapsed: 34.819000


Eratosthenes counted:   50365183
Time elapsed: 34.114000

Bloody close, eh? (:

This post has been edited by Tellus: 19 April 2010 - 09:14 AM

[Tellus]

Oh and I also increased the speed of Era (A LOT!):

#include<stdio.h>
#define MAXN  900000000
#include <stdlib.h>
#include <time.h>
#include <math.h>
const unsigned long n=MAXN, sqrtLimit=sqrt(n);
char sieveE[MAXN];
char sieveA[MAXN];

void EratosthenesSieve (unsigned long n);
void AtkinSieve(unsigned long Limit);

int main(void) {
   
   
    AtkinSieve( n ) ;
    printf("\n");
    EratosthenesSieve( n );
    system("pause");    
    return 0;
}


void EratosthenesSieve (unsigned long n)
{ clock_t start = clock();
  unsigned long j, i, br=0 , doubleI;
  sieveE[1]=1;

  for(i=3; i<=sqrtLimit ; i+=2)
       if(!sieveE[i]) 
            for( j=i*i; j<= n ; j+=i<<1 )
              sieveE[j] = 1;
       
                
  for(i=3; i<=n ; i+=2)
  if(!sieveE[i])br++;
printf("\nEratosthenes counted: %10lu", br+1);
printf("\nTime elapsed: %f\n", ((double)clock() - start) / CLOCKS_PER_SEC);
}


void AtkinSieve(unsigned long Limit)
{           clock_t start = clock();
            unsigned long sqrtLimit=sqrt(Limit);
            unsigned long xsqd = 1;
            unsigned long xadd = 3;
            unsigned long LimitFour = Limit/4;            
            unsigned long LimitThree = Limit/3;
            unsigned long n;
            unsigned long fl; //flag
            unsigned long fll; //flag
            unsigned long ysqd = 1;
            unsigned long yadd = 3; 
            for (unsigned long x = 1; x<=sqrtLimit; x++)

            {       
                
                ysqd=1;
                yadd =3;
                fl=0;
                fll=0;
                for (unsigned long y = 1; y<=sqrtLimit; y++)
                {     if(fl==2 && fll==1)break; //no point checking further
                     
                            
                         if(fl<1)
                         {
                                 if(xsqd <= LimitFour-ysqd/4) //Since we check for Limit >= 4*x*x + y*y there is no point in
                                                        //going further than Sqrt(Limit/4)-ysqd/4 for the xsqd
                                 {
                                              n = (xsqd << 2) + ysqd;  // n = 4*x*x + y*y
                                                                                         
                                              if (n % 12 == 5 || n % 12 == 1)
                                                {
                                                sieveA[n] ^= true;
                                                }
                                             
                        
                                }
                                else fl++; // once this is exceeded no point checking again 
                        } 
                         
                         if(fl<2)
                         {
                                 if(xsqd <= LimitThree-ysqd/3) //Since we check for Limit >= 3*x*x + y*y there is no point in 
                                                      //going further than Sqrt(Limit/3)-ysqd/4 for the xsqd
                                 {
                                             n = (xsqd << 2) - xsqd + ysqd;  // n = 3*x*x + y*y
                                             
                                             
                                             if (n % 12 == 7)
                                                {
                                                sieveA[n] ^= true;
                                                }
                                          
                                 }
                                 else fl++; // once this is exceeded no point checking again 
                        } 
                         
                        
                        
                        if(!fll)
                        {
                                if(xsqd>ysqd/3){
                                //once n goes lower than 0 no point checking here 
                                       if ( (x > y) && (xsqd <= LimitThree+ysqd/3) ) 
                                       //Here we have the system Limit >= 3*x*x-y*y >0
                                              {
                                              n = (xsqd << 2) - xsqd - ysqd;  // n = 3*x*x - y*y
                                
                                
                                              if (n % 12 == 11)
                                                 {
                                                 sieveA[n] ^= true;
                                                 }
                              
                                              }
                              }
                              else fll=1;
                        }
                        
                       
                        ysqd+=yadd;
                        yadd+=2;
                        
                        
                        
                }
                 xsqd+=xadd;
                 xadd+=2;
        }                 
        unsigned long br=0;             
        sieveA[2]=1;sieveA[3]=1;
        for(unsigned long i=5 ; i<=sqrtLimit ; i++){
              if(sieveA[i]){
                           unsigned long j=i*i;
                           for ( unsigned long k = j ; k <= MAXN ; k += j*2 )
                                {
                                sieveA[k] = false;
                                }
                           }
        }           
        for( unsigned long i=3 ; i<=MAXN ; i+=2 ) // without checking for evens
        if(sieveA[i])br++;
        printf("\nAtkin counted: %10lu", br+1);
        printf("\nTime elapsed: %f\n", ((double)clock() - start) / CLOCKS_PER_SEC);
}

Results:

Up to 979999999:
Atkin counted:   49881580
Time elapsed: 34.570000


Eratosthenes counted:   49881580
Time elapsed: 17.400000

That's it.. no more primes for me (will move on to graph trees at last!). Tried to squeeze the max of both algorithms.

This post has been edited by Tellus: 21 April 2010 - 07:58 AM