diff options
Diffstat (limited to 'security/nss/lib/freebl/mpi/mpprime.c')
-rw-r--r-- | security/nss/lib/freebl/mpi/mpprime.c | 599 |
1 files changed, 599 insertions, 0 deletions
diff --git a/security/nss/lib/freebl/mpi/mpprime.c b/security/nss/lib/freebl/mpi/mpprime.c new file mode 100644 index 000000000..58287192e --- /dev/null +++ b/security/nss/lib/freebl/mpi/mpprime.c @@ -0,0 +1,599 @@ +/* + * mpprime.c + * + * Utilities for finding and working with prime and pseudo-prime + * integers + * + * This Source Code Form is subject to the terms of the Mozilla Public + * License, v. 2.0. If a copy of the MPL was not distributed with this + * file, You can obtain one at http://mozilla.org/MPL/2.0/. */ + +#include "mpi-priv.h" +#include "mpprime.h" +#include "mplogic.h" +#include <stdlib.h> +#include <string.h> + +#define SMALL_TABLE 0 /* determines size of hard-wired prime table */ + +#define RANDOM() rand() + +#include "primes.c" /* pull in the prime digit table */ + +/* + Test if any of a given vector of digits divides a. If not, MP_NO + is returned; otherwise, MP_YES is returned and 'which' is set to + the index of the integer in the vector which divided a. + */ +mp_err s_mpp_divp(mp_int *a, const mp_digit *vec, int size, int *which); + +/* {{{ mpp_divis(a, b) */ + +/* + mpp_divis(a, b) + + Returns MP_YES if a is divisible by b, or MP_NO if it is not. + */ + +mp_err +mpp_divis(mp_int *a, mp_int *b) +{ + mp_err res; + mp_int rem; + + if ((res = mp_init(&rem)) != MP_OKAY) + return res; + + if ((res = mp_mod(a, b, &rem)) != MP_OKAY) + goto CLEANUP; + + if (mp_cmp_z(&rem) == 0) + res = MP_YES; + else + res = MP_NO; + +CLEANUP: + mp_clear(&rem); + return res; + +} /* end mpp_divis() */ + +/* }}} */ + +/* {{{ mpp_divis_d(a, d) */ + +/* + mpp_divis_d(a, d) + + Return MP_YES if a is divisible by d, or MP_NO if it is not. + */ + +mp_err +mpp_divis_d(mp_int *a, mp_digit d) +{ + mp_err res; + mp_digit rem; + + ARGCHK(a != NULL, MP_BADARG); + + if (d == 0) + return MP_NO; + + if ((res = mp_mod_d(a, d, &rem)) != MP_OKAY) + return res; + + if (rem == 0) + return MP_YES; + else + return MP_NO; + +} /* end mpp_divis_d() */ + +/* }}} */ + +/* {{{ mpp_random(a) */ + +/* + mpp_random(a) + + Assigns a random value to a. This value is generated using the + standard C library's rand() function, so it should not be used for + cryptographic purposes, but it should be fine for primality testing, + since all we really care about there is good statistical properties. + + As many digits as a currently has are filled with random digits. + */ + +mp_err +mpp_random(mp_int *a) + +{ + mp_digit next = 0; + unsigned int ix, jx; + + ARGCHK(a != NULL, MP_BADARG); + + for (ix = 0; ix < USED(a); ix++) { + for (jx = 0; jx < sizeof(mp_digit); jx++) { + next = (next << CHAR_BIT) | (RANDOM() & UCHAR_MAX); + } + DIGIT(a, ix) = next; + } + + return MP_OKAY; + +} /* end mpp_random() */ + +/* }}} */ + +/* {{{ mpp_random_size(a, prec) */ + +mp_err +mpp_random_size(mp_int *a, mp_size prec) +{ + mp_err res; + + ARGCHK(a != NULL && prec > 0, MP_BADARG); + + if ((res = s_mp_pad(a, prec)) != MP_OKAY) + return res; + + return mpp_random(a); + +} /* end mpp_random_size() */ + +/* }}} */ + +/* {{{ mpp_divis_vector(a, vec, size, which) */ + +/* + mpp_divis_vector(a, vec, size, which) + + Determines if a is divisible by any of the 'size' digits in vec. + Returns MP_YES and sets 'which' to the index of the offending digit, + if it is; returns MP_NO if it is not. + */ + +mp_err +mpp_divis_vector(mp_int *a, const mp_digit *vec, int size, int *which) +{ + ARGCHK(a != NULL && vec != NULL && size > 0, MP_BADARG); + + return s_mpp_divp(a, vec, size, which); + +} /* end mpp_divis_vector() */ + +/* }}} */ + +/* {{{ mpp_divis_primes(a, np) */ + +/* + mpp_divis_primes(a, np) + + Test whether a is divisible by any of the first 'np' primes. If it + is, returns MP_YES and sets *np to the value of the digit that did + it. If not, returns MP_NO. + */ +mp_err +mpp_divis_primes(mp_int *a, mp_digit *np) +{ + int size, which; + mp_err res; + + ARGCHK(a != NULL && np != NULL, MP_BADARG); + + size = (int)*np; + if (size > prime_tab_size) + size = prime_tab_size; + + res = mpp_divis_vector(a, prime_tab, size, &which); + if (res == MP_YES) + *np = prime_tab[which]; + + return res; + +} /* end mpp_divis_primes() */ + +/* }}} */ + +/* {{{ mpp_fermat(a, w) */ + +/* + Using w as a witness, try pseudo-primality testing based on Fermat's + little theorem. If a is prime, and (w, a) = 1, then w^a == w (mod + a). So, we compute z = w^a (mod a) and compare z to w; if they are + equal, the test passes and we return MP_YES. Otherwise, we return + MP_NO. + */ +mp_err +mpp_fermat(mp_int *a, mp_digit w) +{ + mp_int base, test; + mp_err res; + + if ((res = mp_init(&base)) != MP_OKAY) + return res; + + mp_set(&base, w); + + if ((res = mp_init(&test)) != MP_OKAY) + goto TEST; + + /* Compute test = base^a (mod a) */ + if ((res = mp_exptmod(&base, a, a, &test)) != MP_OKAY) + goto CLEANUP; + + if (mp_cmp(&base, &test) == 0) + res = MP_YES; + else + res = MP_NO; + +CLEANUP: + mp_clear(&test); +TEST: + mp_clear(&base); + + return res; + +} /* end mpp_fermat() */ + +/* }}} */ + +/* + Perform the fermat test on each of the primes in a list until + a) one of them shows a is not prime, or + b) the list is exhausted. + Returns: MP_YES if it passes tests. + MP_NO if fermat test reveals it is composite + Some MP error code if some other error occurs. + */ +mp_err +mpp_fermat_list(mp_int *a, const mp_digit *primes, mp_size nPrimes) +{ + mp_err rv = MP_YES; + + while (nPrimes-- > 0 && rv == MP_YES) { + rv = mpp_fermat(a, *primes++); + } + return rv; +} + +/* {{{ mpp_pprime(a, nt) */ + +/* + mpp_pprime(a, nt) + + Performs nt iteration of the Miller-Rabin probabilistic primality + test on a. Returns MP_YES if the tests pass, MP_NO if one fails. + If MP_NO is returned, the number is definitely composite. If MP_YES + is returned, it is probably prime (but that is not guaranteed). + */ + +mp_err +mpp_pprime(mp_int *a, int nt) +{ + mp_err res; + mp_int x, amo, m, z; /* "amo" = "a minus one" */ + int iter; + unsigned int jx; + mp_size b; + + ARGCHK(a != NULL, MP_BADARG); + + MP_DIGITS(&x) = 0; + MP_DIGITS(&amo) = 0; + MP_DIGITS(&m) = 0; + MP_DIGITS(&z) = 0; + + /* Initialize temporaries... */ + MP_CHECKOK(mp_init(&amo)); + /* Compute amo = a - 1 for what follows... */ + MP_CHECKOK(mp_sub_d(a, 1, &amo)); + + b = mp_trailing_zeros(&amo); + if (!b) { /* a was even ? */ + res = MP_NO; + goto CLEANUP; + } + + MP_CHECKOK(mp_init_size(&x, MP_USED(a))); + MP_CHECKOK(mp_init(&z)); + MP_CHECKOK(mp_init(&m)); + MP_CHECKOK(mp_div_2d(&amo, b, &m, 0)); + + /* Do the test nt times... */ + for (iter = 0; iter < nt; iter++) { + + /* Choose a random value for 1 < x < a */ + MP_CHECKOK(s_mp_pad(&x, USED(a))); + mpp_random(&x); + MP_CHECKOK(mp_mod(&x, a, &x)); + if (mp_cmp_d(&x, 1) <= 0) { + iter--; /* don't count this iteration */ + continue; /* choose a new x */ + } + + /* Compute z = (x ** m) mod a */ + MP_CHECKOK(mp_exptmod(&x, &m, a, &z)); + + if (mp_cmp_d(&z, 1) == 0 || mp_cmp(&z, &amo) == 0) { + res = MP_YES; + continue; + } + + res = MP_NO; /* just in case the following for loop never executes. */ + for (jx = 1; jx < b; jx++) { + /* z = z^2 (mod a) */ + MP_CHECKOK(mp_sqrmod(&z, a, &z)); + res = MP_NO; /* previous line set res to MP_YES */ + + if (mp_cmp_d(&z, 1) == 0) { + break; + } + if (mp_cmp(&z, &amo) == 0) { + res = MP_YES; + break; + } + } /* end testing loop */ + + /* If the test passes, we will continue iterating, but a failed + test means the candidate is definitely NOT prime, so we will + immediately break out of this loop + */ + if (res == MP_NO) + break; + + } /* end iterations loop */ + +CLEANUP: + mp_clear(&m); + mp_clear(&z); + mp_clear(&x); + mp_clear(&amo); + return res; + +} /* end mpp_pprime() */ + +/* }}} */ + +/* Produce table of composites from list of primes and trial value. +** trial must be odd. List of primes must not include 2. +** sieve should have dimension >= MAXPRIME/2, where MAXPRIME is largest +** prime in list of primes. After this function is finished, +** if sieve[i] is non-zero, then (trial + 2*i) is composite. +** Each prime used in the sieve costs one division of trial, and eliminates +** one or more values from the search space. (3 eliminates 1/3 of the values +** alone!) Each value left in the search space costs 1 or more modular +** exponentations. So, these divisions are a bargain! +*/ +mp_err +mpp_sieve(mp_int *trial, const mp_digit *primes, mp_size nPrimes, + unsigned char *sieve, mp_size nSieve) +{ + mp_err res; + mp_digit rem; + mp_size ix; + unsigned long offset; + + memset(sieve, 0, nSieve); + + for (ix = 0; ix < nPrimes; ix++) { + mp_digit prime = primes[ix]; + mp_size i; + if ((res = mp_mod_d(trial, prime, &rem)) != MP_OKAY) + return res; + + if (rem == 0) { + offset = 0; + } else { + offset = prime - rem; + } + + for (i = offset; i < nSieve * 2; i += prime) { + if (i % 2 == 0) { + sieve[i / 2] = 1; + } + } + } + + return MP_OKAY; +} + +#define SIEVE_SIZE 32 * 1024 + +mp_err +mpp_make_prime(mp_int *start, mp_size nBits, mp_size strong, + unsigned long *nTries) +{ + mp_digit np; + mp_err res; + unsigned int i = 0; + mp_int trial; + mp_int q; + mp_size num_tests; + unsigned char *sieve; + + ARGCHK(start != 0, MP_BADARG); + ARGCHK(nBits > 16, MP_RANGE); + + sieve = malloc(SIEVE_SIZE); + ARGCHK(sieve != NULL, MP_MEM); + + MP_DIGITS(&trial) = 0; + MP_DIGITS(&q) = 0; + MP_CHECKOK(mp_init(&trial)); + MP_CHECKOK(mp_init(&q)); + /* values originally taken from table 4.4, + * HandBook of Applied Cryptography, augmented by FIPS-186 + * requirements, Table C.2 and C.3 */ + if (nBits >= 2000) { + num_tests = 3; + } else if (nBits >= 1536) { + num_tests = 4; + } else if (nBits >= 1024) { + num_tests = 5; + } else if (nBits >= 550) { + num_tests = 6; + } else if (nBits >= 450) { + num_tests = 7; + } else if (nBits >= 400) { + num_tests = 8; + } else if (nBits >= 350) { + num_tests = 9; + } else if (nBits >= 300) { + num_tests = 10; + } else if (nBits >= 250) { + num_tests = 20; + } else if (nBits >= 200) { + num_tests = 41; + } else if (nBits >= 100) { + num_tests = 38; /* funny anomaly in the FIPS tables, for aux primes, the + * required more iterations for larger aux primes */ + } else + num_tests = 50; + + if (strong) + --nBits; + MP_CHECKOK(mpl_set_bit(start, nBits - 1, 1)); + MP_CHECKOK(mpl_set_bit(start, 0, 1)); + for (i = mpl_significant_bits(start) - 1; i >= nBits; --i) { + MP_CHECKOK(mpl_set_bit(start, i, 0)); + } + /* start sieveing with prime value of 3. */ + MP_CHECKOK(mpp_sieve(start, prime_tab + 1, prime_tab_size - 1, + sieve, SIEVE_SIZE)); + +#ifdef DEBUG_SIEVE + res = 0; + for (i = 0; i < SIEVE_SIZE; ++i) { + if (!sieve[i]) + ++res; + } + fprintf(stderr, "sieve found %d potential primes.\n", res); +#define FPUTC(x, y) fputc(x, y) +#else +#define FPUTC(x, y) +#endif + + res = MP_NO; + for (i = 0; i < SIEVE_SIZE; ++i) { + if (sieve[i]) /* this number is composite */ + continue; + MP_CHECKOK(mp_add_d(start, 2 * i, &trial)); + FPUTC('.', stderr); + /* run a Fermat test */ + res = mpp_fermat(&trial, 2); + if (res != MP_OKAY) { + if (res == MP_NO) + continue; /* was composite */ + goto CLEANUP; + } + + FPUTC('+', stderr); + /* If that passed, run some Miller-Rabin tests */ + res = mpp_pprime(&trial, num_tests); + if (res != MP_OKAY) { + if (res == MP_NO) + continue; /* was composite */ + goto CLEANUP; + } + FPUTC('!', stderr); + + if (!strong) + break; /* success !! */ + + /* At this point, we have strong evidence that our candidate + is itself prime. If we want a strong prime, we need now + to test q = 2p + 1 for primality... + */ + MP_CHECKOK(mp_mul_2(&trial, &q)); + MP_CHECKOK(mp_add_d(&q, 1, &q)); + + /* Test q for small prime divisors ... */ + np = prime_tab_size; + res = mpp_divis_primes(&q, &np); + if (res == MP_YES) { /* is composite */ + mp_clear(&q); + continue; + } + if (res != MP_NO) + goto CLEANUP; + + /* And test with Fermat, as with its parent ... */ + res = mpp_fermat(&q, 2); + if (res != MP_YES) { + mp_clear(&q); + if (res == MP_NO) + continue; /* was composite */ + goto CLEANUP; + } + + /* And test with Miller-Rabin, as with its parent ... */ + res = mpp_pprime(&q, num_tests); + if (res != MP_YES) { + mp_clear(&q); + if (res == MP_NO) + continue; /* was composite */ + goto CLEANUP; + } + + /* If it passed, we've got a winner */ + mp_exch(&q, &trial); + mp_clear(&q); + break; + + } /* end of loop through sieved values */ + if (res == MP_YES) + mp_exch(&trial, start); +CLEANUP: + mp_clear(&trial); + mp_clear(&q); + if (nTries) + *nTries += i; + if (sieve != NULL) { + memset(sieve, 0, SIEVE_SIZE); + free(sieve); + } + return res; +} + +/*========================================================================*/ +/*------------------------------------------------------------------------*/ +/* Static functions visible only to the library internally */ + +/* {{{ s_mpp_divp(a, vec, size, which) */ + +/* + Test for divisibility by members of a vector of digits. Returns + MP_NO if a is not divisible by any of them; returns MP_YES and sets + 'which' to the index of the offender, if it is. Will stop on the + first digit against which a is divisible. + */ + +mp_err +s_mpp_divp(mp_int *a, const mp_digit *vec, int size, int *which) +{ + mp_err res; + mp_digit rem; + + int ix; + + for (ix = 0; ix < size; ix++) { + if ((res = mp_mod_d(a, vec[ix], &rem)) != MP_OKAY) + return res; + + if (rem == 0) { + if (which) + *which = ix; + return MP_YES; + } + } + + return MP_NO; + +} /* end s_mpp_divp() */ + +/* }}} */ + +/*------------------------------------------------------------------------*/ +/* HERE THERE BE DRAGONS */ |