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 }