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