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

   1 /* -*-c-*-
   2  *
   3  * $Id: limlee.c,v 1.1 2000/07/09 21:30:58 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.1  2000/07/09 21:30:58  mdw
  34  * Lim-Lee prime generation.
  35  *
  36  */
  37
  38 /*----- Header files ------------------------------------------------------*/
  39
  40 #include <mLib/alloc.h>
  41 #include <mLib/dstr.h>
  42
  43 #include "limlee.h"
  44 #include "mpmul.h"
  45 #include "mprand.h"
  46 #include "pgen.h"
  47 #include "primorial.h"
  48 #include "rabin.h"
  49
  50 /*----- Main code ---------------------------------------------------------*/
  51
  52 /* --- @limlee@ --- *
  53  *
  54  * Arguments:   @const char *name@ = pointer to name root
  55  *              @mp *d@ = pointer to destination integer
  56  *              @mp *newp@ = how to generate factor primes
  57  *              @unsigned ql@ = size of individual factors
  58  *              @unsigned pl@ = size of large prime
  59  *              @grand *r@ = a random number source
  60  *              @unsigned on@ = number of outer attempts to make
  61  *              @pgen_proc *oev@ = outer event handler function
  62  *              @void *oec@ = argument for the outer event handler
  63  *              @pgen_proc *iev@ = inner event handler function
  64  *              @void *iec@ = argument for the inner event handler
  65  *              @size_t *nf@, @mp ***f@ = output array for factors
  66  *
  67  * Returns:     A Lim-Lee prime, or null if generation failed.
  68  *
  69  * Use:         Generates Lim-Lee primes.  A Lim-Lee prime %$p$% is one which
  70  *              satisfies %$p = 2 \prod_i q_i + 1$%, where all of the %$q_i$%
  71  *              are large enough to resist square-root discrete log
  72  *              algorithms.
  73  *
  74  *              If we succeed, and @f@ is non-null, we write the array of
  75  *              factors chosen to @f@ for the benefit of the caller.
  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 static int comb_next(octet *c, unsigned n, unsigned r)
  85 {
  86   unsigned g = 0;
  87
  88   /* --- How the algorithm works --- *
  89    *
  90    * Set bits start at the end and work their way towards the start.
  91    * Excepting bits already at the start, we scan for the lowest set bit, and
  92    * move it one place nearer the start.  A group of bits at the start are
  93    * counted and reset just below the `moved' bit.  If there is no moved bit
  94    * then we're done.
  95    */
  96
  97   /* --- Count the group at the start --- */
  98
  99   for (; *c; c++) {
 100     g++;
 101     *c = 0;
 102   }
 103   if (g == r)
 104     return (0);
 105
 106   /* --- Move the next bit down one --- *
 107    *
 108    * There must be one, because otherwise we'd have counted %$r$% bits
 109    * earlier.
 110    */
 111
 112   for (; !*c; c++)
 113     ;
 114   *c = 0;
 115   g++;
 116   for (; g; g--)
 117     *--c = 1;
 118   return (1);
 119 }
 120
 121 mp *limlee(const char *name, mp *d, mp *newp,
 122            unsigned ql, unsigned pl, grand *r,
 123            unsigned on, pgen_proc *oev, void *oec,
 124            pgen_proc *iev, void *iec,
 125            size_t *nf, mp ***f)
 126 {
 127   dstr dn = DSTR_INIT;
 128   unsigned qql;
 129   mp *qq = 0;
 130   unsigned nn;
 131   unsigned mm;
 132   mp **v;
 133   octet *c;
 134   unsigned i;
 135   unsigned long seq = 0;
 136   pgen_event ev;
 137   unsigned ntest;
 138   rabin rb;
 139   pgen_filterctx pf;
 140
 141   /* --- First of all, decide on a number of factors to make --- */
 142
 143   nn = pl/ql;
 144   qql = pl%ql;
 145   if (!nn)
 146     return (0);
 147   else if (qql && nn > 1) {
 148     nn--;
 149     qql += ql;
 150   }
 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   mm = nn * 3 + 5;
 159   if (mm < 25)
 160     mm = 25;
 161
 162   /* --- Now allocate the working memory --- */
 163
 164   primorial_setup();
 165   v = xmalloc(mm * sizeof(mp *));
 166   c = xmalloc(mm);
 167
 168   /* --- Initialize everything and try to find a prime --- */
 169
 170   ev.name = name;
 171   ev.m = 0;
 172   ev.steps = on;
 173   ev.tests = ntest = rabin_iters(pl);
 174   ev.r = r;
 175
 176   if (oev && oev(PGEN_BEGIN, &ev, oec) == PGEN_ABORT)
 177     goto fail;
 178
 179   if (qql) {
 180     dstr_putf(&dn, "%s [+]", name);
 181     qq = mprand(d, qql, r, 1);
 182     pf.step = 2;
 183     qq = pgen(dn.buf, qq, qq, iev, iec,
 184               0, pgen_filter, &pf, rabin_iters(qql), pgen_test, &rb);
 185   }
 186
 187 again:
 188   comb_init(c, mm, nn);
 189   for (i = 0; i < mm; i++)
 190     v[i] = 0;
 191
 192   /* --- The main combinations loop --- */
 193
 194   do {
 195     mpmul mmul = MPMUL_INIT;
 196
 197     /* --- Multiply a bunch of primes together --- */
 198
 199     if (qq) {
 200       mpmul_add(&mmul, qq);
 201     }
 202     for (i = 0; i < mm; i++) {
 203       if (!c[i])
 204         continue;
 205       if (!v[i]) {
 206         mp *z;
 207
 208         DRESET(&dn);
 209         dstr_putf(&dn, "%s [%lu] = ", name, seq++);
 210         z = mprand(newp, ql, ev.r, 1);
 211         z = pgen(dn.buf, z, z, iev, iec,
 212                  0, pgen_filter, &pf, rabin_iters(ql), pgen_test, &rb);
 213         v[i] = z;
 214       }
 215       mpmul_add(&mmul, v[i]);
 216     }
 217
 218     /* --- Now do some testing --- */
 219
 220     {
 221       mp *p = mpmul_done(&mmul);
 222       mp *g = newp;
 223       int rc;
 224
 225       /* --- Check for small factors --- */
 226
 227       p = mp_lsl(p, p, 1);
 228       p = mp_add(p, p, MP_ONE);
 229       mp_gcd(&g, 0, 0, p, primorial);
 230       if (MP_CMP(g, !=, MP_ONE)) {
 231         mp_drop(g);
 232         mp_drop(p);
 233         continue;
 234       }
 235       mp_drop(g);
 236
 237       /* --- Send an event out --- */
 238
 239       ev.m = p;
 240       if (oev && oev(PGEN_TRY, &ev, oec) == PGEN_ABORT) {
 241         mp_drop(p);
 242         goto fail;
 243       }
 244
 245       /* --- Do the Rabin testing --- */
 246
 247       rabin_create(&rb, p);
 248       g = MP_NEW;
 249       do {
 250         g = mprand_range(g, p, ev.r, 1);
 251         rc = rabin_test(&rb, g);
 252         if (rc == PGEN_PASS) {
 253           ev.tests--;
 254           if (!ev.tests)
 255             rc = PGEN_DONE;
 256         }
 257         if (oev &&oev(rc, &ev, oec) == PGEN_ABORT)
 258           rc = PGEN_ABORT;
 259       } while (rc == PGEN_PASS);
 260
 261       rabin_destroy(&rb);
 262       mp_drop(g);
 263       if (rc == PGEN_DONE)
 264         d = p;
 265       else
 266         mp_drop(p);
 267       if (rc == PGEN_ABORT)
 268         goto fail;
 269       if (rc == PGEN_DONE)
 270         goto done;
 271       ev.tests = ntest;
 272       ev.m = 0;
 273     }
 274   } while (comb_next(c, mm, nn));
 275
 276   /* --- That failed --- */
 277
 278   if (ev.steps) {
 279     ev.steps--;
 280     if (!ev.steps) {
 281       if (oev)
 282         oev(PGEN_ABORT, &ev, &oec);
 283       goto fail;
 284     }
 285   }
 286
 287   for (i = 0; i < mm; i++)
 288     mp_drop(v[i]);
 289   goto again;
 290
 291   /* --- We did it! --- */
 292
 293 done: {
 294     mp **vv = 0;
 295     if (f) {
 296       if (qq)
 297         nn++;
 298       *nf = nn;
 299       *f = vv = xmalloc(nn * sizeof(mp *));
 300     }
 301
 302     for (i = 0; i < mm; i++) {
 303       if (c[i] && vv)
 304         *vv++ = v[i];
 305       else if (v[i])
 306         mp_drop(v[i]);
 307     }
 308     if (qq) {
 309       if (vv)
 310         *vv++ = qq;
 311       else
 312         mp_drop(qq);
 313     }
 314     xfree(v);
 315     xfree(c);
 316     dstr_destroy(&dn);
 317     return (d);
 318   }
 319
 320   /* --- We blew it --- */
 321
 322 fail:
 323   for (i = 0; i < mm; i++)
 324     mp_drop(v[i]);
 325   if (qq)
 326     mp_drop(qq);
 327   xfree(v);
 328   xfree(c);
 329   dstr_destroy(&dn);
 330   return (0);
 331 }
 332