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