tommath/bn_mp_prime_is_prime.c

321 lines
9.7 KiB
C
Raw Blame History

This file contains ambiguous Unicode characters

This file contains Unicode characters that might be confused with other characters. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.

#include "tommath_private.h"
#ifdef BN_MP_PRIME_IS_PRIME_C
/* LibTomMath, multiple-precision integer library -- Tom St Denis
*
* LibTomMath is a library that provides multiple-precision
* integer arithmetic as well as number theoretic functionality.
*
* The library was designed directly after the MPI library by
* Michael Fromberger but has been written from scratch with
* additional optimizations in place.
*
* The library is free for all purposes without any express
* guarantee it works.
*/
// portable integer log of two with small footprint
static unsigned int floor_ilog2(int value)
{
unsigned int r = 0;
while ((value >>= 1) != 0) {
r++;
}
return r;
}
int mp_prime_is_prime(const mp_int *a, int t, int *result)
{
mp_int b;
int ix, err, res, p_max = 0, size_a, len;
unsigned int fips_rand, mask;
/* default to no */
*result = MP_NO;
/* valid value of t? */
if (t > PRIME_SIZE) {
puts("t > PRIME_SIZE");
return MP_VAL;
}
/* Some shortcuts */
/* N > 3 */
if (a->used == 1) {
if (a->dp[0] == 0 || a->dp[0] == 1) {
*result = 0;
return MP_OKAY;
}
if (a->dp[0] == 2) {
*result = 1;
return MP_OKAY;
}
}
/* N must be odd */
if (mp_iseven(a) == MP_YES) {
*result = 0;
return MP_OKAY;
}
/* N is not a perfect square: floor(sqrt(N))^2 != N */
if ((err = mp_is_square(a, &res)) != MP_OKAY) {
return err;
}
if (res != 0) {
*result = 0;
return MP_OKAY;
}
/* is the input equal to one of the primes in the table? */
for (ix = 0; ix < PRIME_SIZE; ix++) {
if (mp_cmp_d(a, ltm_prime_tab[ix]) == MP_EQ) {
*result = 1;
return MP_OKAY;
}
}
/* first perform trial division */
if ((err = mp_prime_is_divisible(a, &res)) != MP_OKAY) {
return err;
}
/* return if it was trivially divisible */
if (res == MP_YES) {
return MP_OKAY;
}
/*
Run the Miller-Rabin test with base 2 for the BPSW test.
*/
if ((err = mp_init_set(&b,2)) != MP_OKAY) {
return err;
}
if ((err = mp_prime_miller_rabin(a, &b, &res)) != MP_OKAY) {
goto LBL_B;
}
if (res == MP_NO) {
goto LBL_B;
}
/*
Rumours have it that Mathematica does a second M-R test with base 3.
Other rumours have it that their strong L-S test is slightly different.
It does not hurt, though, beside a bit of extra runtime.
*/
b.dp[0]++;
if ((err = mp_prime_miller_rabin(a, &b, &res)) != MP_OKAY) {
goto LBL_B;
}
if (res == MP_NO) {
goto LBL_B;
}
#ifdef MP_8BIT
if(t >= 0 && t < 8) {
t = 8;
}
#else
// commented out for testing purposes
//#ifdef LTM_USE_STRONG_LUCAS_SELFRIDGE_TEST
if ((err = mp_prime_strong_lucas_selfridge(a, &res)) != MP_OKAY) {
goto LBL_B;
}
if (res == MP_NO) {
goto LBL_B;
}
//#endif
// commented out for testing purposes
//#ifdef LTM_USE_FROBENIUS_UNDERWOOD_TEST
if ((err = mp_prime_frobenius_underwood(a, &res)) != MP_OKAY) {
goto LBL_B;
}
if (res == MP_NO) {
goto LBL_B;
}
//#endif
#endif
/*
abs(t) extra rounds of M-R to extend the range of primes it can find if t < 0.
Only recommended if the input range is known to be < 3317044064679887385961981
It uses the bases for a deterministic M-R test if input < 3317044064679887385961981
The caller has to check the size.
Not for cryptographic use because with known bases strong M-R pseudoprimes can
be constructed. Use at least one MM-R test with a random base (t >= 1).
The 1119 bit large number
80383745745363949125707961434194210813883768828755814583748891752229742737653\
33652186502336163960045457915042023603208766569966760987284043965408232928738\
79185086916685732826776177102938969773947016708230428687109997439976544144845\
34115587245063340927902227529622941498423068816854043264575340183297861112989\
60644845216191652872597534901
has been constructed by F. Arnault (F. Arnault, "Rabin-Miller primality test:
composite numbers which pass it.", Mathematics of Computation, 1995, 64. Jg.,
Nr. 209, S. 355-361), is a semiprime with the two factors
40095821663949960541830645208454685300518816604113250877450620473800321707011\
96242716223191597219733582163165085358166969145233813917169287527980445796800\
452592031836601
20047910831974980270915322604227342650259408302056625438725310236900160853505\
98121358111595798609866791081582542679083484572616906958584643763990222898400\
226296015918301
and it is a strong pseudoprime to all forty-six prime M-R bases up to 200
It does not fail the strong Bailley-PSP test as implemented here, it is just
given as an example, if not the reason to use the BPSW-test instead of M-R-tests
with a sequence of primes 2...n.
*/
if (t < 0) {
t = -t;
/*
Sorenson, Jonathan; Webster, Jonathan (2015).
"Strong Pseudoprimes to Twelve Prime Bases".
*/
/* 318665857834031151167461 */
if ((err = mp_read_radix(&b, "437ae92817f9fc85b7e5", 16)) != MP_OKAY) {
goto LBL_B;
}
if (mp_cmp(a,&b) == MP_LT) {
p_max = 12;
}
/* 3317044064679887385961981 */
if ((err = mp_read_radix(&b, "2be6951adc5b22410a5fd", 16)) != MP_OKAY) {
goto LBL_B;
}
if (mp_cmp(a,&b) == MP_LT) {
p_max = 13;
}
// for compatibility with the current API (well, compatible within a sign's width)
if (p_max < t) {
p_max = t;
}
if(p_max > PRIME_SIZE) {
err = MP_VAL;
goto LBL_B;
}
/* we did bases 2 and 3 already, skip them */
for (ix = 2; ix < p_max; ix++) {
mp_set(&b,ltm_prime_tab[ix]);
if ((err = mp_prime_miller_rabin(a, &b, &res)) != MP_OKAY) {
goto LBL_B;
}
if (res == MP_NO) {
goto LBL_B;
}
}
}
/*
Do "t" M-R tests with random bases between 3 and "a".
See Fips 186.4 p. 126ff
*/
else if (t > 0) {
// The mp_digit's have a defined bit-size but the size of the
// array a.dp is a simple 'int' and this library can not assume full
// compliance to the current C-standard (ISO/IEC 9899:2011) because
// it gets used for small embeded processors, too. Some of those MCUs
// have compilers that one cannot call standard compliant by any means.
// Hence the ugly type-fiddling in the following code.
size_a = mp_count_bits(a);
mask = (1u << floor_ilog2(size_a)) - 1u;
/*
Assuming the General Rieman hypothesis (never thought to write that in a
comment) the upper bound can be lowered to 2*(log a)^2.
E. Bach, “Explicit bounds for primality testing and related problems,”
Math. Comp. 55 (1990), 355380.
size_a = (size_a/10) * 7;
len = 2 * (size_a * size_a);
E.g.: a number of size 2^2048 would be reduced to the upper limit
floor(2048/10)*7 = 1428
2 * 1428^2 = 4078368
(would have been ~4030331.9962 with floats and natural log instead)
That number is smaller than 2^28, the default bit-size of mp_digit.
*/
/*
How many tests, you might ask? Dana Jacobsen of Math::Prime::Util fame
does exactly 1. In words: one. Look at the end of _GMP_is_prime() in
Math-Prime-Util-GMP-0.50/primality.c if you do not believe it.
The function mp_rand() goes to some length to use a cryptographically
good PRNG. That also means that the chance to always get the same base
in the loop is non-zero, although very low.
If the BPSW test and/or the addtional Frobenious test have been
performed instead of just the Miller-Rabin test with the bases 2 and 3,
a single extra test should suffice, so such a very unlikely event
will not do much harm.
To preemptivly answer the dangling question: no, a witness does not
need to be prime.
*/
for (ix = 0; ix < t; ix++) {
// mp_rand() guarantees the first digit to be non-zero
if ((err = mp_rand(&b, 1)) != MP_OKAY) {
goto LBL_B;
}
// Reduce digit before casting because mp_digit might be bigger than
// an unsigned int and "mask" on the other side is most probably not.
fips_rand = (unsigned int) (b.dp[0] & (mp_digit) mask);
#ifdef MP_8BIT
// One 8-bit digit is too small, so concatenate two if the size of
// unsigned int allows for it.
if( (sizeof(unsigned int) * CHAR_BIT)/2 >= (sizeof(mp_digit) * CHAR_BIT) ) {
if ((err = mp_rand(&b, 1)) != MP_OKAY) {
goto LBL_B;
}
fips_rand <<= sizeof(mp_digit) * CHAR_BIT;
fips_rand |= (unsigned int) b.dp[0];
}
#endif
// Ceil, because small numbers have a right to live, too,
len = (int) ( ((fips_rand & mask) + DIGIT_BIT) / DIGIT_BIT);
// Unlikely.
if(len < 0){
ix--;
continue;
}
if ((err = mp_rand(&b, len)) != MP_OKAY) {
goto LBL_B;
}
// Although the chance for b <= 3 is miniscule, try again.
if(mp_cmp_d(&b,3) != MP_GT) {
ix--;
continue;
}
if ((err = mp_prime_miller_rabin(a, &b, &res)) != MP_OKAY) {
goto LBL_B;
}
if (res == MP_NO) {
goto LBL_B;
}
}
}
/* passed the test */
*result = MP_YES;
LBL_B:
mp_clear(&b);
return err;
}
#endif
/* ref: $Format:%D$ */
/* git commit: $Format:%H$ */
/* commit time: $Format:%ai$ */