sys/src/cmd/primes.c

   1 #include <u.h>
   2 #include <libc.h>
   3 #include <bio.h>
   4
   5 double big = 9.007199254740992e15;
   6
   7 int pt[] =
   8 {
   9           2,  3,  5,  7, 11, 13, 17, 19, 23, 29,
  10          31, 37, 41, 43, 47, 53, 59, 61, 67, 71,
  11          73, 79, 83, 89, 97,101,103,107,109,113,
  12         127,131,137,139,149,151,157,163,167,173,
  13         179,181,191,193,197,199,211,223,227,229,
  14 };
  15 double wheel[] =
  16 {
  17         10, 2, 4, 2, 4, 6, 2, 6, 4, 2,
  18          4, 6, 6, 2, 6, 4, 2, 6, 4, 6,
  19          8, 4, 2, 4, 2, 4, 8, 6, 4, 6,
  20          2, 4, 6, 2, 6, 6, 4, 2, 4, 6,
  21          2, 6, 4, 2, 4, 2,10, 2,
  22 };
  23 uchar table[1000];
  24 uchar bittab[] =
  25 {
  26         1, 2, 4, 8, 16, 32, 64, 128,
  27 };
  28
  29 enum {
  30         ptsiz   = nelem(pt),
  31         whsiz   = nelem(wheel),
  32         tabsiz  = nelem(table),
  33         tsiz8   = tabsiz*8,
  34 };
  35
  36 void    mark(double nn, long k);
  37
  38 void
  39 usage(void)
  40 {
  41         fprint(2, "usage: %s [start [finish]]\n", argv0);
  42         exits("limits");
  43 }
  44
  45 void
  46 ouch(void)
  47 {
  48         fprint(2, "limits exceeded\n");
  49         exits("limits");
  50 }
  51
  52 void
  53 main(int argc, char *argv[])
  54 {
  55         int i;
  56         double k, temp, v, limit, nn;
  57         char *l;
  58         Biobuf bin;
  59
  60         ARGBEGIN{
  61         default:
  62                 usage();
  63                 break;
  64         }ARGEND;
  65
  66         nn = 0;
  67         limit = big;
  68         switch (argc) {
  69         case 0:
  70                 Binit(&bin, 0, OREAD);
  71                 while ((l = Brdline(&bin, '\n')) != nil) {
  72                         if (*l == '\n')
  73                                 continue;
  74                         nn = atof(l);
  75                         if(nn < 0)
  76                                 sysfatal("negative start");
  77                         break;
  78                 }
  79                 Bterm(&bin);
  80                 break;
  81         case 2:
  82                 limit = atof(argv[1]);
  83                 if(limit < nn)
  84                         exits(0);
  85                 if(limit > big)
  86                         ouch();
  87                 /* fallthrough */
  88         case 1:
  89                 nn = atof(argv[0]);
  90                 break;
  91         default:
  92                 usage();
  93                 break;
  94         }
  95
  96         if(nn < 0 || nn > big)
  97                 ouch();
  98         if(nn == 0)
  99                 nn = 1;
 100
 101         if(nn < 230) {
 102                 for(i=0; i<ptsiz; i++) {
 103                         if(pt[i] < nn)
 104                                 continue;
 105                         if(pt[i] > limit)
 106                                 exits(0);
 107                         print("%d\n", pt[i]);
 108 //                      if(limit >= big)
 109 //                              exits(0);
 110                 }
 111                 nn = 230;
 112         }
 113
 114         modf(nn/2, &temp);
 115         nn = 2*temp + 1;
 116 /*
 117  *      clear the sieve table.
 118  */
 119         for(;;) {
 120                 for(i = 0; i < tabsiz; i++)
 121                         table[i] = 0;
 122 /*
 123  *      run the sieve.
 124  */
 125                 v = sqrt(nn+tsiz8);
 126                 mark(nn, 3);
 127                 mark(nn, 5);
 128                 mark(nn, 7);
 129                 for(i = 0, k = 11; k <= v; k += wheel[i]) {
 130                         mark(nn, k);
 131                         i++;
 132                         if(i >= whsiz)
 133                                 i = 0;
 134                 }
 135 /*
 136  *      now get the primes from the table
 137  *      and print them.
 138  */
 139                 for(i = 0; i < tsiz8; i += 2) {
 140                         if(table[i>>3] & bittab[i&07])
 141                                 continue;
 142                         temp = nn + i;
 143                         if(temp > limit)
 144                                 exits(0);
 145                         print("%.0f\n", temp);
 146 //                      if(limit >= big)
 147 //                              exits(0);
 148                 }
 149                 nn += tsiz8;
 150         }
 151 }
 152
 153 void
 154 mark(double nn, long k)
 155 {
 156         double t1;
 157         long j;
 158
 159         modf(nn/k, &t1);
 160         j = k*t1 - nn;
 161         if(j < 0)
 162                 j += k;
 163         for(; j < tsiz8; j += k)
 164                 table[j>>3] |= bittab[j&07];
 165 }