git.distorted.org.uk Git - u/mdw/catacomb/blob - limlee.c

   1 /* -*-c-*-
   2  *
   3  * $Id: limlee.c,v 1.5 2000/08/18 19:16:51 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.5  2000/08/18 19:16:51  mdw
  34  * New stepper interface for constructing Lim-Lee primes.
  35  *
  36  * Revision 1.4  2000/08/15 21:45:05  mdw
  37  * Use the new trial division equipment in pfilt.  This gives a 10%
  38  * performance improvement in dsa-gen.t.
  39  *
  40  * Revision 1.3  2000/07/29 09:58:32  mdw
  41  * (limlee): Bug fix.  Old versions didn't set the filter step if @ql@ was
  42  * an exact divisor of @pl@.
  43  *
  44  * Revision 1.2  2000/07/26 18:00:00  mdw
  45  * No footer line!
  46  *
  47  * Revision 1.1  2000/07/09 21:30:58  mdw
  48  * Lim-Lee prime generation.
  49  *
  50  */
  51
  52 /*----- Header files ------------------------------------------------------*/
  53
  54 #include <mLib/alloc.h>
  55 #include <mLib/dstr.h>
  56
  57 #include "limlee.h"
  58 #include "mpmul.h"
  59 #include "mprand.h"
  60 #include "pgen.h"
  61 #include "rabin.h"
  62
  63 /*----- Stepping through combinations -------------------------------------*/
  64
  65 /* --- @comb_init@ --- *
  66  *
  67  * Arguments:   @octet *c@ = pointer to byte-flag array
  68  *              @unsigned n@ = number of items in the array
  69  *              @unsigned r@ = number of desired items
  70  *
  71  * Returns:     ---
  72  *
  73  * Use:         Initializes a byte-flag array which, under the control of
  74  *              @comb_next@, will step through all combinations of @r@ chosen
  75  *              elements.
  76  */
  77
  78 static void comb_init(octet *c, unsigned n, unsigned r)
  79 {
  80   memset(c, 0, n - r);
  81   memset(c + (n - r), 1, r);
  82 }
  83
  84 /* --- @comb_next@ --- *
  85  *
  86  * Arguments:   @octet *c@ = pointer to byte-flag array
  87  *              @unsigned n@ = number of items in the array
  88  *              @unsigned r@ = number of desired items
  89  *
  90  * Returns:     Nonzero if another combination was returned, zero if we've
  91  *              reached the end.
  92  *
  93  * Use:         Steps on to the next combination in sequence.
  94  */
  95
  96 static int comb_next(octet *c, unsigned n, unsigned r)
  97 {
  98   unsigned g = 0;
  99
 100   /* --- How the algorithm works --- *
 101    *
 102    * Set bits start at the end and work their way towards the start.
 103    * Excepting bits already at the start, we scan for the lowest set bit, and
 104    * move it one place nearer the start.  A group of bits at the start are
 105    * counted and reset just below the `moved' bit.  If there is no moved bit
 106    * then we're done.
 107    */
 108
 109   /* --- Count the group at the start --- */
 110
 111   for (; *c; c++) {
 112     g++;
 113     *c = 0;
 114   }
 115   if (g == r)
 116     return (0);
 117
 118   /* --- Move the next bit down one --- *
 119    *
 120    * There must be one, because otherwise we'd have counted %$r$% bits
 121    * earlier.
 122    */
 123
 124   for (; !*c; c++)
 125     ;
 126   *c = 0;
 127   g++;
 128   for (; g; g--)
 129     *--c = 1;
 130   return (1);
 131 }
 132
 133 /*----- Default prime generator -------------------------------------------*/
 134
 135 static void llgen(limlee_factor *f, unsigned pl, limlee_stepctx *l)
 136 {
 137   pgen_filterctx pf;
 138   rabin r;
 139   mp *p;
 140
 141 again:
 142   p = mprand(l->newp, pl, l->r, 1);
 143   pf.step = 2;
 144   p = pgen(l->d.buf, p, p, l->iev, l->iec, 0, pgen_filter, &pf,
 145            rabin_iters(pl), pgen_test, &r);
 146   if (!p)
 147     goto again;
 148   f->p = p;
 149 }
 150
 151 static void llfree(limlee_factor *f, limlee_stepctx *l)
 152 {
 153   if (f->p)
 154     mp_drop(f->p);
 155 }
 156
 157 static const limlee_primeops primeops_simple = { llgen, llfree };
 158
 159 /*----- Lim-Lee stepper ---------------------------------------------------*/
 160
 161 /* --- @init@ --- *
 162  *
 163  * Arguments:   @pgen_event *ev@ = pointer to event block
 164  *              @limlee_stepctx *l@ = pointer to Lim-Lee context
 165  *
 166  * Returns:     A @PGEN@ result code.
 167  *
 168  * Use:         Initializes the stepper.
 169  */
 170
 171 static int init(pgen_event *ev, limlee_stepctx *l)
 172 {
 173   size_t i;
 174   unsigned qql;
 175
 176   /* --- First of all, decide on a number of factors to make --- */
 177
 178   l->nf = l->pl / l->ql;
 179   qql = l->pl % l->ql;
 180   if (!l->nf)
 181     return (PGEN_ABORT);
 182   else if (qql && l->nf > 1) {
 183     l->nf--;
 184     qql += l->ql;
 185   }
 186
 187   /* --- Now decide on how many primes I'll actually generate --- *
 188    *
 189    * The formula %$m = \max(3 n + 5, 25)$% comes from GPG's prime generation
 190    * library.
 191    */
 192
 193   l->poolsz = l->nf * 3 + 5;
 194   if (l->poolsz < 25)
 195     l->poolsz = 25;
 196
 197   /* --- Allocate and initialize the various tables --- */
 198
 199   l->c = xmalloc(l->poolsz);
 200   l->v = xmalloc(l->poolsz * sizeof(limlee_factor));
 201   comb_init(l->c, l->poolsz, l->nf);
 202   for (i = 0; i < l->poolsz; i++)
 203     l->v[i].p = 0;
 204
 205   /* --- Other bits of initialization --- */
 206
 207   l->seq = 0;
 208   l->r = ev->r;
 209   dstr_create(&l->d);
 210   if (!l->pops) {
 211     l->pops = &primeops_simple;
 212     l->pc = 0;
 213   }
 214
 215   /* --- Find a big prime --- */
 216
 217   if (!qql)
 218     l->qq.p = 0;
 219   else {
 220     dstr_putf(&l->d, "%s*", ev->name);
 221     l->pops->pgen(&l->qq, qql, l);
 222   }
 223
 224   return (PGEN_TRY);
 225 }
 226
 227 /* --- @next@ --- *
 228  *
 229  * Arguments:   @int rq@ = request which triggered this call
 230  *              @pgen_event *ev@ = pointer to event block
 231  *              @limlee_stepctx *l@ = pointer to Lim-Lee context
 232  *
 233  * Returns:     A @PGEN@ result code.
 234  *
 235  * Use:         Initializes the stepper.
 236  */
 237
 238 static int next(int rq, pgen_event *ev, limlee_stepctx *l)
 239 {
 240   mp *p;
 241   int rc;
 242
 243   if (ev->m)
 244     mp_drop(ev->m);
 245   l->r = ev->r;
 246
 247   for (;;) {
 248     size_t i;
 249     mpmul mm = MPMUL_INIT;
 250
 251     /* --- Step on to next combination --- */
 252
 253     if (rq == PGEN_TRY && !comb_next(l->c, l->poolsz, l->nf)) {
 254       for (i = 0; i < l->poolsz; i++) {
 255         l->pops->pfree(&l->v[i], l);
 256         l->v[i].p = 0;
 257       }
 258     }
 259     rq = PGEN_TRY; /* For next time through */
 260
 261     /* --- Gather up some factors --- */
 262
 263     if (l->qq.p)
 264       mpmul_add(&mm, l->qq.p);
 265     for (i = 0; i < l->poolsz; i++) {
 266       if (!l->c[i])
 267         continue;
 268       if (!l->v[i].p) {
 269         DRESET(&l->d);
 270         dstr_putf(&l->d, "%s_%lu", ev->name, l->seq++);
 271         l->pops->pgen(&l->v[i], l->ql, l);
 272       }
 273       mpmul_add(&mm, l->v[i].p);
 274     }
 275
 276     /* --- Check it for small factors --- */
 277
 278     p = mpmul_done(&mm);
 279     p = mp_lsl(p, p, 1);
 280     p->v[0] |= 1;
 281     if ((rc = pfilt_smallfactor(p)) != PGEN_FAIL)
 282       break;
 283     mp_drop(p);
 284   }
 285
 286   ev->m = p;
 287   return (rc);
 288 }
 289
 290 /* --- @done@ --- *
 291  *
 292  * Arguments:   @pgen_event *ev@ = pointer to event block
 293  *              @limlee_stepctx *l@ = pointer to Lim-Lee context
 294  *
 295  * Returns:     A @PGEN@ result code.
 296  *
 297  * Use:         Finalizes the stepper.  The output values in the context
 298  *              take on their final results; other resources are discarded.
 299  */
 300
 301 static int done(pgen_event *ev, limlee_stepctx *l)
 302 {
 303   size_t i, j;
 304   limlee_factor *v;
 305
 306   /* --- If an output vector of factors is wanted, produce one --- */
 307
 308   if (!(l->f & LIMLEE_KEEPFACTORS))
 309     v = 0;
 310   else {
 311     if (l->qq.p)
 312       l->nf++;
 313     v = xmalloc(l->nf * sizeof(limlee_factor));
 314   }
 315
 316   for (i = 0, j = 0; i < l->poolsz; i++) {
 317     if (v && l->c[i])
 318       v[j++] = l->v[i];
 319     else if (l->v[i].p)
 320       l->pops->pfree(&l->v[i], l);
 321   }
 322
 323   if (l->qq.p) {
 324     if (v)
 325       v[j++] = l->qq;
 326     else
 327       l->pops->pfree(&l->qq, l);
 328   }
 329
 330   xfree(l->v);
 331   l->v = v;
 332
 333   /* --- Free other resources --- */
 334
 335   xfree(l->c);
 336   dstr_destroy(&l->d);
 337
 338   /* --- Done --- */
 339
 340   return (PGEN_DONE);
 341 }
 342
 343 /* --- @limlee_step@ --- */
 344
 345 int limlee_step(int rq, pgen_event *ev, void *p)
 346 {
 347   limlee_stepctx *l = p;
 348   int rc;
 349
 350   switch (rq) {
 351     case PGEN_BEGIN:
 352       if ((rc = init(ev, l)) != PGEN_TRY)
 353         return (rc);
 354     case PGEN_TRY:
 355       return (next(rq, ev, l));
 356     case PGEN_DONE:
 357       return (done(ev, l));
 358   }
 359   return (PGEN_ABORT);
 360 }
 361
 362 /*----- Main code ---------------------------------------------------------*/
 363
 364 /* --- @limlee@ --- *
 365  *
 366  * Arguments:   @const char *name@ = pointer to name root
 367  *              @mp *d@ = pointer to destination integer
 368  *              @mp *newp@ = how to generate factor primes
 369  *              @unsigned ql@ = size of individual factors
 370  *              @unsigned pl@ = size of large prime
 371  *              @grand *r@ = a random number source
 372  *              @unsigned on@ = number of outer attempts to make
 373  *              @pgen_proc *oev@ = outer event handler function
 374  *              @void *oec@ = argument for the outer event handler
 375  *              @pgen_proc *iev@ = inner event handler function
 376  *              @void *iec@ = argument for the inner event handler
 377  *              @size_t *nf@, @mp ***f@ = output array for factors
 378  *
 379  * Returns:     A Lim-Lee prime, or null if generation failed.
 380  *
 381  * Use:         Generates Lim-Lee primes.  A Lim-Lee prime %$p$% is one which
 382  *              satisfies %$p = 2 \prod_i q_i + 1$%, where all of the %$q_i$%
 383  *              are large enough to resist square-root discrete log
 384  *              algorithms.
 385  *
 386  *              If we succeed, and @f@ is non-null, we write the array of
 387  *              factors chosen to @f@ for the benefit of the caller.
 388  */
 389
 390 mp *limlee(const char *name, mp *d, mp *newp,
 391            unsigned ql, unsigned pl, grand *r,
 392            unsigned on, pgen_proc *oev, void *oec,
 393            pgen_proc *iev, void *iec,
 394            size_t *nf, mp ***f)
 395 {
 396 #ifdef notdef
 397   dstr dn = DSTR_INIT;
 398   unsigned qql;
 399   mp *qq = 0;
 400   unsigned nn;
 401   unsigned mm;
 402   mp **v;
 403   octet *c;
 404   unsigned i;
 405   unsigned long seq = 0;
 406   pgen_event ev;
 407   unsigned ntest;
 408   rabin rb;
 409   pgen_filterctx pf;
 410
 411   /* --- First of all, decide on a number of factors to make --- */
 412
 413   nn = pl/ql;
 414   qql = pl%ql;
 415   if (!nn)
 416     return (0);
 417   else if (qql && nn > 1) {
 418     nn--;
 419     qql += ql;
 420   }
 421
 422   /* --- Now decide on how many primes I'll actually generate --- *
 423    *
 424    * The formula %$m = \max(3 n + 5, 25)$% comes from GPG's prime generation
 425    * library.
 426    */
 427
 428   mm = nn * 3 + 5;
 429   if (mm < 25)
 430     mm = 25;
 431
 432   /* --- Now allocate the working memory --- */
 433
 434   v = xmalloc(mm * sizeof(mp *));
 435   c = xmalloc(mm);
 436
 437   /* --- Initialize everything and try to find a prime --- */
 438
 439   ev.name = name;
 440   ev.m = 0;
 441   ev.steps = on;
 442   ev.tests = ntest = rabin_iters(pl);
 443   ev.r = r;
 444
 445   if (oev && oev(PGEN_BEGIN, &ev, oec) == PGEN_ABORT)
 446     goto fail;
 447
 448   pf.step = 2;
 449   if (qql) {
 450     dstr_putf(&dn, "%s*", name);
 451     qq = mprand(d, qql, r, 1);
 452     qq = pgen(dn.buf, qq, qq, iev, iec,
 453               0, pgen_filter, &pf, rabin_iters(qql), pgen_test, &rb);
 454   }
 455
 456 again:
 457   comb_init(c, mm, nn);
 458   for (i = 0; i < mm; i++)
 459     v[i] = 0;
 460
 461   /* --- The main combinations loop --- */
 462
 463   do {
 464     mpmul mmul = MPMUL_INIT;
 465
 466     /* --- Multiply a bunch of primes together --- */
 467
 468     if (qq) {
 469       mpmul_add(&mmul, qq);
 470     }
 471     for (i = 0; i < mm; i++) {
 472       if (!c[i])
 473         continue;
 474       if (!v[i]) {
 475         mp *z;
 476
 477         DRESET(&dn);
 478         dstr_putf(&dn, "%s_%lu] = ", name, seq++);
 479         z = mprand(newp, ql, ev.r, 1);
 480         z = pgen(dn.buf, z, z, iev, iec,
 481                  0, pgen_filter, &pf, rabin_iters(ql), pgen_test, &rb);
 482         v[i] = z;
 483       }
 484       mpmul_add(&mmul, v[i]);
 485     }
 486
 487     /* --- Now do some testing --- */
 488
 489     {
 490       mp *p = mpmul_done(&mmul);
 491       mp *g;
 492       int rc;
 493
 494       /* --- Check for small factors --- */
 495
 496       p = mp_lsl(p, p, 1);
 497       p = mp_add(p, p, MP_ONE);
 498       rc = pfilt_smallfactor(p);
 499       if (rc == PGEN_FAIL) {
 500         mp_drop(p);
 501         continue;
 502       }
 503
 504       /* --- Send an event out --- */
 505
 506       ev.m = p;
 507       if (oev && oev(PGEN_TRY, &ev, oec) == PGEN_ABORT) {
 508         mp_drop(p);
 509         goto fail;
 510       }
 511
 512       /* --- Do the Rabin testing --- */
 513
 514       rabin_create(&rb, p);
 515       g = MP_NEW;
 516       do {
 517         g = mprand_range(g, p, ev.r, 1);
 518         rc = rabin_test(&rb, g);
 519         if (rc == PGEN_PASS) {
 520           ev.tests--;
 521           if (!ev.tests)
 522             rc = PGEN_DONE;
 523         }
 524         if (oev &&oev(rc, &ev, oec) == PGEN_ABORT)
 525           rc = PGEN_ABORT;
 526       } while (rc == PGEN_PASS);
 527
 528       rabin_destroy(&rb);
 529       mp_drop(g);
 530       if (rc == PGEN_DONE)
 531         d = p;
 532       else
 533         mp_drop(p);
 534       if (rc == PGEN_ABORT)
 535         goto fail;
 536       if (rc == PGEN_DONE)
 537         goto done;
 538       ev.tests = ntest;
 539       ev.m = 0;
 540     }
 541   } while (comb_next(c, mm, nn));
 542
 543   /* --- That failed --- */
 544
 545   if (ev.steps) {
 546     ev.steps--;
 547     if (!ev.steps) {
 548       if (oev)
 549         oev(PGEN_ABORT, &ev, &oec);
 550       goto fail;
 551     }
 552   }
 553
 554   for (i = 0; i < mm; i++)
 555     mp_drop(v[i]);
 556   goto again;
 557
 558   /* --- We did it! --- */
 559
 560 done: {
 561     mp **vv = 0;
 562     if (f) {
 563       if (qq)
 564         nn++;
 565       *nf = nn;
 566       *f = vv = xmalloc(nn * sizeof(mp *));
 567     }
 568
 569     for (i = 0; i < mm; i++) {
 570       if (c[i] && vv)
 571         *vv++ = v[i];
 572       else if (v[i])
 573         mp_drop(v[i]);
 574     }
 575     if (qq) {
 576       if (vv)
 577         *vv++ = qq;
 578       else
 579         mp_drop(qq);
 580     }
 581     xfree(v);
 582     xfree(c);
 583     dstr_destroy(&dn);
 584     return (d);
 585   }
 586
 587   /* --- We blew it --- */
 588
 589 fail:
 590   for (i = 0; i < mm; i++)
 591     mp_drop(v[i]);
 592   if (qq)
 593     mp_drop(qq);
 594   xfree(v);
 595   xfree(c);
 596   dstr_destroy(&dn);
 597   return (0);
 598 #else
 599   limlee_stepctx l;
 600   rabin rr;
 601
 602   l.f = 0; if (f) l.f |= LIMLEE_KEEPFACTORS;
 603   l.newp = newp;
 604   l.pl = pl; l.ql = ql;
 605   l.pops = 0;
 606   l.iev = iev;
 607   l.iec = iec;
 608
 609   d = pgen(name, d, 0, oev, oec, on, limlee_step, &l,
 610            rabin_iters(pl), pgen_test, &rr);
 611
 612   if (f) {
 613     mp **v;
 614     size_t i;
 615     v = xmalloc(l.nf * sizeof(mp *));
 616     for (i = 0; i < l.nf; i++)
 617       v[i] = l.v[i].p;
 618     xfree(l.v);
 619     *f = v;
 620     *nf = l.nf;
 621   }
 622
 623   return (d);
 624 #endif
 625 }
 626
 627 /*----- That's all, folks -------------------------------------------------*/