diff --git a/demo/test.c b/demo/test.c
index 16fef5570..662fa7e05 100644
--- a/demo/test.c
+++ b/demo/test.c
@@ -824,6 +824,165 @@ static int test_mp_prime_rand(void)
return EXIT_FAILURE;
}
+/* Some small pseudoprimes to test the individual implementations */
+
+/* Miller-Rabin base 2 */
+static const uint32_t SPSP_2[] = {
+ 2047, 3277, 4033, 4681, 8321, 15841, 29341, 42799,
+ 49141, 52633, 65281, 74665, 80581, 85489, 88357, 90751
+};
+
+/* Miller-Rabin base 3 */
+static const uint32_t SPSP_3[] = {
+ 121, 703, 1891, 3281, 8401, 8911, 10585, 12403, 16531,
+ 18721, 19345, 23521, 31621, 44287, 47197, 55969, 63139,
+ 74593, 79003, 82513, 87913, 88573, 97567
+};
+
+/* SPSP to all bases < 100 */
+/* still needs computing
+static const char *SPSP_2_100_LARGE[] = {
+ "",
+ "",
+ "",
+ "",
+ ""
+};
+*/
+
+/* Extra strong Lucas test with Baillie's parameters Q = 1, P = 3 */
+static const uint32_t ESLPSP[] = {
+ 989, 3239, 5777, 10877, 27971, 29681, 30739, 31631, 39059, 72389,
+ 73919, 75077, 100127, 113573, 125249, 137549, 137801, 153931, 155819,
+ 161027, 162133, 189419, 218321, 231703, 249331, 370229, 429479, 430127,
+ 459191, 473891, 480689, 600059, 621781, 632249, 635627
+};
+
+/*
+ Almost extra strong Lucas test with Baillie's parameters Q = 1, P = 3
+ Only those that are not in ESLPSP.
+ */
+static const uint32_t AESLPSP[] = {
+ 10469, 154697, 233659, 472453, 629693, 852389, 1091093, 1560437,
+ 1620673, 1813601, 1969109, 2415739, 2595329, 2756837, 3721549,
+ 4269341, 5192309, 7045433, 7226669, 7265561
+};
+
+/* Some randomly choosen 200 decimal digits large primes (https://primes.utm.edu/lists/small/small2.html) */
+static const char *medium_primes[10] = {
+ "C8Ckh0vviS3HUPdB1NSrSm+gOodw/f1aQ5+aaH1W6RMB0jVkO6lTaL54O3o7U5BSGUFGxm5gAvisbJamasuLZS8g3ZsJ2JM4Vtn9cQZRfkP6b8V",
+ "64xDN9FqLBiovZ/9q/EPm0DONpIfn5MbJKHa+IjT0fjAzkg34FpAmad+CwhcpKaiTbZEpErut+DhpVyiQfqBFrgcGnGhhIrMF/XkyY3aVx6E96B",
+ "8cyuMlENm0vh/eWwgHUpDKqmLyCSsRQZRWvbHpA2jHDZv1EhHkVhceg3OFRZn/aXRBnbdtsc2xO6sWh9KZ5Mo7u9rJgBJMVtDnu094MCExj1YvB",
+ "BRFZFsYjSz45un8qptnuSqEsy9wV0BzbMpVAB1TrwImENOVIc1cASZNQ/mXG2xtazqgn/juVzFo91XLx9PtIlkcK0L2T6fBNgy8Lc7dSVoKQ+XP",
+ "Ez/mDl+to2gm69+VdIHI9Q7vaO3DuIdLVT69myM3HYwVBE+G24KffAOUAp3FGrSOU+LtERMiIYIEtxPI7n/DRJtmL2i0+REwGpTMge2d2EpabfB",
+ "5+Uz1gPFjZJ/nNdEOmOaMouJSGzygo42qz7xOwXn/moSUvBpPjo4twRGbK0+qaeU/RI8yYYxXr3OBP4w+/jgL3mN9GiENDM5LtEKMiQrZ9jIVEb",
+ "AQ5nD1+G1grv41s/XlK+0YTGyZgr/88PzdQJ8QT9tavisTgyG6k8/80A4HQhnFndskHNAaB2EW5fE7KH3kk7m89s8JnVqkJyGZWSfs1+JlmHLPf",
+ "3F19vPmM0Ih89KZ04Xmd62QB9F6E2sztT10A7Kcqc44eKvsNHh+JY6Z6gJXkbWg1Iw7xr29QAhEF/o1YAgfutQtpdzHkex06Yd71kPsaZdKXiC5",
+ "2fIcJ1t/VYCColXGs+ji/txNMEXn2FXdowLzlo7QKqzAWHdAbwtltSO5qpSp3OUiEOGUUi3hbyw3iQRE8nFJaikJ89Wdox6vpPtIsc3QRjexMnv",
+ "8aOicQ5gIbFCarFUgSgzh40LpuZ0jjK1u48/YT+C0h1dAQ8CIEgZjHZT+5/7cCRGmJlo+XCp7S41MSQ2ZNRSJh2texRYtvAXBAZfR8A8twl316P"
+};
+
+#define ARR_LENGTH(a) ((int)(sizeof((a))/sizeof((a)[0])))
+
+static int test_mp_prime_miller_rabin(void)
+{
+ mp_int a, b, c;
+ bool result;
+ int i;
+ DOR(mp_init_multi(&a, &b, &c, NULL));
+
+ /* SPSP to base 2 */
+ mp_set(&b, 2u);
+ for (i = 0; i < ARR_LENGTH(SPSP_2); i++) {
+ result = false;
+ mp_set_u32(&a, SPSP_2[i]);
+ DO(mp_prime_miller_rabin(&a, &b, &result));
+ EXPECT(result == true);
+ }
+
+ /* Some larger primes to check for false negatives */
+ for (i = 0; i < 10; i++) {
+ result = false;
+ DO(mp_read_radix(&a, medium_primes[i], 64));
+ DO(mp_prime_miller_rabin(&a, &b, &result));
+ EXPECT(result == true);
+ }
+ /* Some semi-primes */
+ for (i = 0; i < 5; i += 2) {
+ result = false;
+ DO(mp_read_radix(&a, medium_primes[i], 64));
+ DO(mp_read_radix(&c, medium_primes[i+1], 64));
+ DO(mp_mul(&a, &c, &a));
+ DO(mp_prime_miller_rabin(&a, &b, &result));
+ EXPECT(result == false);
+ }
+
+ /* SPSP to base 3 */
+ mp_set(&b, 3u);
+ for (i = 0; i < ARR_LENGTH(SPSP_3); i++) {
+ result = false;
+ mp_set_u32(&a, SPSP_3[i]);
+ DO(mp_prime_miller_rabin(&a, &b, &result));
+ EXPECT(result == true);
+ }
+
+ mp_clear_multi(&a, &b, &c, NULL);
+ return EXIT_SUCCESS;
+LBL_ERR:
+ mp_clear_multi(&a, &b, &c, NULL);
+ return EXIT_FAILURE;
+}
+
+
+static int test_mp_prime_extra_strong_lucas(void)
+{
+ mp_int a, b;
+ bool result;
+ int i;
+
+ DOR(mp_init_multi(&a, &b, NULL));
+
+ /* Check Extra Strong pseudoprimes */
+ for (i = 0; i < ARR_LENGTH(ESLPSP); i++) {
+ result = false;
+ mp_set_u32(&a, ESLPSP[i]);
+ DO(mp_prime_extra_strong_lucas(&a, &result));
+ EXPECT(result == true);
+ }
+
+ /* Check Almost Extra Strong pseudoprimes (not in ESLPSP) */
+ for (i = 0; i < ARR_LENGTH(AESLPSP); i++) {
+ result = false;
+ mp_set_u32(&a, AESLPSP[i]);
+ DO(mp_prime_extra_strong_lucas(&a, &result));
+ EXPECT(result == false);
+ }
+
+ /* Some larger primes to check for false negatives */
+ for (i = 0; i < 10; i++) {
+ result = false;
+ DO(mp_read_radix(&a, medium_primes[i], 64));
+ DO(mp_prime_extra_strong_lucas(&a, &result));
+ EXPECT(result == true);
+ }
+
+ /* Some semi-primes */
+ for (i = 0; i < 5; i++) {
+ result = false;
+ DO(mp_read_radix(&a, medium_primes[i], 64));
+ DO(mp_read_radix(&a, medium_primes[i+1], 64));
+ DO(mp_mul(&a, &b, &a));
+ DO(mp_prime_extra_strong_lucas(&a, &result));
+ EXPECT(result == false);
+ }
+
+ mp_clear_multi(&a, &b, NULL);
+ return EXIT_SUCCESS;
+LBL_ERR:
+ mp_clear_multi(&a, &b, NULL);
+ return EXIT_FAILURE;
+}
+
static int test_mp_prime_is_prime(void)
{
int ix;
@@ -893,7 +1052,7 @@ static int test_mp_prime_is_prime(void)
DO(mp_read_radix(&a,
"FFFFFFFFFFFFFFFFC90FDAA22168C234C4C6628B80DC1CD129024E088A67CC74020BBEA63B139B22514A08798E3404DDEF9519B3CD3A431B302B0A6DF25F14374FE1356D6D51C245E485B576625E7EC6F44C42E9A63A3620FFFFFFFFFFFFFFFF",
16));
- DO(mp_prime_strong_lucas_selfridge(&a, &cnt));
+ DO(mp_prime_extra_strong_lucas(&a, &cnt));
/* large problem */
EXPECT(cnt);
if ((e != MP_OKAY) || !cnt) {
@@ -1358,7 +1517,7 @@ static int test_mp_log_n(void)
mp_int a;
mp_digit d;
int base, lb, size;
- const int max_base = MP_MIN(INT_MAX, MP_DIGIT_MAX);
+ const mp_digit max_base = MP_MIN(INT_MAX, MP_DIGIT_MAX);
DOR(mp_init(&a));
@@ -2236,6 +2395,8 @@ static int unit_tests(int argc, char **argv)
T1(mp_montgomery_reduce, MP_MONTGOMERY_REDUCE),
T1(mp_root_n, MP_ROOT_N),
T1(mp_or, MP_OR),
+ T1(mp_prime_extra_strong_lucas, MP_PRIME_EXTRA_STRONG_LUCAS),
+ T1(mp_prime_miller_rabin, MP_PRIME_MILLER_RABIN),
T1(mp_prime_is_prime, MP_PRIME_IS_PRIME),
T1(mp_prime_next_prime, MP_PRIME_NEXT_PRIME),
T1(mp_prime_rand, MP_PRIME_RAND),
diff --git a/doc/bn.tex b/doc/bn.tex
index 566b3be32..5cc7bc62c 100644
--- a/doc/bn.tex
+++ b/doc/bn.tex
@@ -2014,15 +2014,6 @@ \subsection{Example}
\chapter{Prime Numbers}
-\section{Fermat Test}
-\index{mp\_prime\_fermat}
-\begin{alltt}
-mp_err mp_prime_fermat (const mp_int *a, const mp_int *b, int *result)
-\end{alltt}
-Performs a Fermat primality test to the base $b$. That is it computes $b^a \mbox{ mod }a$ and
-tests whether the value is equal to $b$ or not. If the values are equal then $a$ is probably prime
-and $result$ is set to one. Otherwise $result$ is set to zero.
-
\section{Miller--Rabin Test}
\index{mp\_prime\_miller\_rabin}
\begin{alltt}
@@ -2032,9 +2023,6 @@ \section{Miller--Rabin Test}
test and is very hard to fool (besides with Carmichael numbers). If $a$ passes the test (therefore
is probably prime) $result$ is set to one. Otherwise $result$ is set to zero.
-Note that it is suggested that you use the Miller--Rabin test instead of the Fermat test since all
-of the failures of Miller--Rabin are a subset of the failures of the Fermat test.
-
\subsection{Required Number of Tests}
Generally to ensure a number is very likely to be prime you have to perform the Miller--Rabin with
at least a half--dozen or so unique bases. However, it has been proven that the probability of
@@ -2229,12 +2217,12 @@ \subsection{Required Number of Tests}
See also table C.1 in FIPS 186-4.
-\section{Strong Lucas--Selfridge Test}
-\index{mp\_prime\_strong\_lucas\_selfridge}
+\section{Extra Strong Lucas Test}
+\index{mp\_prime\_extra\_strong\_lucas}
\begin{alltt}
-mp_err mp_prime_strong_lucas_selfridge(const mp_int *a, bool *result)
+mp_err mp_prime_extra_strong_lucas(const mp_int *a, bool *result)
\end{alltt}
-Performs a strong Lucas--Selfridge test. The strong Lucas--Selfridge test together with the
+Performs a extra strong Lucas test. The extra strong Lucas test together with the
Rabin--Miller test with bases $2$ and $3$ resemble the BPSW test. The single internal use is a
compile--time option in \texttt{mp\_prime\_is\_prime} and can be excluded from the Libtommath build
if not needed.
@@ -2246,8 +2234,7 @@ \section{Frobenius (Underwood) Test}
\end{alltt}
Performs the variant of the Frobenius test as described by Paul Underwood. It can be included at
build--time if the preprocessor macro \texttt{LTM\_USE\_FROBENIUS\_TEST} is defined and will be
-used
-instead of the Lucas--Selfridge test.
+used instead of the extra strong Lucas test.
It returns \texttt{MP\_ITER} if the number of iterations is exhausted, assumes a composite as the
input and sets \texttt{result} accordingly. This will reduce the set of available pseudoprimes by a
diff --git a/libtommath_VS2008.vcproj b/libtommath_VS2008.vcproj
index e27aa9860..a39776cd9 100644
--- a/libtommath_VS2008.vcproj
+++ b/libtommath_VS2008.vcproj
@@ -609,7 +609,7 @@
>
= 0; i--) {
+ if ((err = mp_mul(&Vk, &Vk1, &tmp1)) != MP_OKAY) goto LTM_ERR;
+ if ((err = mp_add(&tmp1, N, &tmp1)) != MP_OKAY) goto LTM_ERR;
+ if ((err = mp_sub_d(&tmp1, P, &tmp1)) != MP_OKAY) goto LTM_ERR;
+ if ((err = mp_mod(&tmp1, N, &tmp1)) != MP_OKAY) goto LTM_ERR;
+ if (s_mp_get_bit(&s, i) == 1) {
+ if ((err = mp_copy(&tmp1, &Vk)) != MP_OKAY) goto LTM_ERR;
+ if ((err = mp_sqr(&Vk1, &Vk1)) != MP_OKAY) goto LTM_ERR;
+ if ((err = mp_add(&Vk1, &Nm2, &Vk1)) != MP_OKAY) goto LTM_ERR;
+ if ((err = mp_mod(&Vk1, N, &Vk1)) != MP_OKAY) goto LTM_ERR;
+ } else {
+ if ((err = mp_copy(&tmp1, &Vk1)) != MP_OKAY) goto LTM_ERR;
+ if ((err = mp_sqr(&Vk, &Vk)) != MP_OKAY) goto LTM_ERR;
+ if ((err = mp_add(&Vk, &Nm2, &Vk)) != MP_OKAY) goto LTM_ERR;
+ if ((err = mp_mod(&Vk, N, &Vk)) != MP_OKAY) goto LTM_ERR;
+ }
+ }
+ /* Check Vk = \pm 2 mod N */
+ if ((mp_cmp_d(&Vk, 2u) == MP_EQ) || (mp_cmp(&Vk, &Nm2) == MP_EQ)) {
+ /* Compute Uk as described in [2] and check for Uk == 0 */
+ if ((err = mp_mul_d(&Vk, P, &tmp1)) != MP_OKAY) goto LTM_ERR;
+ if ((err = mp_mul_2(&Vk1, &Vk1)) != MP_OKAY) goto LTM_ERR;
+ if ((err = mp_sub(&tmp1, &Vk1, &tmp1)) != MP_OKAY) goto LTM_ERR;
+ tmp1.sign = MP_ZPOS;
+ if ((err = mp_mod(&tmp1, N, &Vk1)) != MP_OKAY) goto LTM_ERR;
+ if (mp_iszero(&Vk1)) {
+ *result = true;
+ goto LTM_ERR;
+ }
+ }
+ /* Check for V_{2^t k} = 0 mod N for some t 0 <= t < r - 1 ([3])*/
+ for (i = 0; i < (r - 1); i++) {
+ if (mp_iszero(&Vk)) {
+ *result = true;
+ goto LTM_ERR;
+ }
+ /* Loop: V(2k) = V(k)^2 - 2 and four minus two is two again */
+ if (mp_cmp_d(&Vk, 2u) == MP_EQ) {
+ goto LTM_ERR;
+ }
+ if ((err = mp_sqr(&Vk, &Vk)) != MP_OKAY) goto LTM_ERR;
+ if ((err = mp_sub_d(&Vk, 2u, &Vk)) != MP_OKAY) goto LTM_ERR;
+ if ((err = mp_mod(&Vk, N, &Vk)) != MP_OKAY) goto LTM_ERR;
+ }
+
+LTM_ERR:
+ mp_clear_multi(&Dz, &s, &Vk, &Vk1, &tmp1, &Nm2, NULL);
+ return err;
+}
+#endif
+
+
+
+#endif
diff --git a/mp_prime_fermat.c b/mp_prime_fermat.c
deleted file mode 100644
index ac8116fef..000000000
--- a/mp_prime_fermat.c
+++ /dev/null
@@ -1,41 +0,0 @@
-#include "tommath_private.h"
-#ifdef MP_PRIME_FERMAT_C
-/* LibTomMath, multiple-precision integer library -- Tom St Denis */
-/* SPDX-License-Identifier: Unlicense */
-
-/* performs one Fermat test.
- *
- * If "a" were prime then b**a == b (mod a) since the order of
- * the multiplicative sub-group would be phi(a) = a-1. That means
- * it would be the same as b**(a mod (a-1)) == b**1 == b (mod a).
- *
- * Sets result to 1 if the congruence holds, or zero otherwise.
- */
-mp_err mp_prime_fermat(const mp_int *a, const mp_int *b, bool *result)
-{
- mp_int t;
- mp_err err;
-
- /* ensure b > 1 */
- if (mp_cmp_d(b, 1uL) != MP_GT) {
- return MP_VAL;
- }
-
- /* init t */
- if ((err = mp_init(&t)) != MP_OKAY) {
- return err;
- }
-
- /* compute t = b**a mod a */
- if ((err = mp_exptmod(b, a, a, &t)) != MP_OKAY) {
- goto LBL_ERR;
- }
-
- /* is it equal to b? */
- *result = mp_cmp(&t, b) == MP_EQ;
-
-LBL_ERR:
- mp_clear(&t);
- return err;
-}
-#endif
diff --git a/mp_prime_is_prime.c b/mp_prime_is_prime.c
index bb24f5944..8e78a5497 100644
--- a/mp_prime_is_prime.c
+++ b/mp_prime_is_prime.c
@@ -107,7 +107,7 @@ mp_err mp_prime_is_prime(const mp_int *a, int t, bool *result)
goto LBL_B;
}
#else
- if ((err = mp_prime_strong_lucas_selfridge(a, &res)) != MP_OKAY) {
+ if ((err = mp_prime_extra_strong_lucas(a, &res)) != MP_OKAY) {
goto LBL_B;
}
if (!res) {
diff --git a/sources.cmake b/sources.cmake
index 797d461e2..465f56946 100644
--- a/sources.cmake
+++ b/sources.cmake
@@ -76,7 +76,7 @@ mp_neg.c
mp_or.c
mp_pack.c
mp_pack_count.c
-mp_prime_fermat.c
+mp_prime_extra_strong_lucas.c
mp_prime_frobenius_underwood.c
mp_prime_is_prime.c
mp_prime_miller_rabin.c
diff --git a/tommath.def b/tommath.def
index 1af4e9e79..1f79f0621 100644
--- a/tommath.def
+++ b/tommath.def
@@ -79,7 +79,7 @@ EXPORTS
mp_or
mp_pack
mp_pack_count
- mp_prime_fermat
+ mp_prime_extra_strong_lucas
mp_prime_frobenius_underwood
mp_prime_is_prime
mp_prime_miller_rabin
diff --git a/tommath.h b/tommath.h
index 44c92b22a..ae2e4470e 100644
--- a/tommath.h
+++ b/tommath.h
@@ -499,11 +499,6 @@ mp_err mp_hash(mp_int *a, mp_hval *hash) MP_WUR;
/* ---> Primes <--- */
-/* performs one Fermat test of "a" using base "b".
- * Sets result to 0 if composite or 1 if probable prime
- */
-mp_err mp_prime_fermat(const mp_int *a, const mp_int *b, bool *result) MP_WUR;
-
/* performs one Miller-Rabin test of "a" using base "b".
* Sets result to 0 if composite or 1 if probable prime
*/
@@ -514,10 +509,11 @@ mp_err mp_prime_miller_rabin(const mp_int *a, const mp_int *b, bool *result) MP_
*/
int mp_prime_rabin_miller_trials(int size) MP_WUR;
-/* performs one strong Lucas-Selfridge test of "a".
+
+/* performs one extra strong Lucas-Selfridge test of "a".
* Sets result to 0 if composite or 1 if probable prime
*/
-mp_err mp_prime_strong_lucas_selfridge(const mp_int *a, bool *result) MP_WUR;
+mp_err mp_prime_extra_strong_lucas(const mp_int *a, bool *result) MP_WUR;
/* performs one Frobenius test of "a" as described by Paul Underwood.
* Sets result to 0 if composite or 1 if probable prime
diff --git a/tommath_class.h b/tommath_class.h
index becbcb193..84599b07f 100644
--- a/tommath_class.h
+++ b/tommath_class.h
@@ -85,7 +85,7 @@
# define MP_OR_C
# define MP_PACK_C
# define MP_PACK_COUNT_C
-# define MP_PRIME_FERMAT_C
+# define MP_PRIME_EXTRA_STRONG_LUCAS_C
# define MP_PRIME_FROBENIUS_UNDERWOOD_C
# define MP_PRIME_IS_PRIME_C
# define MP_PRIME_MILLER_RABIN_C
@@ -582,12 +582,27 @@
# define MP_COUNT_BITS_C
#endif
-#if defined(MP_PRIME_FERMAT_C)
-# define MP_CLEAR_C
+#if defined(MP_PRIME_EXTRA_STRONG_LUCAS_C)
+# define MP_ADD_C
+# define MP_ADD_D_C
+# define MP_CLEAR_MULTI_C
# define MP_CMP_C
# define MP_CMP_D_C
-# define MP_EXPTMOD_C
-# define MP_INIT_C
+# define MP_CNT_LSB_C
+# define MP_COPY_C
+# define MP_COUNT_BITS_C
+# define MP_DIV_2D_C
+# define MP_INIT_MULTI_C
+# define MP_KRONECKER_C
+# define MP_MOD_C
+# define MP_MUL_2_C
+# define MP_MUL_C
+# define MP_MUL_D_C
+# define MP_SET_C
+# define MP_SET_U32_C
+# define MP_SUB_C
+# define MP_SUB_D_C
+# define S_MP_GET_BIT_C
#endif
#if defined(MP_PRIME_FROBENIUS_UNDERWOOD_C)
@@ -619,8 +634,8 @@
# define MP_DIV_2D_C
# define MP_INIT_SET_C
# define MP_IS_SQUARE_C
+# define MP_PRIME_EXTRA_STRONG_LUCAS_C
# define MP_PRIME_MILLER_RABIN_C
-# define MP_PRIME_STRONG_LUCAS_SELFRIDGE_C
# define MP_RAND_C
# define MP_READ_RADIX_C
# define MP_SET_C