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

   1 /* -*-c-*-
   2  *
   3  * $Id: limlee.c,v 1.8 2001/02/03 11:59:07 mdw Exp $
   4  *
   5  * Generate Lim-Lee primes
   6  *
   7  * (c) 2000 Straylight/Edgeware
   8  */
   9
  10 /*----- Licensing notice --------------------------------------------------*
  11  *
  12  * This file is part of Catacomb.
  13  *
  14  * Catacomb is free software; you can redistribute it and/or modify
  15  * it under the terms of the GNU Library General Public License as
  16  * published by the Free Software Foundation; either version 2 of the
  17  * License, or (at your option) any later version.
  18  *
  19  * Catacomb is distributed in the hope that it will be useful,
  20  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  21  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  22  * GNU Library General Public License for more details.
  23  *
  24  * You should have received a copy of the GNU Library General Public
  25  * License along with Catacomb; if not, write to the Free
  26  * Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
  27  * MA 02111-1307, USA.
  28  */
  29
  30 /*----- Revision history --------------------------------------------------*
  31  *
  32  * $Log: limlee.c,v $
  33  * Revision 1.8  2001/02/03 11:59:07  mdw
  34  * Don't use the @pgen@ random number generator for generating primes: it's
  35  * only for testing them.  Use a caller-supplied one instead.
  36  *
  37  * Revision 1.7  2001/01/25 21:40:44  mdw
  38  * Remove dead code now that the new stepper structure is trustworthy.
  39  *
  40  * Revision 1.6  2001/01/25 21:16:20  mdw
  41  * Boring cosmetic stuff.
  42  *
  43  * Revision 1.5  2000/08/18 19:16:51  mdw
  44  * New stepper interface for constructing Lim-Lee primes.
  45  *
  46  * Revision 1.4  2000/08/15 21:45:05  mdw
  47  * Use the new trial division equipment in pfilt.  This gives a 10%
  48  * performance improvement in dsa-gen.t.
  49  *
  50  * Revision 1.3  2000/07/29 09:58:32  mdw
  51  * (limlee): Bug fix.  Old versions didn't set the filter step if @ql@ was
  52  * an exact divisor of @pl@.
  53  *
  54  * Revision 1.2  2000/07/26 18:00:00  mdw
  55  * No footer line!
  56  *
  57  * Revision 1.1  2000/07/09 21:30:58  mdw
  58  * Lim-Lee prime generation.
  59  *
  60  */
  61
  62 /*----- Header files ------------------------------------------------------*/
  63
  64 #include <mLib/alloc.h>
  65 #include <mLib/dstr.h>
  66
  67 #include "limlee.h"
  68 #include "mpmul.h"
  69 #include "mprand.h"
  70 #include "pgen.h"
  71 #include "rabin.h"
  72
  73 /*----- Stepping through combinations -------------------------------------*/
  74
  75 /* --- @comb_init@ --- *
  76  *
  77  * Arguments:   @octet *c@ = pointer to byte-flag array
  78  *              @unsigned n@ = number of items in the array
  79  *              @unsigned r@ = number of desired items
  80  *
  81  * Returns:     ---
  82  *
  83  * Use:         Initializes a byte-flag array which, under the control of
  84  *              @comb_next@, will step through all combinations of @r@ chosen
  85  *              elements.
  86  */
  87
  88 static void comb_init(octet *c, unsigned n, unsigned r)
  89 {
  90   memset(c, 0, n - r);
  91   memset(c + (n - r), 1, r);
  92 }
  93
  94 /* --- @comb_next@ --- *
  95  *
  96  * Arguments:   @octet *c@ = pointer to byte-flag array
  97  *              @unsigned n@ = number of items in the array
  98  *              @unsigned r@ = number of desired items
  99  *
 100  * Returns:     Nonzero if another combination was returned, zero if we've
 101  *              reached the end.
 102  *
 103  * Use:         Steps on to the next combination in sequence.
 104  */
 105
 106 static int comb_next(octet *c, unsigned n, unsigned r)
 107 {
 108   unsigned g = 0;
 109
 110   /* --- How the algorithm works --- *
 111    *
 112    * Set bits start at the end and work their way towards the start.
 113    * Excepting bits already at the start, we scan for the lowest set bit, and
 114    * move it one place nearer the start.  A group of bits at the start are
 115    * counted and reset just below the `moved' bit.  If there is no moved bit
 116    * then we're done.
 117    */
 118
 119   /* --- Count the group at the start --- */
 120
 121   for (; *c; c++) {
 122     g++;
 123     *c = 0;
 124   }
 125   if (g == r)
 126     return (0);
 127
 128   /* --- Move the next bit down one --- *
 129    *
 130    * There must be one, because otherwise we'd have counted %$r$% bits
 131    * earlier.
 132    */
 133
 134   for (; !*c; c++)
 135     ;
 136   *c = 0;
 137   g++;
 138   for (; g; g--)
 139     *--c = 1;
 140   return (1);
 141 }
 142
 143 /*----- Default prime generator -------------------------------------------*/
 144
 145 static void llgen(limlee_factor *f, unsigned pl, limlee_stepctx *l)
 146 {
 147   pgen_filterctx pf;
 148   rabin r;
 149   mp *p;
 150
 151 again:
 152   p = mprand(l->newp, pl, l->r, 1);
 153   pf.step = 2;
 154   p = pgen(l->d.buf, p, p, l->iev, l->iec, 0, pgen_filter, &pf,
 155            rabin_iters(pl), pgen_test, &r);
 156   if (!p)
 157     goto again;
 158   f->p = p;
 159 }
 160
 161 static void llfree(limlee_factor *f, limlee_stepctx *l)
 162 {
 163   mp_drop(f->p);
 164 }
 165
 166 static const limlee_primeops primeops_simple = { llgen, llfree };
 167
 168 /*----- Lim-Lee stepper ---------------------------------------------------*/
 169
 170 /* --- @init@ --- *
 171  *
 172  * Arguments:   @pgen_event *ev@ = pointer to event block
 173  *              @limlee_stepctx *l@ = pointer to Lim-Lee context
 174  *
 175  * Returns:     A @PGEN@ result code.
 176  *
 177  * Use:         Initializes the stepper.
 178  */
 179
 180 static int init(pgen_event *ev, limlee_stepctx *l)
 181 {
 182   size_t i;
 183   unsigned qql;
 184
 185   /* --- First of all, decide on a number of factors to make --- */
 186
 187   l->nf = l->pl / l->ql;
 188   qql = l->pl % l->ql;
 189   if (!l->nf)
 190     return (PGEN_ABORT);
 191   else if (qql && l->nf > 1) {
 192     l->nf--;
 193     qql += l->ql;
 194   }
 195
 196   /* --- Now decide on how many primes I'll actually generate --- *
 197    *
 198    * The formula %$m = \max(3 n + 5, 25)$% comes from GPG's prime generation
 199    * library.
 200    */
 201
 202   l->poolsz = l->nf * 3 + 5;
 203   if (l->poolsz < 25)
 204     l->poolsz = 25;
 205
 206   /* --- Allocate and initialize the various tables --- */
 207
 208   l->c = xmalloc(l->poolsz);
 209   l->v = xmalloc(l->poolsz * sizeof(limlee_factor));
 210   comb_init(l->c, l->poolsz, l->nf);
 211   for (i = 0; i < l->poolsz; i++)
 212     l->v[i].p = 0;
 213
 214   /* --- Other bits of initialization --- */
 215
 216   l->seq = 0;
 217   dstr_create(&l->d);
 218   if (!l->pops) {
 219     l->pops = &primeops_simple;
 220     l->pc = 0;
 221   }
 222
 223   /* --- Find a big prime --- */
 224
 225   if (!qql)
 226     l->qq.p = 0;
 227   else {
 228     dstr_putf(&l->d, "%s*", ev->name);
 229     l->pops->pgen(&l->qq, qql, l);
 230   }
 231
 232   return (PGEN_TRY);
 233 }
 234
 235 /* --- @next@ --- *
 236  *
 237  * Arguments:   @int rq@ = request which triggered this call
 238  *              @pgen_event *ev@ = pointer to event block
 239  *              @limlee_stepctx *l@ = pointer to Lim-Lee context
 240  *
 241  * Returns:     A @PGEN@ result code.
 242  *
 243  * Use:         Initializes the stepper.
 244  */
 245
 246 static int next(int rq, pgen_event *ev, limlee_stepctx *l)
 247 {
 248   mp *p;
 249   int rc;
 250
 251   mp_drop(ev->m);
 252
 253   for (;;) {
 254     size_t i;
 255     mpmul mm = MPMUL_INIT;
 256
 257     /* --- Step on to next combination --- */
 258
 259     if (rq == PGEN_TRY && !comb_next(l->c, l->poolsz, l->nf)) {
 260       for (i = 0; i < l->poolsz; i++) {
 261         l->pops->pfree(&l->v[i], l);
 262         l->v[i].p = 0;
 263       }
 264     }
 265     rq = PGEN_TRY; /* For next time through */
 266
 267     /* --- Gather up some factors --- */
 268
 269     if (l->qq.p)
 270       mpmul_add(&mm, l->qq.p);
 271     for (i = 0; i < l->poolsz; i++) {
 272       if (!l->c[i])
 273         continue;
 274       if (!l->v[i].p) {
 275         DRESET(&l->d);
 276         dstr_putf(&l->d, "%s_%lu", ev->name, l->seq++);
 277         l->pops->pgen(&l->v[i], l->ql, l);
 278       }
 279       mpmul_add(&mm, l->v[i].p);
 280     }
 281
 282     /* --- Check it for small factors --- */
 283
 284     p = mpmul_done(&mm);
 285     p = mp_lsl(p, p, 1);
 286     p->v[0] |= 1;
 287     if ((rc = pfilt_smallfactor(p)) != PGEN_FAIL)
 288       break;
 289     mp_drop(p);
 290   }
 291
 292   ev->m = p;
 293   return (rc);
 294 }
 295
 296 /* --- @done@ --- *
 297  *
 298  * Arguments:   @pgen_event *ev@ = pointer to event block
 299  *              @limlee_stepctx *l@ = pointer to Lim-Lee context
 300  *
 301  * Returns:     A @PGEN@ result code.
 302  *
 303  * Use:         Finalizes the stepper.  The output values in the context
 304  *              take on their final results; other resources are discarded.
 305  */
 306
 307 static int done(pgen_event *ev, limlee_stepctx *l)
 308 {
 309   size_t i, j;
 310   limlee_factor *v;
 311
 312   /* --- If an output vector of factors is wanted, produce one --- */
 313
 314   if (!(l->f & LIMLEE_KEEPFACTORS))
 315     v = 0;
 316   else {
 317     if (l->qq.p)
 318       l->nf++;
 319     v = xmalloc(l->nf * sizeof(limlee_factor));
 320   }
 321
 322   for (i = 0, j = 0; i < l->poolsz; i++) {
 323     if (v && l->c[i])
 324       v[j++] = l->v[i];
 325     else if (l->v[i].p)
 326       l->pops->pfree(&l->v[i], l);
 327   }
 328
 329   if (l->qq.p) {
 330     if (v)
 331       v[j++] = l->qq;
 332     else
 333       l->pops->pfree(&l->qq, l);
 334   }
 335
 336   xfree(l->v);
 337   l->v = v;
 338
 339   /* --- Free other resources --- */
 340
 341   xfree(l->c);
 342   dstr_destroy(&l->d);
 343
 344   /* --- Done --- */
 345
 346   return (PGEN_DONE);
 347 }
 348
 349 /* --- @limlee_step@ --- */
 350
 351 int limlee_step(int rq, pgen_event *ev, void *p)
 352 {
 353   limlee_stepctx *l = p;
 354   int rc;
 355
 356   switch (rq) {
 357     case PGEN_BEGIN:
 358       if ((rc = init(ev, l)) != PGEN_TRY)
 359         return (rc);
 360     case PGEN_TRY:
 361       return (next(rq, ev, l));
 362     case PGEN_DONE:
 363       return (done(ev, l));
 364   }
 365   return (PGEN_ABORT);
 366 }
 367
 368 /*----- Main code ---------------------------------------------------------*/
 369
 370 /* --- @limlee@ --- *
 371  *
 372  * Arguments:   @const char *name@ = pointer to name root
 373  *              @mp *d@ = pointer to destination integer
 374  *              @mp *newp@ = how to generate factor primes
 375  *              @unsigned ql@ = size of individual factors
 376  *              @unsigned pl@ = size of large prime
 377  *              @grand *r@ = a random number source
 378  *              @unsigned on@ = number of outer attempts to make
 379  *              @pgen_proc *oev@ = outer event handler function
 380  *              @void *oec@ = argument for the outer event handler
 381  *              @pgen_proc *iev@ = inner event handler function
 382  *              @void *iec@ = argument for the inner event handler
 383  *              @size_t *nf@, @mp ***f@ = output array for factors
 384  *
 385  * Returns:     A Lim-Lee prime, or null if generation failed.
 386  *
 387  * Use:         Generates Lim-Lee primes.  A Lim-Lee prime %$p$% is one which
 388  *              satisfies %$p = 2 \prod_i q_i + 1$%, where all of the %$q_i$%
 389  *              are large enough to resist square-root discrete log
 390  *              algorithms.
 391  *
 392  *              If we succeed, and @f@ is non-null, we write the array of
 393  *              factors chosen to @f@ for the benefit of the caller.
 394  */
 395
 396 mp *limlee(const char *name, mp *d, mp *newp,
 397            unsigned ql, unsigned pl, grand *r,
 398            unsigned on, pgen_proc *oev, void *oec,
 399            pgen_proc *iev, void *iec,
 400            size_t *nf, mp ***f)
 401 {
 402   limlee_stepctx l;
 403   rabin rr;
 404
 405   l.f = 0; if (f) l.f |= LIMLEE_KEEPFACTORS;
 406   l.newp = newp;
 407   l.pl = pl; l.ql = ql;
 408   l.pops = 0;
 409   l.iev = iev;
 410   l.iec = iec;
 411   l.r = r;
 412
 413   d = pgen(name, d, 0, oev, oec, on, limlee_step, &l,
 414            rabin_iters(pl), pgen_test, &rr);
 415
 416   if (f) {
 417     mp **v;
 418     size_t i;
 419     v = xmalloc(l.nf * sizeof(mp *));
 420     for (i = 0; i < l.nf; i++)
 421       v[i] = l.v[i].p;
 422     xfree(l.v);
 423     *f = v;
 424     *nf = l.nf;
 425   }
 426
 427   return (d);
 428 }
 429
 430 /*----- That's all, folks -------------------------------------------------*/