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