primes.c - 9base - revived minimalist port of Plan 9 userland to Unix
 (HTM) git clone git://git.suckless.org/9base
 (DIR) Log
 (DIR) Files
 (DIR) Refs
 (DIR) README
 (DIR) LICENSE
       ---
       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 }