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