harvey/sys/src/cmd/primes.c

/*
 * This file is part of the UCB release of Plan 9. It is subject to the license
 * terms in the LICENSE file found in the top-level directory of this
 * distribution and at http://akaros.cs.berkeley.edu/files/Plan9License. No
 * part of the UCB release of Plan 9, including this file, may be copied,
 * modified, propagated, or distributed except according to the terms contained
 * in the LICENSE file.
 */

#include <u.h>
#include <libc.h>
#include <bio.h>

double big = 9.007199254740992e15;

int pt[] =
{
	  2,  3,  5,  7, 11, 13, 17, 19, 23, 29,
	 31, 37, 41, 43, 47, 53, 59, 61, 67, 71,
	 73, 79, 83, 89, 97,101,103,107,109,113,
	127,131,137,139,149,151,157,163,167,173,
	179,181,191,193,197,199,211,223,227,229,
};
double wheel[] =
{
	10, 2, 4, 2, 4, 6, 2, 6, 4, 2,
	 4, 6, 6, 2, 6, 4, 2, 6, 4, 6,
	 8, 4, 2, 4, 2, 4, 8, 6, 4, 6,
	 2, 4, 6, 2, 6, 6, 4, 2, 4, 6,
	 2, 6, 4, 2, 4, 2,10, 2,
};
uint8_t table[1000];
uint8_t bittab[] =
{
	1, 2, 4, 8, 16, 32, 64, 128,
};

enum {
	ptsiz	= nelem(pt),
	whsiz	= nelem(wheel),
	tabsiz	= nelem(table),
	tsiz8	= tabsiz*8,
};

void	mark(double nn, int32_t k);

void
usage(void)
{
	fprint(2, "usage: %s [start [finish]]\n", argv0);
	exits("limits");
}

void
ouch(void)
{
	fprint(2, "limits exceeded\n");
	exits("limits");
}

void
main(int argc, char *argv[])
{
	int i;
	double k, temp, v, limit, nn;
	char *l;
	Biobuf bin;

	ARGBEGIN{
	default:
		usage();
		break;
	}ARGEND;

	nn = 0;
	limit = big;
	switch (argc) {
	case 0:
		Binit(&bin, 0, OREAD);
		while ((l = Brdline(&bin, '\n')) != nil) {
			if (*l == '\n')
				continue;
			nn = atof(l);
			if(nn < 0)
				sysfatal("negative start");
			break;
		}
		Bterm(&bin);
		break;
	case 2:
		limit = atof(argv[1]);
		if(limit < nn)
			exits(0);
		if(limit > big)
			ouch();
		/* fallthrough */
	case 1:
		nn = atof(argv[0]);
		break;
	default:
		usage();
		break;
	}

	if(nn < 0 || nn > big)
		ouch();
	if(nn == 0)
		nn = 1;

	if(nn < 230) {
		for(i=0; i<ptsiz; i++) {
			if(pt[i] < nn)
				continue;
			if(pt[i] > limit)
				exits(0);
			print("%d\n", pt[i]);
//			if(limit >= big)
//				exits(0);
		}
		nn = 230;
	}

	modf(nn/2, &temp);
	nn = 2*temp + 1;
/*
 *	clear the sieve table.
 */
	for(;;) {
		for(i = 0; i < tabsiz; i++)
			table[i] = 0;
/*
 *	run the sieve.
 */
		v = sqrt(nn+tsiz8);
		mark(nn, 3);
		mark(nn, 5);
		mark(nn, 7);
		for(i = 0, k = 11; k <= v; k += wheel[i]) {
			mark(nn, k);
			i++;
			if(i >= whsiz)
				i = 0;
		}
/*
 *	now get the primes from the table
 *	and print them.
 */
		for(i = 0; i < tsiz8; i += 2) {
			if(table[i>>3] & bittab[i&07])
				continue;
			temp = nn + i;
			if(temp > limit)
				exits(0);
			print("%.0f\n", temp);
//			if(limit >= big)
//				exits(0);
		}
		nn += tsiz8;
	}
}

void
mark(double nn, int32_t k)
{
	double t1;
	int32_t j;

	modf(nn/k, &t1);
	j = k*t1 - nn;
	if(j < 0)
		j += k;
	for(; j < tsiz8; j += k)
		table[j>>3] |= bittab[j&07];
}