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 again:
 118   p = mprand(l->newp, pl, l->r, 1);
 119   pf.step = 2;
 120   p = pgen(l->u.s.name, p, p, l->iev, l->iec, 0, pgen_filter, &pf,
 121            rabin_iters(pl), pgen_test, &r);
 122   if (!p)
 123     goto again;
 124   f->p = p;
 125 }
 126
 127 static void llfree(limlee_factor *f, limlee_stepctx *l)
 128 {
 129   mp_drop(f->p);
 130 }
 131
 132 static const limlee_primeops primeops_simple = { llgen, llfree };
 133
 134 /*----- Lim-Lee stepper ---------------------------------------------------*/
 135
 136 /* --- @init@ --- *
 137  *
 138  * Arguments:   @pgen_event *ev@ = pointer to event block
 139  *              @limlee_stepctx *l@ = pointer to Lim-Lee context
 140  *
 141  * Returns:     A @PGEN@ result code.
 142  *
 143  * Use:         Initializes the stepper.
 144  */
 145
 146 static int init(pgen_event *ev, limlee_stepctx *l)
 147 {
 148   size_t i;
 149
 150   /* --- First of all, decide on a number of factors to make --- */
 151
 152   l->nf = l->pl / l->ql;
 153   if (l->nf < 2) return (PGEN_ABORT);
 154   l->nf--;
 155
 156   /* --- Now decide on how many primes I'll actually generate --- *
 157    *
 158    * The formula %$m = \max(3 n + 5, 25)$% comes from GPG's prime generation
 159    * library.
 160    */
 161
 162   l->poolsz = l->nf * 3 + 5;
 163   if (l->poolsz < 25)
 164     l->poolsz = 25;
 165
 166   /* --- Allocate and initialize the various tables --- */
 167
 168   l->c = xmalloc(l->poolsz);
 169   l->v = xmalloc(l->poolsz * sizeof(limlee_factor));
 170   comb_init(l->c, l->poolsz, l->nf);
 171   for (i = 0; i < l->poolsz; i++)
 172     l->v[i].p = 0;
 173
 174   /* --- Other bits of initialization --- */
 175
 176   l->seq = 0;
 177   if (!l->pops) {
 178     l->pops = &primeops_simple;
 179     l->pc = 0;
 180   }
 181
 182   /* --- Find a big prime later --- */
 183
 184   l->qq.p = 0;
 185
 186   return (PGEN_TRY);
 187 }
 188
 189 /* --- @next@ --- *
 190  *
 191  * Arguments:   @int rq@ = request which triggered this call
 192  *              @pgen_event *ev@ = pointer to event block
 193  *              @limlee_stepctx *l@ = pointer to Lim-Lee context
 194  *
 195  * Returns:     A @PGEN@ result code.
 196  *
 197  * Use:         Initializes the stepper.
 198  */
 199
 200 static int next(int rq, pgen_event *ev, limlee_stepctx *l)
 201 {
 202   dstr d = DSTR_INIT;
 203   mp *p;
 204   int rc;
 205   int dist;
 206   unsigned nb;
 207
 208   mp_drop(ev->m);
 209
 210   for (;;) {
 211     size_t i;
 212     mpmul mm = MPMUL_INIT;
 213
 214     /* --- Step on to next combination --- */
 215
 216     if (rq == PGEN_TRY && !comb_next(l->c, l->poolsz, l->nf)) {
 217       for (i = 0; i < l->poolsz; i++) {
 218         l->pops->pfree(&l->v[i], l);
 219         l->v[i].p = 0;
 220       }
 221     }
 222     rq = PGEN_TRY; /* For next time through */
 223
 224     /* --- If the large factor is performing badly, make a new one --- */
 225
 226     if (l->qq.p) {
 227       dist = l->u.s.disp < 0 ? -l->u.s.disp : l->u.s.disp;
 228       if (dist && dist > l->u.s.steps/dist) {
 229         l->pops->pfree(&l->qq, l);
 230         l->qq.p = 0;
 231       }
 232     }
 233
 234     /* --- Gather up some factors --- */
 235
 236     if (l->qq.p) mpmul_add(&mm, l->qq.p);
 237     for (i = 0; i < l->poolsz; i++) {
 238       if (!l->c[i])
 239         continue;
 240       if (!l->v[i].p) {
 241         DRESET(&d);
 242         dstr_putf(&d, "%s_%lu", ev->name, l->seq++);
 243         l->u.s.name = d.buf;
 244         l->pops->pgen(&l->v[i], l->ql, l);
 245       }
 246       mpmul_add(&mm, l->v[i].p);
 247     }
 248
 249     /* --- Check on the large factor --- */
 250
 251     p = mpmul_done(&mm);
 252     if (!l->qq.p) {
 253       DRESET(&d);
 254       dstr_putf(&d, "%s*_%lu", ev->name, l->seq++);
 255       l->u.s.name = d.buf;
 256       l->pops->pgen(&l->qq, l->pl - mp_bits(p), l);
 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);
 278   }
 279
 280   ev->m = p;
 281   DDESTROY(&d);
 282   return (rc);
 283 }
 284
 285 /* --- @done@ --- *
 286  *
 287  * Arguments:   @pgen_event *ev@ = pointer to event block
 288  *              @limlee_stepctx *l@ = pointer to Lim-Lee context
 289  *
 290  * Returns:     A @PGEN@ result code.
 291  *
 292  * Use:         Finalizes the stepper.  The output values in the context
 293  *              take on their final results; other resources are discarded.
 294  */
 295
 296 static int done(pgen_event *ev, limlee_stepctx *l)
 297 {
 298   size_t i, j;
 299   limlee_factor *v;
 300
 301   /* --- If an output vector of factors is wanted, produce one --- */
 302
 303   if (!(l->f & LIMLEE_KEEPFACTORS))
 304     v = 0;
 305   else {
 306     if (l->qq.p)
 307       l->nf++;
 308     v = xmalloc(l->nf * sizeof(limlee_factor));
 309   }
 310
 311   for (i = 0, j = 0; i < l->poolsz; i++) {
 312     if (v && l->c[i])
 313       v[j++] = l->v[i];
 314     else if (l->v[i].p)
 315       l->pops->pfree(&l->v[i], l);
 316   }
 317
 318   if (l->qq.p) {
 319     if (v)
 320       v[j++] = l->qq;
 321     else
 322       l->pops->pfree(&l->qq, l);
 323   }
 324
 325   xfree(l->v);
 326   l->v = v;
 327
 328   /* --- Free other resources --- */
 329
 330   xfree(l->c);
 331
 332   /* --- Done --- */
 333
 334   return (PGEN_DONE);
 335 }
 336
 337 /* --- @limlee_step@ --- */
 338
 339 int limlee_step(int rq, pgen_event *ev, void *p)
 340 {
 341   limlee_stepctx *l = p;
 342   int rc;
 343
 344   switch (rq) {
 345     case PGEN_BEGIN:
 346       if ((rc = init(ev, l)) != PGEN_TRY)
 347         return (rc);
 348     case PGEN_TRY:
 349       return (next(rq, ev, l));
 350     case PGEN_DONE:
 351       return (done(ev, l));
 352   }
 353   return (PGEN_ABORT);
 354 }
 355
 356 /*----- Main code ---------------------------------------------------------*/
 357
 358 /* --- @limlee@ --- *
 359  *
 360  * Arguments:   @const char *name@ = pointer to name root
 361  *              @mp *d@ = pointer to destination integer
 362  *              @mp *newp@ = how to generate factor primes
 363  *              @unsigned ql@ = size of individual factors
 364  *              @unsigned pl@ = size of large prime
 365  *              @grand *r@ = a random number source
 366  *              @unsigned on@ = number of outer attempts to make
 367  *              @pgen_proc *oev@ = outer event handler function
 368  *              @void *oec@ = argument for the outer event handler
 369  *              @pgen_proc *iev@ = inner event handler function
 370  *              @void *iec@ = argument for the inner event handler
 371  *              @size_t *nf@, @mp ***f@ = output array for factors
 372  *
 373  * Returns:     A Lim-Lee prime, or null if generation failed.
 374  *
 375  * Use:         Generates Lim-Lee primes.  A Lim-Lee prime %$p$% is one which
 376  *              satisfies %$p = 2 \prod_i q_i + 1$%, where all of the %$q_i$%
 377  *              are large enough to resist square-root discrete log
 378  *              algorithms.
 379  *
 380  *              If we succeed, and @f@ is non-null, we write the array of
 381  *              factors chosen to @f@ for the benefit of the caller.
 382  */
 383
 384 mp *limlee(const char *name, mp *d, mp *newp,
 385            unsigned ql, unsigned pl, grand *r,
 386            unsigned on, pgen_proc *oev, void *oec,
 387            pgen_proc *iev, void *iec,
 388            size_t *nf, mp ***f)
 389 {
 390   limlee_stepctx l;
 391   rabin rr;
 392
 393   l.f = 0; if (f) l.f |= LIMLEE_KEEPFACTORS;
 394   l.newp = newp;
 395   l.pl = pl; l.ql = ql;
 396   l.pops = 0;
 397   l.iev = iev;
 398   l.iec = iec;
 399   l.r = r;
 400
 401   d = pgen(name, d, 0, oev, oec, on, limlee_step, &l,
 402            rabin_iters(pl), pgen_test, &rr);
 403
 404   if (d && f) {
 405     mp **v;
 406     size_t i;
 407     v = xmalloc(l.nf * sizeof(mp *));
 408     for (i = 0; i < l.nf; i++)
 409       v[i] = l.v[i].p;
 410     xfree(l.v);
 411     *f = v;
 412     *nf = l.nf;
 413   }
 414
 415   return (d);
 416 }
 417
 418 /*----- That's all, folks -------------------------------------------------*/