/* 
 * Copyright (c) 2003 Gianni Tedesco.
 * Released under the terms of the GNU GPL v2
 *
 * A program to print all prime numbers up to N using the sieve of
 * Eratosthenes. Requires N/8 bytes of memory to calculate. This
 * implementation uses a bit-array making it less sensitive to cache
 * timings and can generate the first 10,000,000 primes in just over a
 * second on my G4 667MHz. About 10 times faster than a byte array.
 */
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#ifndef BITS_PER_LONG
#define BITS_PER_LONG (sizeof(unsigned long)*8)
#endif

/* return number of longs needed to store (max) bits */
static inline unsigned int
slen(unsigned int max)
{
	unsigned int q, r;

	q = max / BITS_PER_LONG;
	r = max % BITS_PER_LONG;

	if ( r )
		q++;

	return q;
}

/* Requires GCC dynamic arrays */
void sieve(unsigned int max)
{
	unsigned long comps[slen(max)];
	unsigned int i, j, iscomp;

	memset(comps, 0, sizeof(comps));

	/* One is prime */
	printf("1\n");

	/* Start from two, the smallest prime from which composites are made */
	for(i=2; (i*i) < max; i++) {
		iscomp = comps[i/BITS_PER_LONG] & (1<<(i%BITS_PER_LONG));
		if ( iscomp )
			continue;

		/* This number is prime */
		printf("%i\n", i);

		/* If this number is prime, then mark all multiples as
		 * not-prime because they are composite
		 */
		for(j=i*2; j<max; j+=i)
			comps[j/BITS_PER_LONG] |= 1<<(j%BITS_PER_LONG);
	}

	/* Print all remaining primes which cannot produce composites
	 * within the requested range
	 */
	for(; i<max; i++) {
		iscomp = comps[i/BITS_PER_LONG] & (1<<(i%BITS_PER_LONG));
		if ( iscomp )
			continue;
		printf("%i\n", i);
	}
}

int main(int argc, char **argv)
{
	unsigned long int val = 32;
	char *end;

	if ( argc > 1 ) {
		val = strtoul(argv[1], &end, 0);
		if ( end==NULL || end==argv[1] || *end!='\0' ) {
			fprintf(stderr, "Usage: %s [max]\n", argv[0]);
			return 1;
		}
	}

	sieve(val);
	return 0;
}