primes.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 }