1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/security/nss/lib/freebl/mpi/mpprime.c Wed Dec 31 06:09:35 2014 +0100
1.3 @@ -0,0 +1,584 @@
1.4 +/*
1.5 + * mpprime.c
1.6 + *
1.7 + * Utilities for finding and working with prime and pseudo-prime
1.8 + * integers
1.9 + *
1.10 + * This Source Code Form is subject to the terms of the Mozilla Public
1.11 + * License, v. 2.0. If a copy of the MPL was not distributed with this
1.12 + * file, You can obtain one at http://mozilla.org/MPL/2.0/. */
1.13 +
1.14 +#include "mpi-priv.h"
1.15 +#include "mpprime.h"
1.16 +#include "mplogic.h"
1.17 +#include <stdlib.h>
1.18 +#include <string.h>
1.19 +
1.20 +#define SMALL_TABLE 0 /* determines size of hard-wired prime table */
1.21 +
1.22 +#define RANDOM() rand()
1.23 +
1.24 +#include "primes.c" /* pull in the prime digit table */
1.25 +
1.26 +/*
1.27 + Test if any of a given vector of digits divides a. If not, MP_NO
1.28 + is returned; otherwise, MP_YES is returned and 'which' is set to
1.29 + the index of the integer in the vector which divided a.
1.30 + */
1.31 +mp_err s_mpp_divp(mp_int *a, const mp_digit *vec, int size, int *which);
1.32 +
1.33 +/* {{{ mpp_divis(a, b) */
1.34 +
1.35 +/*
1.36 + mpp_divis(a, b)
1.37 +
1.38 + Returns MP_YES if a is divisible by b, or MP_NO if it is not.
1.39 + */
1.40 +
1.41 +mp_err mpp_divis(mp_int *a, mp_int *b)
1.42 +{
1.43 + mp_err res;
1.44 + mp_int rem;
1.45 +
1.46 + if((res = mp_init(&rem)) != MP_OKAY)
1.47 + return res;
1.48 +
1.49 + if((res = mp_mod(a, b, &rem)) != MP_OKAY)
1.50 + goto CLEANUP;
1.51 +
1.52 + if(mp_cmp_z(&rem) == 0)
1.53 + res = MP_YES;
1.54 + else
1.55 + res = MP_NO;
1.56 +
1.57 +CLEANUP:
1.58 + mp_clear(&rem);
1.59 + return res;
1.60 +
1.61 +} /* end mpp_divis() */
1.62 +
1.63 +/* }}} */
1.64 +
1.65 +/* {{{ mpp_divis_d(a, d) */
1.66 +
1.67 +/*
1.68 + mpp_divis_d(a, d)
1.69 +
1.70 + Return MP_YES if a is divisible by d, or MP_NO if it is not.
1.71 + */
1.72 +
1.73 +mp_err mpp_divis_d(mp_int *a, mp_digit d)
1.74 +{
1.75 + mp_err res;
1.76 + mp_digit rem;
1.77 +
1.78 + ARGCHK(a != NULL, MP_BADARG);
1.79 +
1.80 + if(d == 0)
1.81 + return MP_NO;
1.82 +
1.83 + if((res = mp_mod_d(a, d, &rem)) != MP_OKAY)
1.84 + return res;
1.85 +
1.86 + if(rem == 0)
1.87 + return MP_YES;
1.88 + else
1.89 + return MP_NO;
1.90 +
1.91 +} /* end mpp_divis_d() */
1.92 +
1.93 +/* }}} */
1.94 +
1.95 +/* {{{ mpp_random(a) */
1.96 +
1.97 +/*
1.98 + mpp_random(a)
1.99 +
1.100 + Assigns a random value to a. This value is generated using the
1.101 + standard C library's rand() function, so it should not be used for
1.102 + cryptographic purposes, but it should be fine for primality testing,
1.103 + since all we really care about there is good statistical properties.
1.104 +
1.105 + As many digits as a currently has are filled with random digits.
1.106 + */
1.107 +
1.108 +mp_err mpp_random(mp_int *a)
1.109 +
1.110 +{
1.111 + mp_digit next = 0;
1.112 + unsigned int ix, jx;
1.113 +
1.114 + ARGCHK(a != NULL, MP_BADARG);
1.115 +
1.116 + for(ix = 0; ix < USED(a); ix++) {
1.117 + for(jx = 0; jx < sizeof(mp_digit); jx++) {
1.118 + next = (next << CHAR_BIT) | (RANDOM() & UCHAR_MAX);
1.119 + }
1.120 + DIGIT(a, ix) = next;
1.121 + }
1.122 +
1.123 + return MP_OKAY;
1.124 +
1.125 +} /* end mpp_random() */
1.126 +
1.127 +/* }}} */
1.128 +
1.129 +/* {{{ mpp_random_size(a, prec) */
1.130 +
1.131 +mp_err mpp_random_size(mp_int *a, mp_size prec)
1.132 +{
1.133 + mp_err res;
1.134 +
1.135 + ARGCHK(a != NULL && prec > 0, MP_BADARG);
1.136 +
1.137 + if((res = s_mp_pad(a, prec)) != MP_OKAY)
1.138 + return res;
1.139 +
1.140 + return mpp_random(a);
1.141 +
1.142 +} /* end mpp_random_size() */
1.143 +
1.144 +/* }}} */
1.145 +
1.146 +/* {{{ mpp_divis_vector(a, vec, size, which) */
1.147 +
1.148 +/*
1.149 + mpp_divis_vector(a, vec, size, which)
1.150 +
1.151 + Determines if a is divisible by any of the 'size' digits in vec.
1.152 + Returns MP_YES and sets 'which' to the index of the offending digit,
1.153 + if it is; returns MP_NO if it is not.
1.154 + */
1.155 +
1.156 +mp_err mpp_divis_vector(mp_int *a, const mp_digit *vec, int size, int *which)
1.157 +{
1.158 + ARGCHK(a != NULL && vec != NULL && size > 0, MP_BADARG);
1.159 +
1.160 + return s_mpp_divp(a, vec, size, which);
1.161 +
1.162 +} /* end mpp_divis_vector() */
1.163 +
1.164 +/* }}} */
1.165 +
1.166 +/* {{{ mpp_divis_primes(a, np) */
1.167 +
1.168 +/*
1.169 + mpp_divis_primes(a, np)
1.170 +
1.171 + Test whether a is divisible by any of the first 'np' primes. If it
1.172 + is, returns MP_YES and sets *np to the value of the digit that did
1.173 + it. If not, returns MP_NO.
1.174 + */
1.175 +mp_err mpp_divis_primes(mp_int *a, mp_digit *np)
1.176 +{
1.177 + int size, which;
1.178 + mp_err res;
1.179 +
1.180 + ARGCHK(a != NULL && np != NULL, MP_BADARG);
1.181 +
1.182 + size = (int)*np;
1.183 + if(size > prime_tab_size)
1.184 + size = prime_tab_size;
1.185 +
1.186 + res = mpp_divis_vector(a, prime_tab, size, &which);
1.187 + if(res == MP_YES)
1.188 + *np = prime_tab[which];
1.189 +
1.190 + return res;
1.191 +
1.192 +} /* end mpp_divis_primes() */
1.193 +
1.194 +/* }}} */
1.195 +
1.196 +/* {{{ mpp_fermat(a, w) */
1.197 +
1.198 +/*
1.199 + Using w as a witness, try pseudo-primality testing based on Fermat's
1.200 + little theorem. If a is prime, and (w, a) = 1, then w^a == w (mod
1.201 + a). So, we compute z = w^a (mod a) and compare z to w; if they are
1.202 + equal, the test passes and we return MP_YES. Otherwise, we return
1.203 + MP_NO.
1.204 + */
1.205 +mp_err mpp_fermat(mp_int *a, mp_digit w)
1.206 +{
1.207 + mp_int base, test;
1.208 + mp_err res;
1.209 +
1.210 + if((res = mp_init(&base)) != MP_OKAY)
1.211 + return res;
1.212 +
1.213 + mp_set(&base, w);
1.214 +
1.215 + if((res = mp_init(&test)) != MP_OKAY)
1.216 + goto TEST;
1.217 +
1.218 + /* Compute test = base^a (mod a) */
1.219 + if((res = mp_exptmod(&base, a, a, &test)) != MP_OKAY)
1.220 + goto CLEANUP;
1.221 +
1.222 +
1.223 + if(mp_cmp(&base, &test) == 0)
1.224 + res = MP_YES;
1.225 + else
1.226 + res = MP_NO;
1.227 +
1.228 + CLEANUP:
1.229 + mp_clear(&test);
1.230 + TEST:
1.231 + mp_clear(&base);
1.232 +
1.233 + return res;
1.234 +
1.235 +} /* end mpp_fermat() */
1.236 +
1.237 +/* }}} */
1.238 +
1.239 +/*
1.240 + Perform the fermat test on each of the primes in a list until
1.241 + a) one of them shows a is not prime, or
1.242 + b) the list is exhausted.
1.243 + Returns: MP_YES if it passes tests.
1.244 + MP_NO if fermat test reveals it is composite
1.245 + Some MP error code if some other error occurs.
1.246 + */
1.247 +mp_err mpp_fermat_list(mp_int *a, const mp_digit *primes, mp_size nPrimes)
1.248 +{
1.249 + mp_err rv = MP_YES;
1.250 +
1.251 + while (nPrimes-- > 0 && rv == MP_YES) {
1.252 + rv = mpp_fermat(a, *primes++);
1.253 + }
1.254 + return rv;
1.255 +}
1.256 +
1.257 +/* {{{ mpp_pprime(a, nt) */
1.258 +
1.259 +/*
1.260 + mpp_pprime(a, nt)
1.261 +
1.262 + Performs nt iteration of the Miller-Rabin probabilistic primality
1.263 + test on a. Returns MP_YES if the tests pass, MP_NO if one fails.
1.264 + If MP_NO is returned, the number is definitely composite. If MP_YES
1.265 + is returned, it is probably prime (but that is not guaranteed).
1.266 + */
1.267 +
1.268 +mp_err mpp_pprime(mp_int *a, int nt)
1.269 +{
1.270 + mp_err res;
1.271 + mp_int x, amo, m, z; /* "amo" = "a minus one" */
1.272 + int iter;
1.273 + unsigned int jx;
1.274 + mp_size b;
1.275 +
1.276 + ARGCHK(a != NULL, MP_BADARG);
1.277 +
1.278 + MP_DIGITS(&x) = 0;
1.279 + MP_DIGITS(&amo) = 0;
1.280 + MP_DIGITS(&m) = 0;
1.281 + MP_DIGITS(&z) = 0;
1.282 +
1.283 + /* Initialize temporaries... */
1.284 + MP_CHECKOK( mp_init(&amo));
1.285 + /* Compute amo = a - 1 for what follows... */
1.286 + MP_CHECKOK( mp_sub_d(a, 1, &amo) );
1.287 +
1.288 + b = mp_trailing_zeros(&amo);
1.289 + if (!b) { /* a was even ? */
1.290 + res = MP_NO;
1.291 + goto CLEANUP;
1.292 + }
1.293 +
1.294 + MP_CHECKOK( mp_init_size(&x, MP_USED(a)) );
1.295 + MP_CHECKOK( mp_init(&z) );
1.296 + MP_CHECKOK( mp_init(&m) );
1.297 + MP_CHECKOK( mp_div_2d(&amo, b, &m, 0) );
1.298 +
1.299 + /* Do the test nt times... */
1.300 + for(iter = 0; iter < nt; iter++) {
1.301 +
1.302 + /* Choose a random value for 1 < x < a */
1.303 + s_mp_pad(&x, USED(a));
1.304 + mpp_random(&x);
1.305 + MP_CHECKOK( mp_mod(&x, a, &x) );
1.306 + if(mp_cmp_d(&x, 1) <= 0) {
1.307 + iter--; /* don't count this iteration */
1.308 + continue; /* choose a new x */
1.309 + }
1.310 +
1.311 + /* Compute z = (x ** m) mod a */
1.312 + MP_CHECKOK( mp_exptmod(&x, &m, a, &z) );
1.313 +
1.314 + if(mp_cmp_d(&z, 1) == 0 || mp_cmp(&z, &amo) == 0) {
1.315 + res = MP_YES;
1.316 + continue;
1.317 + }
1.318 +
1.319 + res = MP_NO; /* just in case the following for loop never executes. */
1.320 + for (jx = 1; jx < b; jx++) {
1.321 + /* z = z^2 (mod a) */
1.322 + MP_CHECKOK( mp_sqrmod(&z, a, &z) );
1.323 + res = MP_NO; /* previous line set res to MP_YES */
1.324 +
1.325 + if(mp_cmp_d(&z, 1) == 0) {
1.326 + break;
1.327 + }
1.328 + if(mp_cmp(&z, &amo) == 0) {
1.329 + res = MP_YES;
1.330 + break;
1.331 + }
1.332 + } /* end testing loop */
1.333 +
1.334 + /* If the test passes, we will continue iterating, but a failed
1.335 + test means the candidate is definitely NOT prime, so we will
1.336 + immediately break out of this loop
1.337 + */
1.338 + if(res == MP_NO)
1.339 + break;
1.340 +
1.341 + } /* end iterations loop */
1.342 +
1.343 +CLEANUP:
1.344 + mp_clear(&m);
1.345 + mp_clear(&z);
1.346 + mp_clear(&x);
1.347 + mp_clear(&amo);
1.348 + return res;
1.349 +
1.350 +} /* end mpp_pprime() */
1.351 +
1.352 +/* }}} */
1.353 +
1.354 +/* Produce table of composites from list of primes and trial value.
1.355 +** trial must be odd. List of primes must not include 2.
1.356 +** sieve should have dimension >= MAXPRIME/2, where MAXPRIME is largest
1.357 +** prime in list of primes. After this function is finished,
1.358 +** if sieve[i] is non-zero, then (trial + 2*i) is composite.
1.359 +** Each prime used in the sieve costs one division of trial, and eliminates
1.360 +** one or more values from the search space. (3 eliminates 1/3 of the values
1.361 +** alone!) Each value left in the search space costs 1 or more modular
1.362 +** exponentations. So, these divisions are a bargain!
1.363 +*/
1.364 +mp_err mpp_sieve(mp_int *trial, const mp_digit *primes, mp_size nPrimes,
1.365 + unsigned char *sieve, mp_size nSieve)
1.366 +{
1.367 + mp_err res;
1.368 + mp_digit rem;
1.369 + mp_size ix;
1.370 + unsigned long offset;
1.371 +
1.372 + memset(sieve, 0, nSieve);
1.373 +
1.374 + for(ix = 0; ix < nPrimes; ix++) {
1.375 + mp_digit prime = primes[ix];
1.376 + mp_size i;
1.377 + if((res = mp_mod_d(trial, prime, &rem)) != MP_OKAY)
1.378 + return res;
1.379 +
1.380 + if (rem == 0) {
1.381 + offset = 0;
1.382 + } else {
1.383 + offset = prime - (rem / 2);
1.384 + }
1.385 + for (i = offset; i < nSieve ; i += prime) {
1.386 + sieve[i] = 1;
1.387 + }
1.388 + }
1.389 +
1.390 + return MP_OKAY;
1.391 +}
1.392 +
1.393 +#define SIEVE_SIZE 32*1024
1.394 +
1.395 +mp_err mpp_make_prime(mp_int *start, mp_size nBits, mp_size strong,
1.396 + unsigned long * nTries)
1.397 +{
1.398 + mp_digit np;
1.399 + mp_err res;
1.400 + int i = 0;
1.401 + mp_int trial;
1.402 + mp_int q;
1.403 + mp_size num_tests;
1.404 + unsigned char *sieve;
1.405 +
1.406 + ARGCHK(start != 0, MP_BADARG);
1.407 + ARGCHK(nBits > 16, MP_RANGE);
1.408 +
1.409 + sieve = malloc(SIEVE_SIZE);
1.410 + ARGCHK(sieve != NULL, MP_MEM);
1.411 +
1.412 + MP_DIGITS(&trial) = 0;
1.413 + MP_DIGITS(&q) = 0;
1.414 + MP_CHECKOK( mp_init(&trial) );
1.415 + MP_CHECKOK( mp_init(&q) );
1.416 + /* values taken from table 4.4, HandBook of Applied Cryptography */
1.417 + if (nBits >= 1300) {
1.418 + num_tests = 2;
1.419 + } else if (nBits >= 850) {
1.420 + num_tests = 3;
1.421 + } else if (nBits >= 650) {
1.422 + num_tests = 4;
1.423 + } else if (nBits >= 550) {
1.424 + num_tests = 5;
1.425 + } else if (nBits >= 450) {
1.426 + num_tests = 6;
1.427 + } else if (nBits >= 400) {
1.428 + num_tests = 7;
1.429 + } else if (nBits >= 350) {
1.430 + num_tests = 8;
1.431 + } else if (nBits >= 300) {
1.432 + num_tests = 9;
1.433 + } else if (nBits >= 250) {
1.434 + num_tests = 12;
1.435 + } else if (nBits >= 200) {
1.436 + num_tests = 15;
1.437 + } else if (nBits >= 150) {
1.438 + num_tests = 18;
1.439 + } else if (nBits >= 100) {
1.440 + num_tests = 27;
1.441 + } else
1.442 + num_tests = 50;
1.443 +
1.444 + if (strong)
1.445 + --nBits;
1.446 + MP_CHECKOK( mpl_set_bit(start, nBits - 1, 1) );
1.447 + MP_CHECKOK( mpl_set_bit(start, 0, 1) );
1.448 + for (i = mpl_significant_bits(start) - 1; i >= nBits; --i) {
1.449 + MP_CHECKOK( mpl_set_bit(start, i, 0) );
1.450 + }
1.451 + /* start sieveing with prime value of 3. */
1.452 + MP_CHECKOK(mpp_sieve(start, prime_tab + 1, prime_tab_size - 1,
1.453 + sieve, SIEVE_SIZE) );
1.454 +
1.455 +#ifdef DEBUG_SIEVE
1.456 + res = 0;
1.457 + for (i = 0; i < SIEVE_SIZE; ++i) {
1.458 + if (!sieve[i])
1.459 + ++res;
1.460 + }
1.461 + fprintf(stderr,"sieve found %d potential primes.\n", res);
1.462 +#define FPUTC(x,y) fputc(x,y)
1.463 +#else
1.464 +#define FPUTC(x,y)
1.465 +#endif
1.466 +
1.467 + res = MP_NO;
1.468 + for(i = 0; i < SIEVE_SIZE; ++i) {
1.469 + if (sieve[i]) /* this number is composite */
1.470 + continue;
1.471 + MP_CHECKOK( mp_add_d(start, 2 * i, &trial) );
1.472 + FPUTC('.', stderr);
1.473 + /* run a Fermat test */
1.474 + res = mpp_fermat(&trial, 2);
1.475 + if (res != MP_OKAY) {
1.476 + if (res == MP_NO)
1.477 + continue; /* was composite */
1.478 + goto CLEANUP;
1.479 + }
1.480 +
1.481 + FPUTC('+', stderr);
1.482 + /* If that passed, run some Miller-Rabin tests */
1.483 + res = mpp_pprime(&trial, num_tests);
1.484 + if (res != MP_OKAY) {
1.485 + if (res == MP_NO)
1.486 + continue; /* was composite */
1.487 + goto CLEANUP;
1.488 + }
1.489 + FPUTC('!', stderr);
1.490 +
1.491 + if (!strong)
1.492 + break; /* success !! */
1.493 +
1.494 + /* At this point, we have strong evidence that our candidate
1.495 + is itself prime. If we want a strong prime, we need now
1.496 + to test q = 2p + 1 for primality...
1.497 + */
1.498 + MP_CHECKOK( mp_mul_2(&trial, &q) );
1.499 + MP_CHECKOK( mp_add_d(&q, 1, &q) );
1.500 +
1.501 + /* Test q for small prime divisors ... */
1.502 + np = prime_tab_size;
1.503 + res = mpp_divis_primes(&q, &np);
1.504 + if (res == MP_YES) { /* is composite */
1.505 + mp_clear(&q);
1.506 + continue;
1.507 + }
1.508 + if (res != MP_NO)
1.509 + goto CLEANUP;
1.510 +
1.511 + /* And test with Fermat, as with its parent ... */
1.512 + res = mpp_fermat(&q, 2);
1.513 + if (res != MP_YES) {
1.514 + mp_clear(&q);
1.515 + if (res == MP_NO)
1.516 + continue; /* was composite */
1.517 + goto CLEANUP;
1.518 + }
1.519 +
1.520 + /* And test with Miller-Rabin, as with its parent ... */
1.521 + res = mpp_pprime(&q, num_tests);
1.522 + if (res != MP_YES) {
1.523 + mp_clear(&q);
1.524 + if (res == MP_NO)
1.525 + continue; /* was composite */
1.526 + goto CLEANUP;
1.527 + }
1.528 +
1.529 + /* If it passed, we've got a winner */
1.530 + mp_exch(&q, &trial);
1.531 + mp_clear(&q);
1.532 + break;
1.533 +
1.534 + } /* end of loop through sieved values */
1.535 + if (res == MP_YES)
1.536 + mp_exch(&trial, start);
1.537 +CLEANUP:
1.538 + mp_clear(&trial);
1.539 + mp_clear(&q);
1.540 + if (nTries)
1.541 + *nTries += i;
1.542 + if (sieve != NULL) {
1.543 + memset(sieve, 0, SIEVE_SIZE);
1.544 + free (sieve);
1.545 + }
1.546 + return res;
1.547 +}
1.548 +
1.549 +/*========================================================================*/
1.550 +/*------------------------------------------------------------------------*/
1.551 +/* Static functions visible only to the library internally */
1.552 +
1.553 +/* {{{ s_mpp_divp(a, vec, size, which) */
1.554 +
1.555 +/*
1.556 + Test for divisibility by members of a vector of digits. Returns
1.557 + MP_NO if a is not divisible by any of them; returns MP_YES and sets
1.558 + 'which' to the index of the offender, if it is. Will stop on the
1.559 + first digit against which a is divisible.
1.560 + */
1.561 +
1.562 +mp_err s_mpp_divp(mp_int *a, const mp_digit *vec, int size, int *which)
1.563 +{
1.564 + mp_err res;
1.565 + mp_digit rem;
1.566 +
1.567 + int ix;
1.568 +
1.569 + for(ix = 0; ix < size; ix++) {
1.570 + if((res = mp_mod_d(a, vec[ix], &rem)) != MP_OKAY)
1.571 + return res;
1.572 +
1.573 + if(rem == 0) {
1.574 + if(which)
1.575 + *which = ix;
1.576 + return MP_YES;
1.577 + }
1.578 + }
1.579 +
1.580 + return MP_NO;
1.581 +
1.582 +} /* end s_mpp_divp() */
1.583 +
1.584 +/* }}} */
1.585 +
1.586 +/*------------------------------------------------------------------------*/
1.587 +/* HERE THERE BE DRAGONS */
