Set Implementation of Sieve of Eratosthenes by James Pate Williams, Jr.

/*
Author:  Pate Williams (c) January 20, 1995
The following is a translation of the Pascal program
sieve found in Pascalgorithms by Edwin D. Reilly and
Francis D. Federighi page 652. This program uses sets
to represent the sieve (see C Programming Language An
Applied Perspective by Lawrence Miller and Alec Qui-
lici pages 160 - 162).
*/
#include "stdafx.h"
#include <math.h>
#include <stdio.h>
#include <time.h>

#define _WORD_SIZE 32
#define _VECT_SIZE 31250000
#define SET_MIN    0
#define SET_MAX    1000000000

typedef long LONG;
typedef long SET[_VECT_SIZE];
typedef LONG ELEMENT;

SET set;

static LONG get_bit_pos(LONG *long_ptr, LONG *bit_ptr,
	ELEMENT element)
{
	*long_ptr = element / _WORD_SIZE;
	*bit_ptr = element % _WORD_SIZE;
	return(element >= SET_MIN && element <= SET_MAX);
}

static void set_bit(ELEMENT element, LONG inset)
{
	LONG bit, word;

	if (get_bit_pos(&word, &bit, element))
		inset ? set[word] |= (01 << bit) :
		set[word] &= ~(01 << bit);
}

static int get_bit(ELEMENT element)
{
	LONG bit, word;

	return(get_bit_pos(&word, &bit, element) ?
		(set[word] >> bit) & 01 : 0);
}

void set_Add(ELEMENT element)
{
	set_bit(element, 1);
}

void set_Del(ELEMENT element)
{
	set_bit(element, 0);
}

int set_Mem(ELEMENT element)
{
	return get_bit(element);
}

void primes(LONG n)
{
	LONG c, i, inc, k;
	double x;

	clock_t now = clock();
	set_Add(2);
	for (i = 3; i <= n; i++)
		if ((i + 1) % 2 == 0)
			set_Add(i);
		else
			set_Del(i);
	c = 3;
	do
	{
		i = c * c;
		inc = c + c;
		while (i <= n)
		{
			set_Del(i);
			i = i + inc;
		}
		c += 2;
		while (set_Mem(c) == 0) c += 1;
	} while (c * c <= n);
	k = 0;
	for (i = 2; i <= n; i++)
		if (set_Mem(i) == 1) k++;
	x = n / log(n) - 5.0;
	x = x + exp(1.0 + 0.15 * log(n) * sqrt(log(n)));
	clock_t later = clock();
	double runtime = (double)(later - now) / CLOCKS_PER_SEC;
	printf("%10ld\t%10ld\t%10.0lf\t%6.4lf\n",
		n, k, x, runtime);
}

int main(void)
{
	LONG n = 10L;

	printf("--------------------------------------------------------\n");
	printf("n\t\tprimes\t\ttheory\t\ttime (s)\n");
	printf("--------------------------------------------------------\n");
	do
	{
		primes(n);
		clock_t later = clock();
		n = 10L * n;
	} while (n < 1000000000);
	printf("--------------------------------------------------------\n");
	return(0);
}