283 lines
7.9 KiB
C
283 lines
7.9 KiB
C
#include "tommath_private.h"
|
|
#ifdef MP_PRIME_IS_PRIME_C
|
|
/* LibTomMath, multiple-precision integer library -- Tom St Denis */
|
|
/* SPDX-License-Identifier: Unlicense */
|
|
|
|
/* portable integer log of two with small footprint */
|
|
static unsigned int s_floor_ilog2(int value)
|
|
{
|
|
unsigned int r = 0;
|
|
while ((value >>= 1) != 0) {
|
|
r++;
|
|
}
|
|
return r;
|
|
}
|
|
|
|
mp_err mp_prime_is_prime(const mp_int *a, int t, bool *result)
|
|
{
|
|
mp_int b;
|
|
int ix;
|
|
bool res;
|
|
mp_err err;
|
|
|
|
/* default to no */
|
|
*result = false;
|
|
|
|
/* Some shortcuts */
|
|
/* N > 3 */
|
|
if (a->used == 1) {
|
|
if ((a->dp[0] == 0u) || (a->dp[0] == 1u)) {
|
|
*result = false;
|
|
return MP_OKAY;
|
|
}
|
|
if (a->dp[0] == 2u) {
|
|
*result = true;
|
|
return MP_OKAY;
|
|
}
|
|
}
|
|
|
|
/* N must be odd */
|
|
if (mp_iseven(a)) {
|
|
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) {
|
|
return MP_OKAY;
|
|
}
|
|
|
|
/* is the input equal to one of the primes in the table? */
|
|
for (ix = 0; ix < MP_PRIME_TAB_SIZE; ix++) {
|
|
if (mp_cmp_d(a, s_mp_prime_tab[ix]) == MP_EQ) {
|
|
*result = true;
|
|
return MP_OKAY;
|
|
}
|
|
}
|
|
/* first perform trial division */
|
|
if ((err = s_mp_prime_is_divisible(a, &res)) != MP_OKAY) {
|
|
return err;
|
|
}
|
|
|
|
/* return if it was trivially divisible */
|
|
if (res) {
|
|
return MP_OKAY;
|
|
}
|
|
|
|
/*
|
|
Run the Miller-Rabin test with base 2 for the BPSW test.
|
|
*/
|
|
if ((err = mp_init_set(&b, 2uL)) != MP_OKAY) {
|
|
return err;
|
|
}
|
|
|
|
if ((err = mp_prime_miller_rabin(a, &b, &res)) != MP_OKAY) {
|
|
goto LBL_B;
|
|
}
|
|
if (!res) {
|
|
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) {
|
|
goto LBL_B;
|
|
}
|
|
|
|
/*
|
|
* Both, the Frobenius-Underwood test and the the Lucas-Selfridge test are quite
|
|
* slow so if speed is an issue, define LTM_USE_ONLY_MR to use M-R tests with
|
|
* bases 2, 3 and t random bases.
|
|
*/
|
|
#ifndef LTM_USE_ONLY_MR
|
|
if (t >= 0) {
|
|
#ifdef LTM_USE_FROBENIUS_TEST
|
|
err = mp_prime_frobenius_underwood(a, &res);
|
|
if ((err != MP_OKAY) && (err != MP_ITER)) {
|
|
goto LBL_B;
|
|
}
|
|
if (!res) {
|
|
goto LBL_B;
|
|
}
|
|
#else
|
|
if ((err = mp_prime_strong_lucas_selfridge(a, &res)) != MP_OKAY) {
|
|
goto LBL_B;
|
|
}
|
|
if (!res) {
|
|
goto LBL_B;
|
|
}
|
|
#endif
|
|
}
|
|
#endif
|
|
|
|
/* run at least one Miller-Rabin test with a random base */
|
|
if (t == 0) {
|
|
t = 1;
|
|
}
|
|
|
|
/*
|
|
Only recommended if the input range is known to be < 3317044064679887385961981
|
|
|
|
It uses the bases necessary for a deterministic M-R test if the input is
|
|
smaller than 3317044064679887385961981
|
|
The caller has to check the size.
|
|
TODO: can be made a bit finer grained but comparing is not free.
|
|
*/
|
|
if (t < 0) {
|
|
int p_max = 0;
|
|
|
|
/*
|
|
Sorenson, Jonathan; Webster, Jonathan (2015).
|
|
"Strong Pseudoprimes to Twelve Prime Bases".
|
|
*/
|
|
/* 0x437ae92817f9fc85b7e5 = 318665857834031151167461 */
|
|
if ((err = mp_read_radix(&b, "437ae92817f9fc85b7e5", 16)) != MP_OKAY) {
|
|
goto LBL_B;
|
|
}
|
|
|
|
if (mp_cmp(a, &b) == MP_LT) {
|
|
p_max = 12;
|
|
} else {
|
|
/* 0x2be6951adc5b22410a5fd = 3317044064679887385961981 */
|
|
if ((err = mp_read_radix(&b, "2be6951adc5b22410a5fd", 16)) != MP_OKAY) {
|
|
goto LBL_B;
|
|
}
|
|
|
|
if (mp_cmp(a, &b) == MP_LT) {
|
|
p_max = 13;
|
|
} else {
|
|
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, s_mp_prime_tab[ix]);
|
|
if ((err = mp_prime_miller_rabin(a, &b, &res)) != MP_OKAY) {
|
|
goto LBL_B;
|
|
}
|
|
if (!res) {
|
|
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) {
|
|
unsigned int mask;
|
|
int size_a;
|
|
|
|
/*
|
|
* 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 << s_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), 355-380.
|
|
|
|
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++) {
|
|
unsigned int fips_rand;
|
|
int len;
|
|
|
|
/* 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);
|
|
if (fips_rand > (unsigned int)(INT_MAX - MP_DIGIT_BIT)) {
|
|
len = INT_MAX / MP_DIGIT_BIT;
|
|
} else {
|
|
len = (((int)fips_rand + MP_DIGIT_BIT) / MP_DIGIT_BIT);
|
|
}
|
|
/* Unlikely. */
|
|
if (len < 0) {
|
|
ix--;
|
|
continue;
|
|
}
|
|
if ((err = mp_rand(&b, len)) != MP_OKAY) {
|
|
goto LBL_B;
|
|
}
|
|
/*
|
|
* That number might got too big and the witness has to be
|
|
* smaller than "a"
|
|
*/
|
|
len = mp_count_bits(&b);
|
|
if (len >= size_a) {
|
|
len = (len - size_a) + 1;
|
|
if ((err = mp_div_2d(&b, len, &b, NULL)) != MP_OKAY) {
|
|
goto LBL_B;
|
|
}
|
|
}
|
|
/* Although the chance for b <= 3 is miniscule, try again. */
|
|
if (mp_cmp_d(&b, 3uL) != MP_GT) {
|
|
ix--;
|
|
continue;
|
|
}
|
|
if ((err = mp_prime_miller_rabin(a, &b, &res)) != MP_OKAY) {
|
|
goto LBL_B;
|
|
}
|
|
if (!res) {
|
|
goto LBL_B;
|
|
}
|
|
}
|
|
}
|
|
|
|
/* passed the test */
|
|
*result = true;
|
|
LBL_B:
|
|
mp_clear(&b);
|
|
return err;
|
|
}
|
|
|
|
#endif
|