Back to index

lightning-sunbird  0.9+nobinonly
mpprime.c
Go to the documentation of this file.
00001 /*
00002  *  mpprime.c
00003  *
00004  *  Utilities for finding and working with prime and pseudo-prime
00005  *  integers
00006  *
00007  * ***** BEGIN LICENSE BLOCK *****
00008  * Version: MPL 1.1/GPL 2.0/LGPL 2.1
00009  *
00010  * The contents of this file are subject to the Mozilla Public License Version
00011  * 1.1 (the "License"); you may not use this file except in compliance with
00012  * the License. You may obtain a copy of the License at
00013  * http://www.mozilla.org/MPL/
00014  *
00015  * Software distributed under the License is distributed on an "AS IS" basis,
00016  * WITHOUT WARRANTY OF ANY KIND, either express or implied. See the License
00017  * for the specific language governing rights and limitations under the
00018  * License.
00019  *
00020  * The Original Code is the MPI Arbitrary Precision Integer Arithmetic library.
00021  *
00022  * The Initial Developer of the Original Code is
00023  * Michael J. Fromberger.
00024  * Portions created by the Initial Developer are Copyright (C) 1997
00025  * the Initial Developer. All Rights Reserved.
00026  *
00027  * Contributor(s):
00028  *   Netscape Communications Corporation
00029  *
00030  * Alternatively, the contents of this file may be used under the terms of
00031  * either the GNU General Public License Version 2 or later (the "GPL"), or
00032  * the GNU Lesser General Public License Version 2.1 or later (the "LGPL"),
00033  * in which case the provisions of the GPL or the LGPL are applicable instead
00034  * of those above. If you wish to allow use of your version of this file only
00035  * under the terms of either the GPL or the LGPL, and not to allow others to
00036  * use your version of this file under the terms of the MPL, indicate your
00037  * decision by deleting the provisions above and replace them with the notice
00038  * and other provisions required by the GPL or the LGPL. If you do not delete
00039  * the provisions above, a recipient may use your version of this file under
00040  * the terms of any one of the MPL, the GPL or the LGPL.
00041  *
00042  * ***** END LICENSE BLOCK ***** */
00043 
00044 #include "mpi-priv.h"
00045 #include "mpprime.h"
00046 #include "mplogic.h"
00047 #include <stdlib.h>
00048 #include <string.h>
00049 
00050 #define SMALL_TABLE 0 /* determines size of hard-wired prime table */
00051 
00052 #define RANDOM() rand()
00053 
00054 #include "primes.c"  /* pull in the prime digit table */
00055 
00056 /* 
00057    Test if any of a given vector of digits divides a.  If not, MP_NO
00058    is returned; otherwise, MP_YES is returned and 'which' is set to
00059    the index of the integer in the vector which divided a.
00060  */
00061 mp_err    s_mpp_divp(mp_int *a, const mp_digit *vec, int size, int *which);
00062 
00063 /* {{{ mpp_divis(a, b) */
00064 
00065 /*
00066   mpp_divis(a, b)
00067 
00068   Returns MP_YES if a is divisible by b, or MP_NO if it is not.
00069  */
00070 
00071 mp_err  mpp_divis(mp_int *a, mp_int *b)
00072 {
00073   mp_err  res;
00074   mp_int  rem;
00075 
00076   if((res = mp_init(&rem)) != MP_OKAY)
00077     return res;
00078 
00079   if((res = mp_mod(a, b, &rem)) != MP_OKAY)
00080     goto CLEANUP;
00081 
00082   if(mp_cmp_z(&rem) == 0)
00083     res = MP_YES;
00084   else
00085     res = MP_NO;
00086 
00087 CLEANUP:
00088   mp_clear(&rem);
00089   return res;
00090 
00091 } /* end mpp_divis() */
00092 
00093 /* }}} */
00094 
00095 /* {{{ mpp_divis_d(a, d) */
00096 
00097 /*
00098   mpp_divis_d(a, d)
00099 
00100   Return MP_YES if a is divisible by d, or MP_NO if it is not.
00101  */
00102 
00103 mp_err  mpp_divis_d(mp_int *a, mp_digit d)
00104 {
00105   mp_err     res;
00106   mp_digit   rem;
00107 
00108   ARGCHK(a != NULL, MP_BADARG);
00109 
00110   if(d == 0)
00111     return MP_NO;
00112 
00113   if((res = mp_mod_d(a, d, &rem)) != MP_OKAY)
00114     return res;
00115 
00116   if(rem == 0)
00117     return MP_YES;
00118   else
00119     return MP_NO;
00120 
00121 } /* end mpp_divis_d() */
00122 
00123 /* }}} */
00124 
00125 /* {{{ mpp_random(a) */
00126 
00127 /*
00128   mpp_random(a)
00129 
00130   Assigns a random value to a.  This value is generated using the
00131   standard C library's rand() function, so it should not be used for
00132   cryptographic purposes, but it should be fine for primality testing,
00133   since all we really care about there is good statistical properties.
00134 
00135   As many digits as a currently has are filled with random digits.
00136  */
00137 
00138 mp_err  mpp_random(mp_int *a)
00139 
00140 {
00141   mp_digit  next = 0;
00142   unsigned int       ix, jx;
00143 
00144   ARGCHK(a != NULL, MP_BADARG);
00145 
00146   for(ix = 0; ix < USED(a); ix++) {
00147     for(jx = 0; jx < sizeof(mp_digit); jx++) {
00148       next = (next << CHAR_BIT) | (RANDOM() & UCHAR_MAX);
00149     }
00150     DIGIT(a, ix) = next;
00151   }
00152 
00153   return MP_OKAY;
00154 
00155 } /* end mpp_random() */
00156 
00157 /* }}} */
00158 
00159 /* {{{ mpp_random_size(a, prec) */
00160 
00161 mp_err  mpp_random_size(mp_int *a, mp_size prec)
00162 {
00163   mp_err   res;
00164 
00165   ARGCHK(a != NULL && prec > 0, MP_BADARG);
00166   
00167   if((res = s_mp_pad(a, prec)) != MP_OKAY)
00168     return res;
00169 
00170   return mpp_random(a);
00171 
00172 } /* end mpp_random_size() */
00173 
00174 /* }}} */
00175 
00176 /* {{{ mpp_divis_vector(a, vec, size, which) */
00177 
00178 /*
00179   mpp_divis_vector(a, vec, size, which)
00180 
00181   Determines if a is divisible by any of the 'size' digits in vec.
00182   Returns MP_YES and sets 'which' to the index of the offending digit,
00183   if it is; returns MP_NO if it is not.
00184  */
00185 
00186 mp_err  mpp_divis_vector(mp_int *a, const mp_digit *vec, int size, int *which)
00187 {
00188   ARGCHK(a != NULL && vec != NULL && size > 0, MP_BADARG);
00189   
00190   return s_mpp_divp(a, vec, size, which);
00191 
00192 } /* end mpp_divis_vector() */
00193 
00194 /* }}} */
00195 
00196 /* {{{ mpp_divis_primes(a, np) */
00197 
00198 /*
00199   mpp_divis_primes(a, np)
00200 
00201   Test whether a is divisible by any of the first 'np' primes.  If it
00202   is, returns MP_YES and sets *np to the value of the digit that did
00203   it.  If not, returns MP_NO.
00204  */
00205 mp_err  mpp_divis_primes(mp_int *a, mp_digit *np)
00206 {
00207   int     size, which;
00208   mp_err  res;
00209 
00210   ARGCHK(a != NULL && np != NULL, MP_BADARG);
00211 
00212   size = (int)*np;
00213   if(size > prime_tab_size)
00214     size = prime_tab_size;
00215 
00216   res = mpp_divis_vector(a, prime_tab, size, &which);
00217   if(res == MP_YES) 
00218     *np = prime_tab[which];
00219 
00220   return res;
00221 
00222 } /* end mpp_divis_primes() */
00223 
00224 /* }}} */
00225 
00226 /* {{{ mpp_fermat(a, w) */
00227 
00228 /*
00229   Using w as a witness, try pseudo-primality testing based on Fermat's
00230   little theorem.  If a is prime, and (w, a) = 1, then w^a == w (mod
00231   a).  So, we compute z = w^a (mod a) and compare z to w; if they are
00232   equal, the test passes and we return MP_YES.  Otherwise, we return
00233   MP_NO.
00234  */
00235 mp_err  mpp_fermat(mp_int *a, mp_digit w)
00236 {
00237   mp_int  base, test;
00238   mp_err  res;
00239   
00240   if((res = mp_init(&base)) != MP_OKAY)
00241     return res;
00242 
00243   mp_set(&base, w);
00244 
00245   if((res = mp_init(&test)) != MP_OKAY)
00246     goto TEST;
00247 
00248   /* Compute test = base^a (mod a) */
00249   if((res = mp_exptmod(&base, a, a, &test)) != MP_OKAY)
00250     goto CLEANUP;
00251 
00252   
00253   if(mp_cmp(&base, &test) == 0)
00254     res = MP_YES;
00255   else
00256     res = MP_NO;
00257 
00258  CLEANUP:
00259   mp_clear(&test);
00260  TEST:
00261   mp_clear(&base);
00262 
00263   return res;
00264 
00265 } /* end mpp_fermat() */
00266 
00267 /* }}} */
00268 
00269 /*
00270   Perform the fermat test on each of the primes in a list until
00271   a) one of them shows a is not prime, or 
00272   b) the list is exhausted.
00273   Returns:  MP_YES if it passes tests.
00274            MP_NO  if fermat test reveals it is composite
00275            Some MP error code if some other error occurs.
00276  */
00277 mp_err mpp_fermat_list(mp_int *a, const mp_digit *primes, mp_size nPrimes)
00278 {
00279   mp_err rv = MP_YES;
00280 
00281   while (nPrimes-- > 0 && rv == MP_YES) {
00282     rv = mpp_fermat(a, *primes++);
00283   }
00284   return rv;
00285 }
00286 
00287 /* {{{ mpp_pprime(a, nt) */
00288 
00289 /*
00290   mpp_pprime(a, nt)
00291 
00292   Performs nt iteration of the Miller-Rabin probabilistic primality
00293   test on a.  Returns MP_YES if the tests pass, MP_NO if one fails.
00294   If MP_NO is returned, the number is definitely composite.  If MP_YES
00295   is returned, it is probably prime (but that is not guaranteed).
00296  */
00297 
00298 mp_err  mpp_pprime(mp_int *a, int nt)
00299 {
00300   mp_err   res;
00301   mp_int   x, amo, m, z;    /* "amo" = "a minus one" */
00302   int      iter;
00303   unsigned int jx;
00304   mp_size  b;
00305 
00306   ARGCHK(a != NULL, MP_BADARG);
00307 
00308   MP_DIGITS(&x) = 0;
00309   MP_DIGITS(&amo) = 0;
00310   MP_DIGITS(&m) = 0;
00311   MP_DIGITS(&z) = 0;
00312 
00313   /* Initialize temporaries... */
00314   MP_CHECKOK( mp_init(&amo));
00315   /* Compute amo = a - 1 for what follows...    */
00316   MP_CHECKOK( mp_sub_d(a, 1, &amo) );
00317 
00318   b = mp_trailing_zeros(&amo);
00319   if (!b) { /* a was even ? */
00320     res = MP_NO;
00321     goto CLEANUP;
00322   }
00323 
00324   MP_CHECKOK( mp_init_size(&x, MP_USED(a)) );
00325   MP_CHECKOK( mp_init(&z) );
00326   MP_CHECKOK( mp_init(&m) );
00327   MP_CHECKOK( mp_div_2d(&amo, b, &m, 0) );
00328 
00329   /* Do the test nt times... */
00330   for(iter = 0; iter < nt; iter++) {
00331 
00332     /* Choose a random value for x < a          */
00333     s_mp_pad(&x, USED(a));
00334     mpp_random(&x);
00335     MP_CHECKOK( mp_mod(&x, a, &x) );
00336 
00337     /* Compute z = (x ** m) mod a               */
00338     MP_CHECKOK( mp_exptmod(&x, &m, a, &z) );
00339     
00340     if(mp_cmp_d(&z, 1) == 0 || mp_cmp(&z, &amo) == 0) {
00341       res = MP_YES;
00342       continue;
00343     }
00344     
00345     res = MP_NO;  /* just in case the following for loop never executes. */
00346     for (jx = 1; jx < b; jx++) {
00347       /* z = z^2 (mod a) */
00348       MP_CHECKOK( mp_sqrmod(&z, a, &z) );
00349       res = MP_NO;   /* previous line set res to MP_YES */
00350 
00351       if(mp_cmp_d(&z, 1) == 0) {
00352        break;
00353       }
00354       if(mp_cmp(&z, &amo) == 0) {
00355        res = MP_YES;
00356        break;
00357       } 
00358     } /* end testing loop */
00359 
00360     /* If the test passes, we will continue iterating, but a failed
00361        test means the candidate is definitely NOT prime, so we will
00362        immediately break out of this loop
00363      */
00364     if(res == MP_NO)
00365       break;
00366 
00367   } /* end iterations loop */
00368   
00369 CLEANUP:
00370   mp_clear(&m);
00371   mp_clear(&z);
00372   mp_clear(&x);
00373   mp_clear(&amo);
00374   return res;
00375 
00376 } /* end mpp_pprime() */
00377 
00378 /* }}} */
00379 
00380 /* Produce table of composites from list of primes and trial value.  
00381 ** trial must be odd. List of primes must not include 2.
00382 ** sieve should have dimension >= MAXPRIME/2, where MAXPRIME is largest 
00383 ** prime in list of primes.  After this function is finished,
00384 ** if sieve[i] is non-zero, then (trial + 2*i) is composite.
00385 ** Each prime used in the sieve costs one division of trial, and eliminates
00386 ** one or more values from the search space. (3 eliminates 1/3 of the values
00387 ** alone!)  Each value left in the search space costs 1 or more modular 
00388 ** exponentations.  So, these divisions are a bargain!
00389 */
00390 mp_err mpp_sieve(mp_int *trial, const mp_digit *primes, mp_size nPrimes, 
00391                unsigned char *sieve, mp_size nSieve)
00392 {
00393   mp_err       res;
00394   mp_digit     rem;
00395   mp_size      ix;
00396   unsigned long offset;
00397 
00398   memset(sieve, 0, nSieve);
00399 
00400   for(ix = 0; ix < nPrimes; ix++) {
00401     mp_digit prime = primes[ix];
00402     mp_size  i;
00403     if((res = mp_mod_d(trial, prime, &rem)) != MP_OKAY) 
00404       return res;
00405 
00406     if (rem == 0) {
00407       offset = 0;
00408     } else {
00409       offset = prime - (rem / 2);
00410     }
00411     for (i = offset; i < nSieve ; i += prime) {
00412       sieve[i] = 1;
00413     }
00414   }
00415 
00416   return MP_OKAY;
00417 }
00418 
00419 #define SIEVE_SIZE 32*1024
00420 
00421 mp_err mpp_make_prime(mp_int *start, mp_size nBits, mp_size strong,
00422                     unsigned long * nTries)
00423 {
00424   mp_digit      np;
00425   mp_err        res;
00426   int           i    = 0;
00427   mp_int        trial;
00428   mp_int        q;
00429   mp_size       num_tests;
00430   unsigned char *sieve;
00431   
00432   sieve = malloc(SIEVE_SIZE);
00433   ARGCHK(sieve != NULL, MP_MEM);
00434 
00435   ARGCHK(start != 0, MP_BADARG);
00436   ARGCHK(nBits > 16, MP_RANGE);
00437 
00438   MP_DIGITS(&trial) = 0;
00439   MP_DIGITS(&q) = 0;
00440   MP_CHECKOK( mp_init(&trial) );
00441   MP_CHECKOK( mp_init(&q)     );
00442   /* values taken from table 4.4, HandBook of Applied Cryptography */
00443   if (nBits >= 1300) {
00444     num_tests = 2;
00445   } else if (nBits >= 850) {
00446     num_tests = 3;
00447   } else if (nBits >= 650) {
00448     num_tests = 4;
00449   } else if (nBits >= 550) {
00450     num_tests = 5;
00451   } else if (nBits >= 450) {
00452     num_tests = 6;
00453   } else if (nBits >= 400) {
00454     num_tests = 7;
00455   } else if (nBits >= 350) {
00456     num_tests = 8;
00457   } else if (nBits >= 300) {
00458     num_tests = 9;
00459   } else if (nBits >= 250) {
00460     num_tests = 12;
00461   } else if (nBits >= 200) {
00462     num_tests = 15;
00463   } else if (nBits >= 150) {
00464     num_tests = 18;
00465   } else if (nBits >= 100) {
00466     num_tests = 27;
00467   } else
00468     num_tests = 50;
00469 
00470   if (strong) 
00471     --nBits;
00472   MP_CHECKOK( mpl_set_bit(start, nBits - 1, 1) );
00473   MP_CHECKOK( mpl_set_bit(start,         0, 1) );
00474   for (i = mpl_significant_bits(start) - 1; i >= nBits; --i) {
00475     MP_CHECKOK( mpl_set_bit(start, i, 0) );
00476   }
00477   /* start sieveing with prime value of 3. */
00478   MP_CHECKOK(mpp_sieve(start, prime_tab + 1, prime_tab_size - 1, 
00479                      sieve, SIEVE_SIZE) );
00480 
00481 #ifdef DEBUG_SIEVE
00482   res = 0;
00483   for (i = 0; i < SIEVE_SIZE; ++i) {
00484     if (!sieve[i])
00485       ++res;
00486   }
00487   fprintf(stderr,"sieve found %d potential primes.\n", res);
00488 #define FPUTC(x,y) fputc(x,y)
00489 #else
00490 #define FPUTC(x,y) 
00491 #endif
00492 
00493   res = MP_NO;
00494   for(i = 0; i < SIEVE_SIZE; ++i) {
00495     if (sieve[i])    /* this number is composite */
00496       continue;
00497     MP_CHECKOK( mp_add_d(start, 2 * i, &trial) );
00498     FPUTC('.', stderr);
00499     /* run a Fermat test */
00500     res = mpp_fermat(&trial, 2);
00501     if (res != MP_OKAY) {
00502       if (res == MP_NO)
00503        continue;     /* was composite */
00504       goto CLEANUP;
00505     }
00506       
00507     FPUTC('+', stderr);
00508     /* If that passed, run some Miller-Rabin tests      */
00509     res = mpp_pprime(&trial, num_tests);
00510     if (res != MP_OKAY) {
00511       if (res == MP_NO)
00512        continue;     /* was composite */
00513       goto CLEANUP;
00514     }
00515     FPUTC('!', stderr);
00516 
00517     if (!strong) 
00518       break;  /* success !! */
00519 
00520     /* At this point, we have strong evidence that our candidate
00521        is itself prime.  If we want a strong prime, we need now
00522        to test q = 2p + 1 for primality...
00523      */
00524     MP_CHECKOK( mp_mul_2(&trial, &q) );
00525     MP_CHECKOK( mp_add_d(&q, 1, &q)  );
00526 
00527     /* Test q for small prime divisors ... */
00528     np = prime_tab_size;
00529     res = mpp_divis_primes(&q, &np);
00530     if (res == MP_YES) { /* is composite */
00531       mp_clear(&q);
00532       continue;
00533     }
00534     if (res != MP_NO) 
00535       goto CLEANUP;
00536 
00537     /* And test with Fermat, as with its parent ... */
00538     res = mpp_fermat(&q, 2);
00539     if (res != MP_YES) {
00540       mp_clear(&q);
00541       if (res == MP_NO)
00542        continue;     /* was composite */
00543       goto CLEANUP;
00544     }
00545 
00546     /* And test with Miller-Rabin, as with its parent ... */
00547     res = mpp_pprime(&q, num_tests);
00548     if (res != MP_YES) {
00549       mp_clear(&q);
00550       if (res == MP_NO)
00551        continue;     /* was composite */
00552       goto CLEANUP;
00553     }
00554 
00555     /* If it passed, we've got a winner */
00556     mp_exch(&q, &trial);
00557     mp_clear(&q);
00558     break;
00559 
00560   } /* end of loop through sieved values */
00561   if (res == MP_YES) 
00562     mp_exch(&trial, start);
00563 CLEANUP:
00564   mp_clear(&trial);
00565   mp_clear(&q);
00566   if (nTries)
00567     *nTries += i;
00568   if (sieve != NULL) {
00569        memset(sieve, 0, SIEVE_SIZE);
00570        free (sieve);
00571   }
00572   return res;
00573 }
00574 
00575 /*========================================================================*/
00576 /*------------------------------------------------------------------------*/
00577 /* Static functions visible only to the library internally                */
00578 
00579 /* {{{ s_mpp_divp(a, vec, size, which) */
00580 
00581 /* 
00582    Test for divisibility by members of a vector of digits.  Returns
00583    MP_NO if a is not divisible by any of them; returns MP_YES and sets
00584    'which' to the index of the offender, if it is.  Will stop on the
00585    first digit against which a is divisible.
00586  */
00587 
00588 mp_err    s_mpp_divp(mp_int *a, const mp_digit *vec, int size, int *which)
00589 {
00590   mp_err    res;
00591   mp_digit  rem;
00592 
00593   int     ix;
00594 
00595   for(ix = 0; ix < size; ix++) {
00596     if((res = mp_mod_d(a, vec[ix], &rem)) != MP_OKAY) 
00597       return res;
00598 
00599     if(rem == 0) {
00600       if(which)
00601        *which = ix;
00602       return MP_YES;
00603     }
00604   }
00605 
00606   return MP_NO;
00607 
00608 } /* end s_mpp_divp() */
00609 
00610 /* }}} */
00611 
00612 /*------------------------------------------------------------------------*/
00613 /* HERE THERE BE DRAGONS                                                  */