/* calcprime - Calculates all the prime numbers up to a compiled-in limit using an Eratosthenes sieve

   By James Stanley

   No license - Use as you wish */

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

/*
//Use this to get the most out of the program
#define LIMIT 8589934591ULL
#define SQRTLIMIT 92682
#define LONG long long
#define FORMAT "%llu\n"
*/


//Use this to check speed, needs 63MB of RAM
//NOTE: Currently 27.0s
#define LIMIT 1000000000
#define SQRTLIMIT 31623
#define LONG long
#define FORMAT "%lu\n"


/*
//Use this for checking accuracy
#define LIMIT 16000000
#define SQRTLIMIT 4000
#define LONG long
#define FORMAT "%lu\n"
*/

#define MEM_SIZE (LIMIT >> 4)+1


#define IS_PRIME(x) (sieve[x >> 4] & biton[(x & 0xf) >> 1])
#define SET_NOTPRIME(x) sieve[x >> 4] &= bitoff[(x & 0xf) >> 1]

unsigned char biton[8]  = { 0x01, 0x02, 0x04, 0x08, 0x10, 0x20, 0x40, 0x80 };
unsigned char bitoff[8] = { 0xfe, 0xfd, 0xfb, 0xf7, 0xef, 0xdf, 0xbf, 0x7f };

/* The following is slower than the above. Wierd. */
/*
#define IS_PRIME(x) (sieve[x >> 4] & shr_biton[x & 0xf])
#define SET_NOTPRIME(x) sieve[x >> 4] &= shr_bitoff[x & 0xf]

unsigned char shr_biton[0x10] = {
	0x01, 0x01, 0x02, 0x02, 0x04, 0x04, 0x08, 0x08, 0x10, 0x10, 0x20, 0x20,
	0x40, 0x40, 0x80, 0x80
};

unsigned char shr_bitoff[0x10] = {
	0xfe, 0xfe, 0xfd, 0xfd, 0xfb, 0xfb, 0xf7, 0xf7, 0xef, 0xef, 0xdf, 0xdf,
	0xbf, 0xbf, 0x7f, 0x7f
};
*/

int main() {
  unsigned char *sieve;
  unsigned LONG i;
	unsigned LONG j;
  unsigned LONG i2;
  unsigned LONG is;
  //time_t start = time(0);
	//unsigned long long c = 0;
	//unsigned long long t = 0;
  
  //fprintf(stderr, "%d\n", (int)clock());
  
  sieve = malloc(MEM_SIZE);
	memset(sieve, 0xff, MEM_SIZE);
  
  //fprintf(stderr, "%d\n", (int)clock());
  
	for(i = 3; i < SQRTLIMIT; i += 2) {
    if(IS_PRIME(i)) {//Only mark off multiples of primes
      i2 = i * 2;
      is = i * i;//We start off marking at it's square
			/* WARNING: the commented out test must be included if
				   LIMIT + 2 * SQRTLIMIT
				 is larger than the largest number that can fit in j!!! */
      for(j = is; (j < LIMIT) /* && (j >= is) */; j += i2) {
				//c++;
				if(IS_PRIME(j)) {
					SET_NOTPRIME(j);
					//t++;
				}
      }
		}
  }

	//printf("%lld total, %lld new, %lld%% required\n", c, t, (t*100)/c);
  
  //fprintf(stderr, "%d\n", (int)clock());
  
  SET_NOTPRIME(0);
  SET_NOTPRIME(1);

  
/*   printf("2\n"); */
/*   for(i = 1; i < LIMIT; i += 2) { */
/*     if(IS_PRIME(i)) printf(FORMAT, i); */
/*   } */

/*   printf("%ds\nHit ENTER...", time(0) - start); */
/*   getchar(); */
  
  
  //fprintf(stderr, "%d\n", (int)clock());
  
  return 0;
}