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

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