Blob
1 #include <u.h>2 #include <libc.h>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)9 double big = 9.007199254740992e15;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 };33 void mark(double nn, long k);34 void ouch(void);36 void37 main(int argc, char *argp[])38 {39 int i;40 double k, temp, v, limit, nn;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;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 }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 table96 * 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 }112 void113 mark(double nn, long k)114 {115 double t1;116 long j;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 }126 void127 ouch(void)128 {129 fprint(2, "limits exceeded\n");130 exits("limits");131 }