GopherProxy

	primes.c - 9base - revived minimalist port of Plan 9 userland to Unix
	git clone git://git.suckless.org/9base
	Log
	Files
	Refs
	README
	LICENSE
	---
	primes.c (2126B)
	---
	1 #include <u.h>
	2 #include <libc.h>
	3
	4 #define ptsiz (sizeof(pt)/sizeof(pt[0]))
	5 #define whsiz (sizeof(wheel)/sizeof(wheel[0]))
	6 #define tabsiz (sizeof(table)/sizeof(table[0]))
	7 #define tsiz8 (tabsiz*8)
	8
	9 double big = 9.007199254740992e15;
	10
	11 int pt[] =
	12 {
	13 2, 3, 5, 7, 11, 13, 17, 19, 23, 29,
	14 31, 37, 41, 43, 47, 53, 59, 61, 67, 71,
	15 73, 79, 83, 89, 97,101,103,107,109,113,
	16 127,131,137,139,149,151,157,163,167,173,
	17 179,181,191,193,197,199,211,223,227,229,
	18 };
	19 double wheel[] =
	20 {
	21 10, 2, 4, 2, 4, 6, 2, 6, 4, 2,
	22 4, 6, 6, 2, 6, 4, 2, 6, 4, 6,
	23 8, 4, 2, 4, 2, 4, 8, 6, 4, 6,
	24 2, 4, 6, 2, 6, 6, 4, 2, 4, 6,
	25 2, 6, 4, 2, 4, 2,10, 2,
	26 };
	27 uchar table[1000];
	28 uchar bittab[] =
	29 {
	30 1, 2, 4, 8, 16, 32, 64, 128,
	31 };
	32
	33 void mark(double nn, long k);
	34 void ouch(void);
	35
	36 void
	37 main(int argc, char *argp[])
	38 {
	39 int i;
	40 double k, temp, v, limit, nn;
	41
	42 if(argc <= 1) {
	43 fprint(2, "usage: primes starting [ending]\n");
	44 exits("usage");
	45 }
	46 nn = atof(argp[1]);
	47 limit = big;
	48 if(argc > 2) {
	49 limit = atof(argp[2]);
	50 if(limit < nn)
	51 exits(0);
	52 if(limit > big)
	53 ouch();
	54 }
	55 if(nn < 0 \|\| nn > big)
	56 ouch();
	57 if(nn == 0)
	58 nn = 1;
	59
	60 if(nn < 230) {
	61 for(i=0; i<ptsiz; i++) {
	62 if(pt[i] < nn)
	63 continue;
	64 if(pt[i] > limit)
	65 exits(0);
	66 print("%d\n", pt[i]);
	67 if(limit >= big)
	68 exits(0);
	69 }
	70 nn = 230;
	71 }
	72
	73 modf(nn/2, &temp);
	74 nn = 2.*temp + 1;
	75 /*
	76 * clear the sieve table.
	77 */
	78 for(;;) {
	79 for(i=0; i<tabsiz; i++)
	80 table[i] = 0;
	81 /*
	82 * run the sieve.
	83 */
	84 v = sqrt(nn+tsiz8);
	85 mark(nn, 3);
	86 mark(nn, 5);
	87 mark(nn, 7);
	88 for(i=0,k=11; k<=v; k+=wheel[i]) {
	89 mark(nn, k);
	90 i++;
	91 if(i >= whsiz)
	92 i = 0;
	93 }
	94 /*
	95 * now get the primes from the table
	96 * and print them.
	97 */
	98 for(i=0; i<tsiz8; i+=2) {
	99 if(table[i>>3] & bittab[i&07])
	100 continue;
	101 temp = nn + i;
	102 if(temp > limit)
	103 exits(0);
	104 print("%.0f\n", temp);
	105 if(limit >= big)
	106 exits(0);
	107 }
	108 nn += tsiz8;
	109 }
	110 }
	111
	112 void
	113 mark(double nn, long k)
	114 {
	115 double t1;
	116 long j;
	117
	118 modf(nn/k, &t1);
	119 j = k*t1 - nn;
	120 if(j < 0)
	121 j += k;
	122 for(; j<tsiz8; j+=k)
	123 table[j>>3] \|= bittab[j&07];
	124 }
	125
	126 void
	127 ouch(void)
	128 {
	129 fprint(2, "limits exceeded\n");
	130 exits("limits");
	131 }