git.distorted.org.uk Git - u/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->d.buf, 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   unsigned qql;
 150
 151   /* --- First of all, decide on a number of factors to make --- */
 152
 153   l->nf = l->pl / l->ql;
 154   qql = l->pl % l->ql;
 155   if (!l->nf)
 156     return (PGEN_ABORT);
 157   else if (qql && l->nf > 1) {
 158     l->nf--;
 159     qql += l->ql;
 160   }
 161
 162   /* --- Now decide on how many primes I'll actually generate --- *
 163    *
 164    * The formula %$m = \max(3 n + 5, 25)$% comes from GPG's prime generation
 165    * library.
 166    */
 167
 168   l->poolsz = l->nf * 3 + 5;
 169   if (l->poolsz < 25)
 170     l->poolsz = 25;
 171
 172   /* --- Allocate and initialize the various tables --- */
 173
 174   l->c = xmalloc(l->poolsz);
 175   l->v = xmalloc(l->poolsz * sizeof(limlee_factor));
 176   comb_init(l->c, l->poolsz, l->nf);
 177   for (i = 0; i < l->poolsz; i++)
 178     l->v[i].p = 0;
 179
 180   /* --- Other bits of initialization --- */
 181
 182   l->seq = 0;
 183   dstr_create(&l->d);
 184   if (!l->pops) {
 185     l->pops = &primeops_simple;
 186     l->pc = 0;
 187   }
 188
 189   /* --- Find a big prime --- */
 190
 191   if (!qql)
 192     l->qq.p = 0;
 193   else {
 194     dstr_putf(&l->d, "%s*", ev->name);
 195     l->pops->pgen(&l->qq, qql, l);
 196   }
 197
 198   return (PGEN_TRY);
 199 }
 200
 201 /* --- @next@ --- *
 202  *
 203  * Arguments:   @int rq@ = request which triggered this call
 204  *              @pgen_event *ev@ = pointer to event block
 205  *              @limlee_stepctx *l@ = pointer to Lim-Lee context
 206  *
 207  * Returns:     A @PGEN@ result code.
 208  *
 209  * Use:         Initializes the stepper.
 210  */
 211
 212 static int next(int rq, pgen_event *ev, limlee_stepctx *l)
 213 {
 214   mp *p;
 215   int rc;
 216
 217   mp_drop(ev->m);
 218
 219   for (;;) {
 220     size_t i;
 221     mpmul mm = MPMUL_INIT;
 222
 223     /* --- Step on to next combination --- */
 224
 225     if (rq == PGEN_TRY && !comb_next(l->c, l->poolsz, l->nf)) {
 226       for (i = 0; i < l->poolsz; i++) {
 227         l->pops->pfree(&l->v[i], l);
 228         l->v[i].p = 0;
 229       }
 230     }
 231     rq = PGEN_TRY; /* For next time through */
 232
 233     /* --- Gather up some factors --- */
 234
 235     if (l->qq.p)
 236       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(&l->d);
 242         dstr_putf(&l->d, "%s_%lu", ev->name, l->seq++);
 243         l->pops->pgen(&l->v[i], l->ql, l);
 244       }
 245       mpmul_add(&mm, l->v[i].p);
 246     }
 247
 248     /* --- Check it for small factors --- */
 249
 250     p = mpmul_done(&mm);
 251     p = mp_lsl(p, p, 1);
 252     p->v[0] |= 1;
 253     if ((rc = pfilt_smallfactor(p)) != PGEN_FAIL)
 254       break;
 255     mp_drop(p);
 256   }
 257
 258   ev->m = p;
 259   return (rc);
 260 }
 261
 262 /* --- @done@ --- *
 263  *
 264  * Arguments:   @pgen_event *ev@ = pointer to event block
 265  *              @limlee_stepctx *l@ = pointer to Lim-Lee context
 266  *
 267  * Returns:     A @PGEN@ result code.
 268  *
 269  * Use:         Finalizes the stepper.  The output values in the context
 270  *              take on their final results; other resources are discarded.
 271  */
 272
 273 static int done(pgen_event *ev, limlee_stepctx *l)
 274 {
 275   size_t i, j;
 276   limlee_factor *v;
 277
 278   /* --- If an output vector of factors is wanted, produce one --- */
 279
 280   if (!(l->f & LIMLEE_KEEPFACTORS))
 281     v = 0;
 282   else {
 283     if (l->qq.p)
 284       l->nf++;
 285     v = xmalloc(l->nf * sizeof(limlee_factor));
 286   }
 287
 288   for (i = 0, j = 0; i < l->poolsz; i++) {
 289     if (v && l->c[i])
 290       v[j++] = l->v[i];
 291     else if (l->v[i].p)
 292       l->pops->pfree(&l->v[i], l);
 293   }
 294
 295   if (l->qq.p) {
 296     if (v)
 297       v[j++] = l->qq;
 298     else
 299       l->pops->pfree(&l->qq, l);
 300   }
 301
 302   xfree(l->v);
 303   l->v = v;
 304
 305   /* --- Free other resources --- */
 306
 307   xfree(l->c);
 308   dstr_destroy(&l->d);
 309
 310   /* --- Done --- */
 311
 312   return (PGEN_DONE);
 313 }
 314
 315 /* --- @limlee_step@ --- */
 316
 317 int limlee_step(int rq, pgen_event *ev, void *p)
 318 {
 319   limlee_stepctx *l = p;
 320   int rc;
 321
 322   switch (rq) {
 323     case PGEN_BEGIN:
 324       if ((rc = init(ev, l)) != PGEN_TRY)
 325         return (rc);
 326     case PGEN_TRY:
 327       return (next(rq, ev, l));
 328     case PGEN_DONE:
 329       return (done(ev, l));
 330   }
 331   return (PGEN_ABORT);
 332 }
 333
 334 /*----- Main code ---------------------------------------------------------*/
 335
 336 /* --- @limlee@ --- *
 337  *
 338  * Arguments:   @const char *name@ = pointer to name root
 339  *              @mp *d@ = pointer to destination integer
 340  *              @mp *newp@ = how to generate factor primes
 341  *              @unsigned ql@ = size of individual factors
 342  *              @unsigned pl@ = size of large prime
 343  *              @grand *r@ = a random number source
 344  *              @unsigned on@ = number of outer attempts to make
 345  *              @pgen_proc *oev@ = outer event handler function
 346  *              @void *oec@ = argument for the outer event handler
 347  *              @pgen_proc *iev@ = inner event handler function
 348  *              @void *iec@ = argument for the inner event handler
 349  *              @size_t *nf@, @mp ***f@ = output array for factors
 350  *
 351  * Returns:     A Lim-Lee prime, or null if generation failed.
 352  *
 353  * Use:         Generates Lim-Lee primes.  A Lim-Lee prime %$p$% is one which
 354  *              satisfies %$p = 2 \prod_i q_i + 1$%, where all of the %$q_i$%
 355  *              are large enough to resist square-root discrete log
 356  *              algorithms.
 357  *
 358  *              If we succeed, and @f@ is non-null, we write the array of
 359  *              factors chosen to @f@ for the benefit of the caller.
 360  */
 361
 362 mp *limlee(const char *name, mp *d, mp *newp,
 363            unsigned ql, unsigned pl, grand *r,
 364            unsigned on, pgen_proc *oev, void *oec,
 365            pgen_proc *iev, void *iec,
 366            size_t *nf, mp ***f)
 367 {
 368   limlee_stepctx l;
 369   rabin rr;
 370
 371   l.f = 0; if (f) l.f |= LIMLEE_KEEPFACTORS;
 372   l.newp = newp;
 373   l.pl = pl; l.ql = ql;
 374   l.pops = 0;
 375   l.iev = iev;
 376   l.iec = iec;
 377   l.r = r;
 378
 379   d = pgen(name, d, 0, oev, oec, on, limlee_step, &l,
 380            rabin_iters(pl), pgen_test, &rr);
 381
 382   if (d && f) {
 383     mp **v;
 384     size_t i;
 385     v = xmalloc(l.nf * sizeof(mp *));
 386     for (i = 0; i < l.nf; i++)
 387       v[i] = l.v[i].p;
 388     xfree(l.v);
 389     *f = v;
 390     *nf = l.nf;
 391   }
 392
 393   return (d);
 394 }
 395
 396 /*----- That's all, folks -------------------------------------------------*/