Blob
1 #include <u.h>2 #include <libc.h>3 #include <bio.h>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)10 double big = 9.007199254740992e15;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 };34 void mark(double nn, long k);35 void ouch(void);36 Biobuf bout;38 void39 main(int argc, char *argp[])40 {41 int i;42 double k, temp, v, limit, nn;44 Binit(&bout, 1, OWRITE);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;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 }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 table100 * 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 }116 void117 mark(double nn, long k)118 {119 double t1;120 long j;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 }130 void131 ouch(void)132 {133 fprint(2, "limits exceeded\n");134 exits("limits");135 }