chiark - git - mdw - catacomb/blob - limlee.c

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