chiark - git - mdw - catacomb/blob - math/primeiter.c

   1 /* -*-c-*-
   2  *
   3  * Iterate over small primes efficiently
   4  *
   5  * (c) 2007 Straylight/Edgeware
   6  */
   7
   8 /*----- Licensing notice --------------------------------------------------*
   9  *
  10  * This file is part of Catacomb.
  11  *
  12  * Catacomb is free software; you can redistribute it and/or modify
  13  * it under the terms of the GNU Library General Public License as
  14  * published by the Free Software Foundation; either version 2 of the
  15  * License, or (at your option) any later version.
  16  *
  17  * Catacomb is distributed in the hope that it will be useful,
  18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  20  * GNU Library General Public License for more details.
  21  *
  22  * You should have received a copy of the GNU Library General Public
  23  * License along with Catacomb; if not, write to the Free
  24  * Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
  25  * MA 02111-1307, USA.
  26  */
  27
  28 /*----- Header files ------------------------------------------------------*/
  29
  30 #include "fibrand.h"
  31 #include "mp.h"
  32 #include "pgen.h"
  33 #include "primeiter.h"
  34 #include "primetab.h"
  35 #include "wheel.h"
  36
  37 /*----- Theory ------------------------------------------------------------*
  38  *
  39  * For small primes, we can just pluck them out of the small primes table.
  40  * For larger primes, we can test them individually, or build a sieve or
  41  * something, but since we don't know when to stop, that could be tricky.
  42  *
  43  * We've built a `wheel', as follows.  Let %$m$% be the product of the first
  44  * %$n$% primes.  There are %$\phi(m)$% integers %$n_i$%, with %$0 < n_i <
  45  * m$% coprime to %$m$%, and any integer %$j > n$% must be congruent to some
  46  * %$n_i$% modulo %$m$%.  The wheel itself doesn't list the %$n_i$%, but
  47  * rather the differences %$\delta_i = n_i - n_{i-1}$% (wrapping around
  48  * appropriately at the ends), so you can just add simple offsets to step
  49  * onwards.  The wheel assumes you start at 1 and move on round.
  50  */
  51
  52 /*----- Main code ---------------------------------------------------------*/
  53
  54 /* --- @wheelsync@ --- *
  55  *
  56  * Arguments:   @primeiter *pi@ = iterator to synchronize
  57  *              @mp *where@ = value to synchronize
  58  *
  59  * Returns:     ---
  60  *
  61  * Use:         Sets up the wheel index to match the given integer.  After
  62  *              this, we can step along the wheel to find candidate primes.
  63  */
  64
  65 static void wheelsync(primeiter *pi, mp *where)
  66 {
  67   mpw w;
  68   mp t;
  69   mpw rr;
  70   mp *r = MP_NEW;
  71   unsigned i, n;
  72
  73   w = WHEELMOD;
  74   mp_build(&t, &w, &w + 1);
  75   mp_div(0, &r, where, &t);
  76   rr = MP_ZEROP(r) ? 0 : r->v[0];
  77
  78   for (i = 0, n = 1; rr > n; n += wheel[i], i++);
  79   w = n - rr;
  80   pi->p = mp_add(MP_NEW, where, &t);
  81   pi->i = i;
  82   pi->r = fibrand_create(0);
  83   MP_DROP(r);
  84 }
  85
  86 /* --- @primeiter_create@ --- *
  87  *
  88  * Arguments:   @primeiter *pi@ = pointer to an iterator structure
  89  *              @mp *start@ = where to start
  90  *
  91  * Returns:     ---
  92  *
  93  * Use:         Initializes a prime iterator.  The first output will be the
  94  *              smallest prime not less than @start@.
  95  */
  96
  97 void primeiter_create(primeiter *pi, mp *start)
  98 {
  99   mpw n;
 100   unsigned l, h, m;
 101
 102   if (!start || MP_CMP(start, <=, MP_TWO))
 103     start = MP_TWO;
 104
 105   if (MP_LEN(start) <= 1) {
 106     n = start->v[0];
 107     if (n <= MAXPRIME) {
 108       l = 0;
 109       h = NPRIME;
 110       for (;;) {
 111         m = l + (h - l)/2;
 112         if (primetab[m] == n) break;
 113         else if (m == l) { m++; break; }
 114         else if (primetab[m] < n) l = m;
 115         else h = m;
 116       }
 117       pi->i = m;
 118       pi->mode = PIMODE_PTAB;
 119       mp_build(&pi->pp, &pi->w, &pi->w + 1);
 120       pi->p = &pi->pp;
 121       return;
 122     }
 123   }
 124
 125   wheelsync(pi, start);
 126   pi->mode = PIMODE_STALL;
 127 }
 128
 129 /* --- @primeiter_destroy@ --- *
 130  *
 131  * Arguments:   @primeiter *pi@ = pointer to iterator structure
 132  *
 133  * Returns:     ---
 134  *
 135  * Use:         Frees up an iterator structure when it's no longer wanted.
 136  */
 137
 138 void primeiter_destroy(primeiter *pi)
 139 {
 140   switch (pi->mode) {
 141     case PIMODE_PTAB:
 142       break;
 143     case PIMODE_STALL:
 144     case PIMODE_WHEEL:
 145       MP_DROP(pi->p);
 146       GR_DESTROY(pi->r);
 147       break;
 148     default:
 149       abort();
 150   }
 151 }
 152
 153 /* --- @primeiter_next@ --- *
 154  *
 155  * Arguments:   @primeiter *pi@ = pointer to an iterator structure
 156  *              @mp *d@ = fake destination
 157  *
 158  * Returns:     The next prime number from the iterator.
 159  *
 160  * Use:         Returns a new prime number.
 161  */
 162
 163 mp *primeiter_next(primeiter *pi, mp *d)
 164 {
 165   mp *p;
 166
 167   switch (pi->mode) {
 168     case PIMODE_PTAB:
 169       pi->w = primetab[pi->i++];
 170       if (pi->i >= NPRIME) {
 171         wheelsync(pi, pi->p);
 172         pi->mode = PIMODE_WHEEL;
 173       }
 174       p = MP_COPY(pi->p);
 175       MP_SPLIT(p);
 176       break;
 177     case PIMODE_STALL:
 178       pi->mode = PIMODE_WHEEL;
 179       goto loop;
 180     case PIMODE_WHEEL:
 181       do {
 182         MP_DEST(pi->p, MP_LEN(pi->p) + 1, pi->p->f);
 183         MPX_UADDN(pi->p->v, pi->p->vl, wheel[pi->i++]);
 184         MP_SHRINK(pi->p);
 185         if (pi->i >= WHEELN) pi->i = 0;
 186       loop:;
 187       } while (!pgen_primep(pi->p, pi->r));
 188       p = MP_COPY(pi->p);
 189       break;
 190     default:
 191       abort();
 192   }
 193   if (d) MP_DROP(d);
 194   return (p);
 195 }
 196
 197 /*----- Test rig ----------------------------------------------------------*/
 198
 199 #ifdef TEST_RIG
 200
 201 #include <mLib/macros.h>
 202 #include <mLib/testrig.h>
 203
 204 static int test(dstr *v)
 205 {
 206   mp *start = *(mp **)v[0].buf;
 207   mp *pp[5], *ret[5];
 208   int i;
 209   primeiter pi;
 210   int ok = 1;
 211
 212   for (i = 0; i < N(pp); i++)
 213     pp[i] = *(mp **)v[i + 1].buf;
 214   primeiter_create(&pi, start);
 215   for (i = 0; i < N(pp); i++) {
 216     ret[i] = primeiter_next(&pi, MP_NEW);
 217     if (!MP_EQ(ret[i], pp[i])) ok = 0;
 218   }
 219   if (!ok) {
 220     fprintf(stderr, "\n*** primeiter test failure:\n***   start = ");
 221     mp_writefile(start, stderr, 10);
 222     for (i = 0; i < N(pp); i++) {
 223       fprintf(stderr, "\n***   p[%d] = ", i);
 224       mp_writefile(ret[i], stderr, 10);
 225       fprintf(stderr, " %s ", MP_EQ(ret[i], pp[i]) ? "==" : "!=");
 226       mp_writefile(pp[i], stderr, 10);
 227     }
 228     fputc('\n', stderr);
 229   }
 230   for (i = 0; i < N(pp); i++) {
 231     MP_DROP(pp[i]);
 232     MP_DROP(ret[i]);
 233   }
 234   primeiter_destroy(&pi);
 235   MP_DROP(start);
 236   assert(mparena_count(MPARENA_GLOBAL) == 0);
 237   return (ok);
 238 }
 239
 240 static test_chunk tests[] = {
 241   { "primeiter", test,
 242     { &type_mp,  &type_mp, &type_mp, &type_mp, &type_mp, &type_mp, } },
 243   { 0 }
 244 };
 245
 246 int main(int argc, char *argv[])
 247 {
 248   test_run(argc, argv, tests, SRCDIR "/t/pgen");
 249   return (0);
 250 }
 251
 252 #endif
 253
 254 /*----- That's all, folks -------------------------------------------------*/