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   mp *p;
 115
 116   p = mprand(l->newp, pl, l->r, 1);
 117   pf.step = 2;
 118   p = pgen(l->u.s.name, p, p, l->iev, l->iec, 0, pgen_filter, &pf,
 119            PGEN_BAILLIEPSWNTESTS, pgen_bailliepswtest, 0);
 120   f->p = p;
 121 }
 122
 123 static void llfree(limlee_factor *f, limlee_stepctx *l)
 124 {
 125   mp_drop(f->p);
 126 }
 127
 128 static const limlee_primeops primeops_simple = { llgen, llfree };
 129
 130 /*----- Lim-Lee stepper ---------------------------------------------------*/
 131
 132 /* --- @init@ --- *
 133  *
 134  * Arguments:   @pgen_event *ev@ = pointer to event block
 135  *              @limlee_stepctx *l@ = pointer to Lim-Lee context
 136  *
 137  * Returns:     A @PGEN@ result code.
 138  *
 139  * Use:         Initializes the stepper.
 140  */
 141
 142 static int init(pgen_event *ev, limlee_stepctx *l)
 143 {
 144   size_t i;
 145
 146   /* --- First of all, decide on a number of factors to make --- */
 147
 148   l->nf = l->pl / l->ql;
 149   if (l->nf < 2) return (PGEN_ABORT);
 150   l->nf--;
 151
 152   /* --- Now decide on how many primes I'll actually generate --- *
 153    *
 154    * The formula %$m = \max(3 n + 5, 25)$% comes from GPG's prime generation
 155    * library.
 156    */
 157
 158   l->poolsz = l->nf * 3 + 5;
 159   if (l->poolsz < 25)
 160     l->poolsz = 25;
 161
 162   /* --- Allocate and initialize the various tables --- */
 163
 164   l->c = xmalloc(l->poolsz);
 165   l->v = xmalloc(l->poolsz * sizeof(limlee_factor));
 166   comb_init(l->c, l->poolsz, l->nf);
 167   for (i = 0; i < l->poolsz; i++)
 168     l->v[i].p = 0;
 169
 170   /* --- Other bits of initialization --- */
 171
 172   l->seq = 0;
 173   if (!l->pops) {
 174     l->pops = &primeops_simple;
 175     l->pc = 0;
 176   }
 177
 178   /* --- Find a big prime later --- */
 179
 180   l->qq.p = 0;
 181
 182   return (PGEN_TRY);
 183 }
 184
 185 /* --- @next@ --- *
 186  *
 187  * Arguments:   @int rq@ = request which triggered this call
 188  *              @pgen_event *ev@ = pointer to event block
 189  *              @limlee_stepctx *l@ = pointer to Lim-Lee context
 190  *
 191  * Returns:     A @PGEN@ result code.
 192  *
 193  * Use:         Initializes the stepper.
 194  */
 195
 196 static int next(int rq, pgen_event *ev, limlee_stepctx *l)
 197 {
 198   dstr d = DSTR_INIT;
 199   mp *p = 0;
 200   int rc;
 201   int dist;
 202   unsigned nb;
 203
 204   mp_drop(ev->m);
 205
 206   for (;;) {
 207     size_t i;
 208     mpmul mm = MPMUL_INIT;
 209
 210     /* --- Step on to next combination --- */
 211
 212     if (rq == PGEN_TRY && !comb_next(l->c, l->poolsz, l->nf)) {
 213       for (i = 0; i < l->poolsz; i++) {
 214         l->pops->pfree(&l->v[i], l);
 215         l->v[i].p = 0;
 216       }
 217     }
 218     rq = PGEN_TRY; /* For next time through */
 219
 220     /* --- If the large factor is performing badly, make a new one --- */
 221
 222     if (l->qq.p) {
 223       dist = l->u.s.disp < 0 ? -l->u.s.disp : l->u.s.disp;
 224       if (dist && dist > l->u.s.steps/3) {
 225         l->pops->pfree(&l->qq, l);
 226         l->qq.p = 0;
 227       }
 228     }
 229
 230     /* --- Gather up some factors --- */
 231
 232     if (l->qq.p) mpmul_add(&mm, l->qq.p);
 233     for (i = 0; i < l->poolsz; i++) {
 234       if (!l->c[i])
 235         continue;
 236       if (!l->v[i].p) {
 237         DRESET(&d);
 238         dstr_putf(&d, "%s_%lu", ev->name, l->seq++);
 239         l->u.s.name = d.buf;
 240         l->pops->pgen(&l->v[i], l->ql, l);
 241         if (!l->v[i].p)
 242           { mp_drop(mpmul_done(&mm)); rc = PGEN_ABORT; goto end; }
 243       }
 244       mpmul_add(&mm, l->v[i].p);
 245     }
 246
 247     /* --- Check on the large factor --- */
 248
 249     p = mpmul_done(&mm);
 250     if (!l->qq.p) {
 251       DRESET(&d);
 252       dstr_putf(&d, "%s*_%lu", ev->name, l->seq++);
 253       l->u.s.name = d.buf;
 254       l->pops->pgen(&l->qq, l->pl - mp_bits(p), l);
 255       if (!l->qq.p) { MP_DROP(p); p = 0; rc = PGEN_ABORT; break; }
 256       l->u.s.steps = l->u.s.disp = 0;
 257       p = mp_mul(p, p, l->qq.p);
 258     }
 259     p = mp_lsl(p, p, 1);
 260     p->v[0] |= 1;
 261
 262     nb = mp_bits(p);
 263     l->u.s.steps++;
 264     if (nb < l->pl) {
 265       l->u.s.disp--;
 266       continue;
 267     } else if (nb > l->pl) {
 268       l->u.s.disp++;
 269       continue;
 270     }
 271
 272     /* --- Check it for small factors --- */
 273
 274     if ((rc = pfilt_smallfactor(p)) != PGEN_FAIL)
 275       break;
 276     MP_DROP(p); p = 0;
 277   }
 278
 279 end:
 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   mp **v;
 393   size_t i;
 394
 395   l.f = 0; if (f) l.f |= LIMLEE_KEEPFACTORS;
 396   l.newp = newp;
 397   l.pl = pl; l.ql = ql;
 398   l.pops = 0;
 399   l.iev = iev;
 400   l.iec = iec;
 401   l.r = r;
 402
 403   d = pgen(name, d, 0, oev, oec, on, limlee_step, &l,
 404            PGEN_BAILLIEPSWNTESTS, pgen_bailliepswtest, &rr);
 405
 406   if (f) {
 407     if (!d) {
 408       for (i = 0; i < l.nf; i++)
 409         if (l.v[i].p) llfree(&l.v[i], &l);
 410     } else {
 411       v = xmalloc(l.nf * sizeof(mp *));
 412       for (i = 0; i < l.nf; i++)
 413         v[i] = l.v[i].p;
 414       xfree(l.v);
 415       *f = v;
 416       *nf = l.nf;
 417     }
 418   }
 419
 420   return (d);
 421 }
 422
 423 /*----- That's all, folks -------------------------------------------------*/