/*
* 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() */
/* }}} */
static mpp_random_fn mpp_random_insecure = &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_insecure)(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)
{
return mpp_pprime_ext_random(a, nt, mpp_random_insecure);
}
mp_err
mpp_pprime_ext_random(mp_int *a,
int nt, mpp_random_fn random)
{
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)));
MP_CHECKOK((*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)
{
return mpp_make_prime_ext_random(start, nBits, strong, mpp_random_insecure);
}
mp_err
mpp_make_prime_ext_random(mp_int *start, mp_size nBits, mp_size strong, mpp_random_fn r
andom)
{
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_ext_random(&trial, num_tests, random);
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_ext_random(&q, num_tests, random);
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 (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 */
