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