tprimes.c - plan9port - [fork] Plan 9 from user space
HTML git clone git://src.adamsgaard.dk/plan9port
DIR Log
DIR Files
DIR Refs
DIR README
DIR LICENSE
---
tprimes.c (2202B)
---
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 }