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).
  59  */
  60
  61 mp *strongprime_setup(const char *name, mp *d, pfilt *f, unsigned nbits,
  62                       grand *r, unsigned n, pgen_proc *event, void *ectx)
  63 {
  64   mp *s, *t, *q;
  65   dstr dn = DSTR_INIT;
  66
  67   mp *rr = d;
  68   pgen_filterctx c;
  69   pgen_jumpctx j;
  70   rabin rb;
  71
  72   /* --- The bitslop parameter --- *
  73    *
  74    * There's quite a lot of prime searching to be done.  The constant
  75    * @BITSLOP@ is a (low) approximation to the base-2 log of the expected
  76    * number of steps to find a prime number.  Experimentation shows that
  77    * numbers around 10 seem to be good.
  78    */
  79
  80 #define BITSLOP 12
  81
  82   /* --- Choose two primes %$s$% and %$t$% of half the required size --- */
  83
  84   assert(((void)"nbits too small in strongprime_setup", nbits/2 > BITSLOP));
  85   nbits = nbits/2 - BITSLOP;
  86   c.step = 1;
  87
  88   rr = mprand(rr, nbits, r, 1);
  89   DRESET(&dn); dstr_putf(&dn, "%s [s]", name);
  90   if ((s = pgen(dn.buf, MP_NEWSEC, rr, event, ectx, n, pgen_filter, &c,
  91            rabin_iters(nbits), pgen_test, &rb)) == 0)
  92     goto fail_s;
  93
  94   rr = mprand(rr, nbits, r, 1);
  95   DRESET(&dn); dstr_putf(&dn, "%s [t]", name);
  96   if ((t = pgen(dn.buf, MP_NEWSEC, rr, event, ectx, n, pgen_filter, &c,
  97                 rabin_iters(nbits), pgen_test, &rb)) == 0)
  98     goto fail_t;
  99
 100   /* --- Choose a suitable value for %$r = 2it + 1$% for some %$i$% --- */
 101
 102   rr = mp_lsl(rr, t, 1);
 103   pfilt_create(&c.f, rr);
 104   rr = mp_lsl(rr, rr, BITSLOP - 1);
 105   rr = mp_add(rr, rr, MP_ONE);
 106   DRESET(&dn); dstr_putf(&dn, "%s [r]", name);
 107   j.j = &c.f;
 108   nbits += BITSLOP;
 109   q = pgen(dn.buf, MP_NEW, rr, event, ectx, n, pgen_jump, &j,
 110            rabin_iters(nbits), pgen_test, &rb);
 111   pfilt_destroy(&c.f);
 112   if (!q)
 113     goto fail_r;
 114
 115   /* --- Select a suitable starting-point for finding %$p$% --- *
 116    *
 117    * This computes %$p_0 = 2(s^{r - 2} \bmod r)s - 1$%.
 118    */
 119
 120   {
 121     mpmont mm;
 122
 123     mpmont_create(&mm, q);
 124     rr = mp_sub(rr, q, MP_TWO);
 125     rr = mpmont_exp(&mm, rr, s, rr);
 126     mpmont_destroy(&mm);
 127     rr = mp_mul(rr, rr, s);
 128     rr = mp_lsl(rr, rr, 1);
 129     rr = mp_sub(rr, rr, MP_ONE);
 130   }
 131
 132   /* --- Now find %$p = p_0 + 2jrs$% for some %$j$% --- */
 133
 134   {
 135     mp *x;
 136     x = mp_mul(MP_NEW, q, s);
 137     x = mp_lsl(x, x, 1);
 138     pfilt_create(f, x);
 139     x = mp_lsl(x, x, BITSLOP - 1);
 140     rr = mp_add(rr, rr, x);
 141     mp_drop(x);
 142   }
 143
 144   /* --- Return the result --- */
 145
 146   mp_drop(q);
 147   mp_drop(t);
 148   mp_drop(s);
 149   dstr_destroy(&dn);
 150   return (rr);
 151
 152   /* --- Tidy up if something failed --- */
 153
 154 fail_r:
 155   mp_drop(t);
 156 fail_t:
 157   mp_drop(s);
 158 fail_s:
 159   mp_drop(rr);
 160   dstr_destroy(&dn);
 161   return (0);
 162
 163 #undef BITSLOP
 164 }
 165
 166 /* --- @strongprime@ --- *
 167  *
 168  * Arguments:   @const char *name@ = pointer to name root
 169  *              @mp *d@ = destination integer
 170  *              @unsigned nbits@ = number of bits wanted
 171  *              @grand *r@ = random number source
 172  *              @unsigned n@ = number of attempts to make
 173  *              @pgen_proc *event@ = event handler function
 174  *              @void *ectx@ = argument for the event handler
 175  *
 176  * Returns:     A `strong' prime, or zero.
 177  *
 178  * Use:         Finds `strong' primes.  A strong prime %$p$% is such that
 179  *
 180  *                * %$p - 1$% has a large prime factor %$r$%,
 181  *                * %$p + 1$% has a large prime factor %$s$%, and
 182  *                * %$r - 1$% has a large prime factor %$t$%.
 183  *
 184  *              The numbers produced may be slightly larger than requested,
 185  *              by a few bits.
 186  */
 187
 188 mp *strongprime(const char *name, mp *d, unsigned nbits, grand *r,
 189                 unsigned n, pgen_proc *event, void *ectx)
 190 {
 191   pfilt f;
 192   pgen_jumpctx j;
 193   rabin rb;
 194
 195   d = strongprime_setup(name, d, &f, nbits, r, n, event, ectx);
 196   j.j = &f;
 197   d = pgen(name, d, d, event, ectx, n, pgen_jump, &j,
 198            rabin_iters(nbits), pgen_test, &rb);
 199   pfilt_destroy(&f);
 200   return (d);
 201 }
 202
 203 /*----- That's all, folks -------------------------------------------------*/