primes.c (2126B)
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 }