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

   1 /* -*-c-*-
   2  *
   3  * Generate Lim-Lee primes
   4  *
   5  * (c) 2000 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 <mLib/alloc.h>
  31 #include <mLib/dstr.h>
  32
  33 #include "limlee.h"
  34 #include "mpmul.h"
  35 #include "mprand.h"
  36 #include "pgen.h"
  37 #include "rabin.h"
  38
  39 /*----- Stepping through combinations -------------------------------------*/
  40
  41 /* --- @comb_init@ --- *
  42  *
  43  * Arguments:   @octet *c@ = pointer to byte-flag array
  44  *              @unsigned n@ = number of items in the array
  45  *              @unsigned r@ = number of desired items
  46  *
  47  * Returns:     ---
  48  *
  49  * Use:         Initializes a byte-flag array which, under the control of
  50  *              @comb_next@, will step through all combinations of @r@ chosen
  51  *              elements.
  52  */
  53
  54 static void comb_init(octet *c, unsigned n, unsigned r)
  55 {
  56   memset(c, 0, n - r);
  57   memset(c + (n - r), 1, r);
  58 }
  59
  60 /* --- @comb_next@ --- *
  61  *
  62  * Arguments:   @octet *c@ = pointer to byte-flag array
  63  *              @unsigned n@ = number of items in the array
  64  *              @unsigned r@ = number of desired items
  65  *
  66  * Returns:     Nonzero if another combination was returned, zero if we've
  67  *              reached the end.
  68  *
  69  * Use:         Steps on to the next combination in sequence.
  70  */
  71
  72 static int comb_next(octet *c, unsigned n, unsigned r)
  73 {
  74   unsigned g = 0;
  75
  76   /* --- How the algorithm works --- *
  77    *
  78    * Set bits start at the end and work their way towards the start.
  79    * Excepting bits already at the start, we scan for the lowest set bit, and
  80    * move it one place nearer the start.  A group of bits at the start are
  81    * counted and reset just below the `moved' bit.  If there is no moved bit
  82    * then we're done.
  83    */
  84
  85   /* --- Count the group at the start --- */
  86
  87   for (; *c; c++) {
  88     g++;
  89     *c = 0;
  90   }
  91   if (g == r)
  92     return (0);
  93
  94   /* --- Move the next bit down one --- *
  95    *
  96    * There must be one, because otherwise we'd have counted %$r$% bits
  97    * earlier.
  98    */
  99
 100   for (; !*c; c++)
 101     ;
 102   *c = 0;
 103   g++;
 104   for (; g; g--)
 105     *--c = 1;
 106   return (1);
 107 }
 108
 109 /*----- Default prime generator -------------------------------------------*/
 110
 111 static void llgen(limlee_factor *f, unsigned pl, limlee_stepctx *l)
 112 {
 113   pgen_filterctx pf;
 114   rabin r;
 115   mp *p;
 116
 117   p = mprand(l->newp, pl, l->r, 1);
 118   pf.step = 2;
 119   p = pgen(l->u.s.name, p, p, l->iev, l->iec, 0, pgen_filter, &pf,
 120            rabin_iters(pl), pgen_test, &r);
 121   f->p = p;
 122 }
 123
 124 static void llfree(limlee_factor *f, limlee_stepctx *l)
 125 {
 126   mp_drop(f->p);
 127 }
 128
 129 static const limlee_primeops primeops_simple = { llgen, llfree };
 130
 131 /*----- Lim-Lee stepper ---------------------------------------------------*/
 132
 133 /* --- @init@ --- *
 134  *
 135  * Arguments:   @pgen_event *ev@ = pointer to event block
 136  *              @limlee_stepctx *l@ = pointer to Lim-Lee context
 137  *
 138  * Returns:     A @PGEN@ result code.
 139  *
 140  * Use:         Initializes the stepper.
 141  */
 142
 143 static int init(pgen_event *ev, limlee_stepctx *l)
 144 {
 145   size_t i;
 146
 147   /* --- First of all, decide on a number of factors to make --- */
 148
 149   l->nf = l->pl / l->ql;
 150   if (l->nf < 2) return (PGEN_ABORT);
 151   l->nf--;
 152
 153   /* --- Now decide on how many primes I'll actually generate --- *
 154    *
 155    * The formula %$m = \max(3 n + 5, 25)$% comes from GPG's prime generation
 156    * library.
 157    */
 158
 159   l->poolsz = l->nf * 3 + 5;
 160   if (l->poolsz < 25)
 161     l->poolsz = 25;
 162
 163   /* --- Allocate and initialize the various tables --- */
 164
 165   l->c = xmalloc(l->poolsz);
 166   l->v = xmalloc(l->poolsz * sizeof(limlee_factor));
 167   comb_init(l->c, l->poolsz, l->nf);
 168   for (i = 0; i < l->poolsz; i++)
 169     l->v[i].p = 0;
 170
 171   /* --- Other bits of initialization --- */
 172
 173   l->seq = 0;
 174   if (!l->pops) {
 175     l->pops = &primeops_simple;
 176     l->pc = 0;
 177   }
 178
 179   /* --- Find a big prime later --- */
 180
 181   l->qq.p = 0;
 182
 183   return (PGEN_TRY);
 184 }
 185
 186 /* --- @next@ --- *
 187  *
 188  * Arguments:   @int rq@ = request which triggered this call
 189  *              @pgen_event *ev@ = pointer to event block
 190  *              @limlee_stepctx *l@ = pointer to Lim-Lee context
 191  *
 192  * Returns:     A @PGEN@ result code.
 193  *
 194  * Use:         Initializes the stepper.
 195  */
 196
 197 static int next(int rq, pgen_event *ev, limlee_stepctx *l)
 198 {
 199   dstr d = DSTR_INIT;
 200   mp *p = 0;
 201   int rc;
 202   int dist;
 203   unsigned nb;
 204
 205   mp_drop(ev->m);
 206
 207   for (;;) {
 208     size_t i;
 209     mpmul mm = MPMUL_INIT;
 210
 211     /* --- Step on to next combination --- */
 212
 213     if (rq == PGEN_TRY && !comb_next(l->c, l->poolsz, l->nf)) {
 214       for (i = 0; i < l->poolsz; i++) {
 215         l->pops->pfree(&l->v[i], l);
 216         l->v[i].p = 0;
 217       }
 218     }
 219     rq = PGEN_TRY; /* For next time through */
 220
 221     /* --- If the large factor is performing badly, make a new one --- */
 222
 223     if (l->qq.p) {
 224       dist = l->u.s.disp < 0 ? -l->u.s.disp : l->u.s.disp;
 225       if (dist && dist > l->u.s.steps/3) {
 226         l->pops->pfree(&l->qq, l);
 227         l->qq.p = 0;
 228       }
 229     }
 230
 231     /* --- Gather up some factors --- */
 232
 233     if (l->qq.p) mpmul_add(&mm, l->qq.p);
 234     for (i = 0; i < l->poolsz; i++) {
 235       if (!l->c[i])
 236         continue;
 237       if (!l->v[i].p) {
 238         DRESET(&d);
 239         dstr_putf(&d, "%s_%lu", ev->name, l->seq++);
 240         l->u.s.name = d.buf;
 241         l->pops->pgen(&l->v[i], l->ql, l);
 242         if (!l->v[i].p)
 243           { mp_drop(mpmul_done(&mm)); rc = PGEN_ABORT; goto end; }
 244       }
 245       mpmul_add(&mm, l->v[i].p);
 246     }
 247
 248     /* --- Check on the large factor --- */
 249
 250     p = mpmul_done(&mm);
 251     if (!l->qq.p) {
 252       DRESET(&d);
 253       dstr_putf(&d, "%s*_%lu", ev->name, l->seq++);
 254       l->u.s.name = d.buf;
 255       l->pops->pgen(&l->qq, l->pl - mp_bits(p), l);
 256       if (!l->qq.p) { MP_DROP(p); p = 0; rc = PGEN_ABORT; break; }
 257       l->u.s.steps = l->u.s.disp = 0;
 258       p = mp_mul(p, p, l->qq.p);
 259     }
 260     p = mp_lsl(p, p, 1);
 261     p->v[0] |= 1;
 262
 263     nb = mp_bits(p);
 264     l->u.s.steps++;
 265     if (nb < l->pl) {
 266       l->u.s.disp--;
 267       continue;
 268     } else if (nb > l->pl) {
 269       l->u.s.disp++;
 270       continue;
 271     }
 272
 273     /* --- Check it for small factors --- */
 274
 275     if ((rc = pfilt_smallfactor(p)) != PGEN_FAIL)
 276       break;
 277     MP_DROP(p); p = 0;
 278   }
 279
 280 end:
 281   ev->m = p;
 282   DDESTROY(&d);
 283   return (rc);
 284 }
 285
 286 /* --- @done@ --- *
 287  *
 288  * Arguments:   @pgen_event *ev@ = pointer to event block
 289  *              @limlee_stepctx *l@ = pointer to Lim-Lee context
 290  *
 291  * Returns:     A @PGEN@ result code.
 292  *
 293  * Use:         Finalizes the stepper.  The output values in the context
 294  *              take on their final results; other resources are discarded.
 295  */
 296
 297 static int done(pgen_event *ev, limlee_stepctx *l)
 298 {
 299   size_t i, j;
 300   limlee_factor *v;
 301
 302   /* --- If an output vector of factors is wanted, produce one --- */
 303
 304   if (!(l->f & LIMLEE_KEEPFACTORS))
 305     v = 0;
 306   else {
 307     if (l->qq.p)
 308       l->nf++;
 309     v = xmalloc(l->nf * sizeof(limlee_factor));
 310   }
 311
 312   for (i = 0, j = 0; i < l->poolsz; i++) {
 313     if (v && l->c[i])
 314       v[j++] = l->v[i];
 315     else if (l->v[i].p)
 316       l->pops->pfree(&l->v[i], l);
 317   }
 318
 319   if (l->qq.p) {
 320     if (v)
 321       v[j++] = l->qq;
 322     else
 323       l->pops->pfree(&l->qq, l);
 324   }
 325
 326   xfree(l->v);
 327   l->v = v;
 328
 329   /* --- Free other resources --- */
 330
 331   xfree(l->c);
 332
 333   /* --- Done --- */
 334
 335   return (PGEN_DONE);
 336 }
 337
 338 /* --- @limlee_step@ --- */
 339
 340 int limlee_step(int rq, pgen_event *ev, void *p)
 341 {
 342   limlee_stepctx *l = p;
 343   int rc;
 344
 345   switch (rq) {
 346     case PGEN_BEGIN:
 347       if ((rc = init(ev, l)) != PGEN_TRY)
 348         return (rc);
 349     case PGEN_TRY:
 350       return (next(rq, ev, l));
 351     case PGEN_DONE:
 352       return (done(ev, l));
 353   }
 354   return (PGEN_ABORT);
 355 }
 356
 357 /*----- Main code ---------------------------------------------------------*/
 358
 359 /* --- @limlee@ --- *
 360  *
 361  * Arguments:   @const char *name@ = pointer to name root
 362  *              @mp *d@ = pointer to destination integer
 363  *              @mp *newp@ = how to generate factor primes
 364  *              @unsigned ql@ = size of individual factors
 365  *              @unsigned pl@ = size of large prime
 366  *              @grand *r@ = a random number source
 367  *              @unsigned on@ = number of outer attempts to make
 368  *              @pgen_proc *oev@ = outer event handler function
 369  *              @void *oec@ = argument for the outer event handler
 370  *              @pgen_proc *iev@ = inner event handler function
 371  *              @void *iec@ = argument for the inner event handler
 372  *              @size_t *nf@, @mp ***f@ = output array for factors
 373  *
 374  * Returns:     A Lim-Lee prime, or null if generation failed.
 375  *
 376  * Use:         Generates Lim-Lee primes.  A Lim-Lee prime %$p$% is one which
 377  *              satisfies %$p = 2 \prod_i q_i + 1$%, where all of the %$q_i$%
 378  *              are large enough to resist square-root discrete log
 379  *              algorithms.
 380  *
 381  *              If we succeed, and @f@ is non-null, we write the array of
 382  *              factors chosen to @f@ for the benefit of the caller.
 383  */
 384
 385 mp *limlee(const char *name, mp *d, mp *newp,
 386            unsigned ql, unsigned pl, grand *r,
 387            unsigned on, pgen_proc *oev, void *oec,
 388            pgen_proc *iev, void *iec,
 389            size_t *nf, mp ***f)
 390 {
 391   limlee_stepctx l;
 392   rabin rr;
 393
 394   l.f = 0; if (f) l.f |= LIMLEE_KEEPFACTORS;
 395   l.newp = newp;
 396   l.pl = pl; l.ql = ql;
 397   l.pops = 0;
 398   l.iev = iev;
 399   l.iec = iec;
 400   l.r = r;
 401
 402   d = pgen(name, d, 0, oev, oec, on, limlee_step, &l,
 403            rabin_iters(pl), pgen_test, &rr);
 404
 405   if (d && f) {
 406     mp **v;
 407     size_t i;
 408     v = xmalloc(l.nf * sizeof(mp *));
 409     for (i = 0; i < l.nf; i++)
 410       v[i] = l.v[i].p;
 411     xfree(l.v);
 412     *f = v;
 413     *nf = l.nf;
 414   }
 415
 416   return (d);
 417 }
 418
 419 /*----- That's all, folks -------------------------------------------------*/