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

   1 /* -*-c-*-
   2  *
   3  * $Id: pfilt.c,v 1.4 2000/10/08 12:14:57 mdw Exp $
   4  *
   5  * Finding and testing prime numbers
   6  *
   7  * (c) 1999 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: pfilt.c,v $
  33  * Revision 1.4  2000/10/08 12:14:57  mdw
  34  * Remove vestiges of @primorial@.
  35  *
  36  * Revision 1.3  2000/08/15 21:44:27  mdw
  37  * (pfilt_smallfactor): New function for doing trial division the hard
  38  * way.
  39  *
  40  * (pfilt_create): Use @mpx_udivn@ for computing residues, for improved
  41  * performance.
  42  *
  43  * Pull the `small prime' test into a separate function, and do it
  44  * properly.
  45  *
  46  * Revision 1.2  2000/06/17 11:54:27  mdw
  47  * Use new MP memory management functions.
  48  *
  49  * Revision 1.1  1999/12/22 15:49:39  mdw
  50  * Renamed from `pgen'.  Reworking for new prime-search system.
  51  *
  52  * Revision 1.3  1999/12/10 23:28:35  mdw
  53  * Track suggested destination changes.
  54  *
  55  * Revision 1.2  1999/11/20 22:23:05  mdw
  56  * Add multiply-and-add function for Diffie-Hellman safe prime generation.
  57  *
  58  * Revision 1.1  1999/11/19 13:17:57  mdw
  59  * Prime number generator and tester.
  60  *
  61  */
  62
  63 /*----- Header files ------------------------------------------------------*/
  64
  65 #include "mp.h"
  66 #include "mpint.h"
  67 #include "pfilt.h"
  68 #include "pgen.h"
  69 #include "primetab.h"
  70
  71 /*----- Main code ---------------------------------------------------------*/
  72
  73 /* --- @smallenough@ --- *
  74  *
  75  * Arguments:   @mp *m@ = integer to test
  76  *
  77  * Returns:     One of the @PGEN@ result codes.
  78  *
  79  * Use:         Assuming that @m@ has been tested by trial division on every
  80  *              prime in the small-primes array, this function will return
  81  *              @PGEN_DONE@ if the number is less than the square of the
  82  *              largest small prime.
  83  */
  84
  85 static int smallenough(mp *m)
  86 {
  87   static mp *max = 0;
  88   int rc = PGEN_TRY;
  89
  90   if (!max) {
  91     max = mp_fromuint(MP_NEW, MAXPRIME);
  92     max = mp_sqr(max, max);
  93     max->a->n--; /* Permanent allocation */
  94   }
  95   if (MP_CMP(m, <, max))
  96     rc = PGEN_DONE;
  97   return (rc);
  98 }
  99
 100 /* --- @pfilt_smallfactor@ --- *
 101  *
 102  * Arguments:   @mp *m@ = integer to test
 103  *
 104  * Returns:     One of the @PGEN@ result codes.
 105  *
 106  * Use:         Tests a number by dividing by a number of small primes.  This
 107  *              is a useful first step if you're testing random primes; for
 108  *              sequential searches, @pfilt_create@ works better.
 109  */
 110
 111 int pfilt_smallfactor(mp *m)
 112 {
 113   int rc = PGEN_TRY;
 114   int i;
 115   size_t sz = MP_LEN(m);
 116   mpw *v = mpalloc(m->a, sz);
 117
 118   /* --- Fill in the residues --- */
 119
 120   for (i = 0; i < NPRIME; i++) {
 121     if (!mpx_udivn(v, v + sz, m->v, m->vl, primetab[i])) {
 122       if (MP_LEN(m) == 1 && m->v[0] == primetab[i])
 123         rc = PGEN_DONE;
 124       else
 125         rc = PGEN_FAIL;
 126     }
 127   }
 128
 129   /* --- Check for small primes --- */
 130
 131   if (rc == PGEN_TRY)
 132     rc = smallenough(m);
 133
 134   /* --- Done --- */
 135
 136   mpfree(m->a, v);
 137   return (rc);
 138 }
 139
 140 /* --- @pfilt_create@ --- *
 141  *
 142  * Arguments:   @pfilt *p@ = pointer to prime filtering context
 143  *              @mp *m@ = pointer to initial number to test
 144  *
 145  * Returns:     One of the @PGEN@ result codes.
 146  *
 147  * Use:         Tests an initial number for primality by computing its
 148  *              residue modulo various small prime numbers.  This is fairly
 149  *              quick, but not particularly certain.  If a @PGEN_TRY@
 150  *              result is returned, perform Rabin-Miller tests to confirm.
 151  */
 152
 153 int pfilt_create(pfilt *p, mp *m)
 154 {
 155   int rc = PGEN_TRY;
 156   int i;
 157   size_t sz = MP_LEN(m);
 158   mpw *v = mpalloc(m->a, sz);
 159
 160   /* --- Take a copy of the number --- */
 161
 162   mp_shrink(m);
 163   p->m = MP_COPY(m);
 164
 165   /* --- Fill in the residues --- */
 166
 167   for (i = 0; i < NPRIME; i++) {
 168     p->r[i] = mpx_udivn(v, v + sz, m->v, m->vl, primetab[i]);
 169     if (!p->r[i] && rc == PGEN_TRY) {
 170       if (MP_LEN(m) == 1 && m->v[0] == primetab[i])
 171         rc = PGEN_DONE;
 172       else
 173         rc = PGEN_FAIL;
 174     }
 175   }
 176
 177   /* --- Check for small primes --- */
 178
 179   if (rc == PGEN_TRY)
 180     rc = smallenough(m);
 181
 182   /* --- Done --- */
 183
 184   mpfree(m->a, v);
 185   return (rc);
 186 }
 187
 188 /* --- @pfilt_destroy@ --- *
 189  *
 190  * Arguments:   @pfilt *p@ = pointer to prime filtering context
 191  *
 192  * Returns:     ---
 193  *
 194  * Use:         Discards a context and all the resources it holds.
 195  */
 196
 197 void pfilt_destroy(pfilt *p)
 198 {
 199   mp_drop(p->m);
 200 }
 201
 202 /* --- @pfilt_step@ --- *
 203  *
 204  * Arguments:   @pfilt *p@ = pointer to prime filtering context
 205  *              @mpw step@ = how much to step the number
 206  *
 207  * Returns:     One of the @PGEN@ result codes.
 208  *
 209  * Use:         Steps a number by a small amount.  Stepping is much faster
 210  *              than initializing with a new number.  The test performed is
 211  *              the same simple one used by @primetab_create@, so @PGEN_TRY@
 212  *              results should be followed up by a Rabin-Miller test.
 213  */
 214
 215 int pfilt_step(pfilt *p, mpw step)
 216 {
 217   int rc = PGEN_TRY;
 218   int i;
 219
 220   /* --- Add the step on to the number --- */
 221
 222   p->m = mp_split(p->m);
 223   mp_ensure(p->m, MP_LEN(p->m) + 1);
 224   mpx_uaddn(p->m->v, p->m->vl, step);
 225   mp_shrink(p->m);
 226
 227   /* --- Update the residue table --- */
 228
 229   for (i = 0; i < NPRIME; i++) {
 230     p->r[i] = (p->r[i] + step) % primetab[i];
 231     if (!p->r[i] && rc == PGEN_TRY) {
 232       if (MP_LEN(p->m) == 1 && p->m->v[0] == primetab[i])
 233         rc = PGEN_DONE;
 234       else
 235         rc = PGEN_FAIL;
 236     }
 237   }
 238
 239   /* --- Check for small primes --- */
 240
 241   if (rc == PGEN_TRY)
 242     rc = smallenough(p->m);
 243
 244   /* --- Done --- */
 245
 246   return (rc);
 247 }
 248
 249 /* --- @pfilt_muladd@ --- *
 250  *
 251  * Arguments:   @pfilt *p@ = destination prime filtering context
 252  *              @const pfilt *q@ = source prime filtering context
 253  *              @mpw m@ = number to multiply by
 254  *              @mpw a@ = number to add
 255  *
 256  * Returns:     One of the @PGEN@ result codes.
 257  *
 258  * Use:         Multiplies the number in a prime filtering context by a
 259  *              small value and then adds a small value.  The destination
 260  *              should either be uninitialized or the same as the source.
 261  *
 262  *              Common things to do include multiplying by 2 and adding 0 to
 263  *              turn a prime into a jump for finding other primes with @q@ as
 264  *              a factor of @p - 1@, or multiplying by 2 and adding 1.
 265  */
 266
 267 int pfilt_muladd(pfilt *p, const pfilt *q, mpw m, mpw a)
 268 {
 269   int rc = PGEN_TRY;
 270   int i;
 271
 272   /* --- Multiply the big number --- */
 273
 274   {
 275     mp *d = mp_new(MP_LEN(q->m) + 2, q->m->f);
 276     mpx_umuln(d->v, d->vl, q->m->v, q->m->vl, m);
 277     mpx_uaddn(d->v, d->vl, a);
 278     if (p == q)
 279       mp_drop(p->m);
 280     mp_shrink(d);
 281     p->m = d;
 282   }
 283
 284   /* --- Gallivant through the residue table --- */
 285
 286   for (i = 0; i < NPRIME; i++) {
 287     p->r[i] = (q->r[i] * m + a) % primetab[i];
 288     if (!p->r[i] && rc == PGEN_TRY) {
 289       if (MP_LEN(p->m) == 1 && p->m->v[0] == primetab[i])
 290         rc = PGEN_DONE;
 291       else
 292         rc = PGEN_FAIL;
 293     }
 294   }
 295
 296   /* --- Check for small primes --- */
 297
 298   if (rc == PGEN_TRY)
 299     rc = smallenough(p->m);
 300
 301   /* --- Finished --- */
 302
 303   return (rc);
 304 }
 305
 306 /* --- @pfilt_jump@ --- *
 307  *
 308  * Arguments:   @pfilt *p@ = pointer to prime filtering context
 309  *              @const pfilt *j@ = pointer to another filtering context
 310  *
 311  * Returns:     One of the @PGEN@ result codes.
 312  *
 313  * Use:         Steps a number by a large amount.  Even so, jumping is much
 314  *              faster than initializing a new number.  The test peformed is
 315  *              the same simple one used by @primetab_create@, so @PGEN_TRY@
 316  *              results should be followed up by a Rabin-Miller test.
 317  *
 318  *              Note that the number stored in the @j@ context is probably
 319  *              better off being even than prime.  The important thing is
 320  *              that all of the residues for the number have already been
 321  *              computed.
 322  */
 323
 324 int pfilt_jump(pfilt *p, const pfilt *j)
 325 {
 326   int rc = PGEN_TRY;
 327   int i;
 328
 329   /* --- Add the step on --- */
 330
 331   p->m = mp_add(p->m, p->m, j->m);
 332
 333   /* --- Update the residue table --- */
 334
 335   for (i = 0; i < NPRIME; i++) {
 336     p->r[i] = p->r[i] + j->r[i];
 337     if (p->r[i] > primetab[i])
 338       p->r[i] -= primetab[i];
 339     if (!p->r[i] && rc == PGEN_TRY) {
 340       if (MP_LEN(p->m) == 1 && p->m->v[0] == primetab[i])
 341         rc = PGEN_DONE;
 342       else
 343         rc = PGEN_FAIL;
 344     }
 345   }
 346
 347   /* --- Check for small primes --- */
 348
 349   if (rc == PGEN_TRY)
 350     rc = smallenough(p->m);
 351
 352   /* --- Done --- */
 353
 354   return (rc);
 355 }
 356
 357 /*----- That's all, folks -------------------------------------------------*/