Blog Entry © Sunday, July 27, 2025, A Bit of Programming Nostalgia Prime Number Related Programs by James Williams, Jr.

/*
  Author:  Pate Williams c 1995

  The following program is a solution to problem 18.15
  in Pascalgorithms by Edwin D. Reilly and Francis D.
  Federighi page 627. The program uses Simpson's rule
  to calculate the number of primes less than or equal
  a given number.
*/

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

typedef double real;

static real f(real x)
{
	return(1.0 / log(x));
}

static real simpson(int n, real a, real b)
{
	int i;
	real evensum, h, oddsum, twoh, x;

	if (n % 2 == 1) n = n - 1;
	h = (b - a) / n;
	twoh = h + h;
	x = h + a;
	oddsum = 0.0;
	for (i = 1; i <= n / 2; i++)
	{
		oddsum += f(x);
		x = twoh + x;
	}
	x = twoh + a;
	evensum = 0.0;
	for (i = 1; i <= n / 2 - 1; i++)
	{
		evensum += f(x);
		x = twoh + x;
	}
	return(h / 3.0 * (f(a) + f(b) + 4.0 * oddsum + 2.0 * evensum));
}

int main(void)
{
	int i, n, Nmaximum = 0, Nminimum = 0, Nstep;

	printf("n = ");			scanf_s("%d", &n);
	printf("N minimum = "); scanf_s("%d", &Nminimum);
	printf("N maximum = "); scanf_s("%d", &Nmaximum);
	printf("N step = ");	scanf_s("%d", &Nstep);
	printf("\n");
	printf("----------------------------------------\n");
	printf("Min\t\tMax\t\tprimes\n");
	printf("----------------------------------------\n");
	for (i = Nminimum; i <= Nmaximum; i += Nstep)
	{
		printf("%8d\t%8d\t%8.0lf\n", Nminimum, i + Nstep,
			simpson(n, Nminimum, i + Nstep));
	}
	printf("----------------------------------------\n");
	return(0);
}

n = 1024
N minimum = 0
N maximum = 10000000
N step = 1000000

----------------------------------------
Min             Max             primes
----------------------------------------
       0         1000000           78551
       0         2000000          148923
       0         3000000          216788
       0         4000000          283122
       0         5000000          348361
       0         6000000          412754
       0         7000000          476461
       0         8000000          539590
       0         9000000          602224
       0        10000000          664424
       0        11000000          726239
----------------------------------------

D:\PrimeCounter\x64\Release\PrimeCounter.exe (process 51884) exited with code 0 (0x0).
Press any key to close this window . . .

/*
  Author:  Pate Williams c 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 <math.h>
#include <stdio.h>

#define _WORD_SIZE 32
#define _VECT_SIZE 524288
#define SET_MIN    0
#define SET_MAX    16777215

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

SET set;

static int get_bit_pos(int* long_ptr, int* 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, int inset)
{
	int bit, word;

	if (get_bit_pos(&word, &bit, element))
	{
		if (inset > 0)
			set[word] |= (01 << bit);
		else
			set[word] &= ~(01 << bit);
	}
}

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

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

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

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

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

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

	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)));
	printf("%8ld\t%8ld\t%8.0lf\n", n, k, x);
}

int main(void)
{
	long n = 100L;

	printf("----------------------------------------\n");
	printf("n\t\tprimes\t\ttheory\n");
	printf("----------------------------------------\n");
	do
	{
		primes((int)n);
		n = 10L * n;
	} while (n < (long)SET_MAX);
	printf("----------------------------------------\n");
	return(0);
}

----------------------------------------
n               primes          theory
----------------------------------------
     100              25              29
    1000             168             181
   10000            1229            1261
  100000            9592            9634
 1000000           78498           78396
10000000          664579          665060
----------------------------------------

D:\Sieve\x64\Release\Sieve.exe (process 60092) exited with code 0 (0x0).
Press any key to close this window . . .