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

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