primebnc -- prime generation benchmark

// primebnc/primebnc -- prime algorithm benchmark

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

//#define MAXCNT		10001
//#define MAXCNT		20002
//#define MAXCNT		30003
//#define MAXCNT		40004
#define MAXCNT		1000000

int opt_f;

typedef unsigned long long val_t;

int tstno;
int inited;
char *reason;

int pcntold;
double bestrate;
double tvprev;

int maxcnt;
val_t pxold[MAXCNT + 1];
val_t pxnow[MAXCNT + 1];

int
prime1(val_t *primes)
{
	val_t i;
	int idx;
	int isprime;
	val_t n = 2;

	idx = 0;
	while (idx <= maxcnt) {
		isprime = 1;
		for (i = 2; i < n; i++) {
			if (n % i == 0) {
				isprime = 0;
				break;
			}
		}

		if (isprime) {
			primes[idx] = n;
			idx++;
		}

		n++;
	}

	tstno = 1;
	reason = "baseline 2 to n";

	return idx;
}

int
prime2(val_t *primes)
{
	val_t i;
	int idx;
	int isprime;
	val_t xsqrt;
	val_t n = 2;
	val_t n2;

	xsqrt = 0;

	idx = 0;
	primes[idx++] = 2;

	while (idx <= maxcnt) {
		// get sqrt(n)
		while ((xsqrt * xsqrt) < n) {
			xsqrt += 1;
		}

		isprime = 1;
		for (i = 2; i <= xsqrt; i++) {
			if (n % i == 0) {
				isprime = 0;
				break;
			}
		}

		if (isprime) {
			primes[idx] = n;
			idx++;
		}

		n2 = n + 1;
		if (n2 < n)
			printf("overflow: idx=%d\n",idx);
		n = n2;
	}

	tstno = 2;
	reason = "2 to sqrt by 2";

	return idx;
}

int
prime3(val_t *primes)
{
	val_t i;
	int idx;
	int isprime;
	val_t xsqrt;
	val_t n;

	xsqrt = 0;

	idx = 0;
	primes[idx++] = 2;
	primes[idx++] = 3;
	n = 5;

	while (idx <= maxcnt) {
		// get sqrt(n)
		while ((xsqrt * xsqrt) < n) {
			xsqrt += 1;
		}

		isprime = 1;
		for (i = 3; i <= xsqrt; i += 2) {
			if (n % i == 0) {
				isprime = 0;
				break;
			}
		}

		if (isprime) {
			primes[idx] = n;
			idx++;
		}

		n += 2;
	}

	tstno = 3;
	reason = "3 to sqrt by 2";

	return idx;
}

int
prime4(val_t *primes)
{
	val_t i;
	int idx;
	int isprime;
	val_t xsqrt;
	val_t n;
	val_t lo;
	val_t hi;

	xsqrt = 0;

	idx = 0;
	primes[idx++] = 2;
	primes[idx++] = 3;
	n = 6;

	while (idx <= maxcnt) {
		lo = n - 1;
		hi = n + 1;

		// get sqrt(n)
		while ((xsqrt * xsqrt) < hi) {
			xsqrt += 1;
		}

		isprime = 3;
		for (i = 3;  i <= xsqrt;  i += 2) {
			if (isprime & 1) {
				if (lo % i == 0)
					isprime &= ~1;
			}
			if (isprime & 2) {
				if (hi % i == 0)
					isprime &= ~2;
			}
			if (! isprime)
				break;
		}

		if (isprime & 1) {
			primes[idx] = lo;
			idx++;
		}

		if (isprime & 2) {
			primes[idx] = hi;
			idx++;
		}

		n += 6;
	}

	tstno = 4;
	reason = "6 to sqrt by 6 (combined 6n-1/6n+1 loops)";

	return idx;
}

int
prime5(val_t *primes)
{
	val_t i;
	int idx;
	int isprime;
	val_t xsqrt;
	val_t n;
	val_t lo;
	val_t hi;

	xsqrt = 0;

	idx = 0;
	primes[idx++] = 2;
	primes[idx++] = 3;
	n = 6;

	while (idx <= maxcnt) {
		lo = n - 1;
		hi = n + 1;

		// get sqrt(n)
		while ((xsqrt * xsqrt) < hi) {
			xsqrt += 1;
		}

		isprime = 1;
		for (i = 3; i <= xsqrt; i += 2) {
			if (lo % i == 0) {
				isprime = 0;
				break;
			}
		}
		if (isprime) {
			primes[idx] = lo;
			idx++;
		}

		isprime = 1;
		for (i = 3; i <= xsqrt; i += 2) {
			if (hi % i == 0) {
				isprime = 0;
				break;
			}
		}
		if (isprime) {
			primes[idx] = hi;
			idx++;
		}

		n += 6;
	}

	tstno = 5;
	reason = "6 to sqrt by 6 (separate 6n-1/6n+1 loops)";

	return idx;
}

int
prime6(val_t *primes)
{
	int cnt;
	int isprime;
	val_t xsqrt;
	val_t n;
	val_t lo;
	val_t hi;
	val_t pval;
	val_t *pptr;
	val_t *pend;

	xsqrt = 0;

	cnt = 0;
	primes[cnt++] = 2;
	primes[cnt++] = 3;
	n = 6;

	while (cnt <= maxcnt) {
		lo = n - 1;
		hi = n + 1;

		// get sqrt(n)
		while ((xsqrt * xsqrt) < hi) {
			xsqrt += 1;
		}

		isprime = 3;
		pptr = primes;
		pend = primes + cnt;
		for (;  pptr < pend;  ++pptr) {
			pval = *pptr;

			// early stop if we exceed square root of number being tested
			if (pval > xsqrt)
				break;

			// test 6n - 1
			if (isprime & 1) {
				if ((lo % pval) == 0)
					isprime &= ~1;
			}

			// test 6n + 1
			if (isprime & 2) {
				if ((hi % pval) == 0)
					isprime &= ~2;
			}

			// bug out if both are non-prime
			if (! isprime)
				break;
		}

		// 6n - 1 is prime
		if (isprime & 1) {
			primes[cnt] = lo;
			cnt++;
		}

		// 6n + 1 is prime
		if (isprime & 2) {
			primes[cnt] = hi;
			cnt++;
		}

		n += 6;
	}

	tstno = 6;
	reason = "6 to sqrt by prime list (combined 6n-1/6n+1 loops)";

	return cnt;
}

int
prime7(val_t *primes)
{
	int cnt;
	int isprime;
	val_t xsqrt;
	val_t n;
	val_t lo;
	val_t hi;
	val_t pval;
	val_t *pptr;
	val_t *pend;

	xsqrt = 0;

	cnt = 0;
	primes[cnt++] = 2;
	primes[cnt++] = 3;
	n = 6;

	while (cnt <= maxcnt) {
		lo = n - 1;
		hi = n + 1;

		// get sqrt(n)
		while ((xsqrt * xsqrt) < hi) {
			xsqrt += 1;
		}

		// check for 6n - 1 is prime
		isprime = 1;
		pptr = primes;
		pend = primes + cnt;
		for (;  pptr < pend;  ++pptr) {
			pval = *pptr;

			// early stop if we exceed square root of number being tested
			if (pval > xsqrt)
				break;

			// test 6n - 1
			if ((lo % pval) == 0) {
				isprime = 0;
				break;
			}
		}
		if (isprime) {
			primes[cnt] = lo;
			cnt++;
		}

		// check for 6n + 1 is prime
		isprime = 1;
		pptr = primes;
		pend = primes + cnt;
		for (;  pptr < pend;  ++pptr) {
			pval = *pptr;

			// early stop if we exceed square root of number being tested
			if (pval > xsqrt)
				break;

			// test 6n + 1
			if ((hi % pval) == 0) {
				isprime = 0;
				break;
			}
		}
		if (isprime) {
			primes[cnt] = hi;
			cnt++;
		}

		n += 6;
	}

	tstno = 7;
	reason = "6 to sqrt by prime list (separate 6n-1/6n+1 loops)";

	return cnt;
}

double
tvgetf(void)
{
	struct timespec ts;
	double sec;

	clock_gettime(CLOCK_REALTIME,&ts);

	sec = ts.tv_nsec;
	sec /= 1e9;
	sec += ts.tv_sec;

	return sec;
}

void
timeit(int (*p)(val_t *))
{
	val_t *pnow;
	val_t *pold;
	int pcntact;
	double tvbeg;
	double tvend;
	double rate;
	double ratio;

	printf("---------------\n");

	pold = pxold;
	pnow = inited ? pxnow : pxold;

	// load up the cache
	for (int i = 0; i < maxcnt; i++)
		pnow[i] = 1;

	tvbeg = tvgetf();
	pcntact = p(pnow);
	tvend = tvgetf();

	tvend -= tvbeg;

	// show prime generation rate
	rate = (double) maxcnt / tvend;
	printf("prime%d: %.9f (%.3f primes/sec",tstno,tvend,rate);
	if (inited) {
		ratio = rate / bestrate;
		printf(" %.3fx %s",ratio,(ratio > 1.0) ? "faster" : "slower");
	}
	printf(")");

	do {
		if (! inited) {
			pcntold = pcntact;
			bestrate = rate;
			printf(" -- %s\n",reason);
			break;
		}

		// show time ratio
		ratio = tvprev / tvend;
		printf(" %.3fx (%s than previous)",
			ratio,(ratio > 1.0) ? "faster" : "slower");

		printf(" -- %s\n",reason);

		for (int i = 0; i < maxcnt; i++) {
			if (pnow[i] != pold[i]) {
				printf("%d: pold=%lld pnow=%lld\n",i,pold[i],pnow[i]);
				break;
			}
		}
	} while (0);

	tvprev = tvend;
	inited = 1;
}

int
main(int argc,char **argv)
{
	char *cp;

	--argc;
	++argv;

	maxcnt = MAXCNT;

	for (;  argc > 0;  --argc, ++argv) {
		cp = *argv;
		if (*cp != '-')
			break;

		switch (cp[1]) {
		case 'f':
			opt_f = 1;
			break;
		case 'N':
			maxcnt = strtol(cp + 2,&cp,10);
			break;
		}
	}

	if (opt_f)
		maxcnt = 40004;

	if (maxcnt > MAXCNT)
		maxcnt = MAXCNT;

	// this takes a whole minute
	if (opt_f)
		timeit(prime1);

	// these are _much_ faster
	timeit(prime2);
	timeit(prime3);
	timeit(prime4);
	timeit(prime5);
	timeit(prime6);
	timeit(prime7);

#if 0
	for (int i = 0; i < maxcnt; i++) {
		printf("%d\n", pxold[i]);
	}
#endif

	return 0;
}