chiark - git - mdw - catacomb/blob - math/strongprime.c

   1 /* -*-c-*-
   2  *
   3  * Generate `strong' prime numbers
   4  *
   5  * (c) 1999 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/dstr.h>
  31
  32 #include "grand.h"
  33 #include "mp.h"
  34 #include "mpmont.h"
  35 #include "mprand.h"
  36 #include "pgen.h"
  37 #include "pfilt.h"
  38 #include "rabin.h"
  39
  40 /*----- Main code ---------------------------------------------------------*/
  41
  42 /* --- @strongprime_setup@ --- *
  43  *
  44  * Arguments:   @const char *name@ = pointer to name root
  45  *              @mp *d@ = destination for search start point
  46  *              @pfilt *f@ = where to store filter jump context
  47  *              @unsigned nbits@ = number of bits wanted
  48  *              @grand *r@ = random number source
  49  *              @unsigned n@ = number of attempts to make
  50  *              @pgen_proc *event@ = event handler function
  51  *              @void *ectx@ = argument for the event handler
  52  *
  53  * Returns:     A starting point for a `strong' prime search, or zero.
  54  *
  55  * Use:         Sets up for a strong prime search, so that primes with
  56  *              particular properties can be found.  It's probably important
  57  *              to note that the number left in the filter context @f@ is
  58  *              congruent to 2 (mod 4); that the jump value is twice the
  59  *              product of two large primes; and that the starting point is
  60  *              at least %$3 \cdot 2^{N-2}$%.  (Hence, if you multiply two
  61  *              such numbers, the product is at least
  62  *
  63  *                      %$9 \cdot 2^{2N-4} > 2^{2N-1}$%
  64  *
  65  *              i.e., it will be (at least) a %$2 N$%-bit value.
  66  */
  67
  68 mp *strongprime_setup(const char *name, mp *d, pfilt *f, unsigned nbits,
  69                       grand *r, unsigned n, pgen_proc *event, void *ectx)
  70 {
  71   mp *s, *t, *q;
  72   dstr dn = DSTR_INIT;
  73   unsigned slop, nb, u, i;
  74
  75   mp *rr = d;
  76   pgen_filterctx c;
  77   pgen_jumpctx j;
  78   rabin rb;
  79
  80   /* --- Figure out how large the smaller primes should be --- *
  81    *
  82    * We want them to be `as large as possible', subject to the constraint
  83    * that we produce a number of the requested size at the end.  This is
  84    * tricky, because the final prime search is going to involve quite large
  85    * jumps from its starting point; the size of the jumps are basically
  86    * determined by our choice here, and if they're too big then we won't find
  87    * a prime in time.
  88    *
  89    * Let's suppose we're trying to make an %$N$%-bit prime.  The expected
  90    * number of steps tends to increase linearly with size, i.e., we need to
  91    * take about %2^k N$% steps for some %$k$%.  If we're jumping by a
  92    * %$J$%-bit quantity each time, from an %$N$%-bit starting point, then we
  93    * will only be able to find a match if %$2^k N 2^{J-1} \le 2^{N-1}$%,
  94    * i.e., if %$J \le N - (k + \log_2 N)$%.
  95    *
  96    * Experimentation shows that taking %$k + \log_2 N = 12$% works well for
  97    * %$N = 1024$%, so %$k = 2$%.  Add a few extra bits for luck.
  98    */
  99
 100   for (i = 1; i && nbits >> i; i <<= 1); assert(i);
 101   for (slop = 6, nb = nbits; nb > 1; i >>= 1) {
 102     u = nb >> i;
 103     if (u) { slop += i; nb = u; }
 104   }
 105   if (nbits/2 <= slop) return (0);
 106
 107   /* --- Choose two primes %$s$% and %$t$% of half the required size --- */
 108
 109   nb = nbits/2 - slop;
 110   c.step = 1;
 111
 112   rr = mprand(rr, nb, r, 1);
 113   DRESET(&dn); dstr_putf(&dn, "%s [s]", name);
 114   if ((s = pgen(dn.buf, MP_NEWSEC, rr, event, ectx, n, pgen_filter, &c,
 115                 rabin_iters(nb), pgen_test, &rb)) == 0)
 116     goto fail_s;
 117
 118   rr = mprand(rr, nb, r, 1);
 119   DRESET(&dn); dstr_putf(&dn, "%s [t]", name);
 120   if ((t = pgen(dn.buf, MP_NEWSEC, rr, event, ectx, n, pgen_filter, &c,
 121                 rabin_iters(nb), pgen_test, &rb)) == 0)
 122     goto fail_t;
 123
 124   /* --- Choose a suitable value for %$r = 2it + 1$% for some %$i$% --- */
 125
 126   rr = mp_lsl(rr, t, 1);
 127   pfilt_create(&c.f, rr);
 128   rr = mp_lsl(rr, rr, slop - 1);
 129   rr = mp_add(rr, rr, MP_ONE);
 130   DRESET(&dn); dstr_putf(&dn, "%s [r]", name);
 131   j.j = &c.f;
 132   q = pgen(dn.buf, MP_NEW, rr, event, ectx, n, pgen_jump, &j,
 133            rabin_iters(nb + slop), pgen_test, &rb);
 134   pfilt_destroy(&c.f);
 135   if (!q)
 136     goto fail_r;
 137
 138   /* --- Select a suitable congruence class for %$p$% --- *
 139    *
 140    * This computes %$p_0 = 2 s (s^{r - 2} \bmod r) - 1$%.
 141    */
 142
 143   {
 144     mpmont mm;
 145
 146     mpmont_create(&mm, q);
 147     rr = mp_sub(rr, q, MP_TWO);
 148     rr = mpmont_exp(&mm, rr, s, rr);
 149     mpmont_destroy(&mm);
 150     rr = mp_mul(rr, rr, s);
 151     rr = mp_lsl(rr, rr, 1);
 152     rr = mp_sub(rr, rr, MP_ONE);
 153   }
 154
 155   /* --- Pick a starting point for the search --- *
 156    *
 157    * Select %$3 \cdot 2^{N-2} < p_1 < 2^N$% at random, only with
 158    * %$p_1 \equiv p_0 \pmod{2 r s}$.
 159    */
 160
 161   {
 162     mp *x, *y;
 163     x = mp_mul(MP_NEW, q, s);
 164     x = mp_lsl(x, x, 1);
 165     pfilt_create(f, x); /* %$2 r s$% */
 166     y = mprand(MP_NEW, nbits, r, 0);
 167     y = mp_setbit(y, y, nbits - 2);
 168     rr = mp_leastcongruent(rr, y, rr, x);
 169     mp_drop(x); mp_drop(y);
 170   }
 171
 172   /* --- Return the result --- */
 173
 174   mp_drop(q);
 175   mp_drop(t);
 176   mp_drop(s);
 177   dstr_destroy(&dn);
 178   return (rr);
 179
 180   /* --- Tidy up if something failed --- */
 181
 182 fail_r:
 183   mp_drop(t);
 184 fail_t:
 185   mp_drop(s);
 186 fail_s:
 187   mp_drop(rr);
 188   dstr_destroy(&dn);
 189   return (0);
 190 }
 191
 192 /* --- @strongprime@ --- *
 193  *
 194  * Arguments:   @const char *name@ = pointer to name root
 195  *              @mp *d@ = destination integer
 196  *              @unsigned nbits@ = number of bits wanted
 197  *              @grand *r@ = random number source
 198  *              @unsigned n@ = number of attempts to make
 199  *              @pgen_proc *event@ = event handler function
 200  *              @void *ectx@ = argument for the event handler
 201  *
 202  * Returns:     A `strong' prime, or zero.
 203  *
 204  * Use:         Finds `strong' primes.  A strong prime %$p$% is such that
 205  *
 206  *                * %$p - 1$% has a large prime factor %$r$%,
 207  *                * %$p + 1$% has a large prime factor %$s$%, and
 208  *                * %$r - 1$% has a large prime factor %$t$%.
 209  */
 210
 211 mp *strongprime(const char *name, mp *d, unsigned nbits, grand *r,
 212                 unsigned n, pgen_proc *event, void *ectx)
 213 {
 214   mp *p;
 215   pfilt f;
 216   pgen_jumpctx j;
 217   rabin rb;
 218
 219   if (d) mp_copy(d);
 220   p = strongprime_setup(name, d, &f, nbits, r, n, event, ectx);
 221   if (!p) { mp_drop(d); return (0); }
 222   j.j = &f;
 223   p = pgen(name, p, p, event, ectx, n, pgen_jump, &j,
 224            rabin_iters(nbits), pgen_test, &rb);
 225   if (mp_bits(p) != nbits) { mp_drop(p); return (0); }
 226   pfilt_destroy(&f);
 227   mp_drop(d);
 228   return (p);
 229 }
 230
 231 /*----- That's all, folks -------------------------------------------------*/