diff --git a/demo/test.c b/demo/test.c
index 939fcc0ed..d182d9a69 100644
--- a/demo/test.c
+++ b/demo/test.c
@@ -883,7 +883,7 @@ static int test_mp_prime_rand(void)
/* test for size */
for (ix = 10; ix < 128; ix++) {
- printf("Testing (not safe-prime): %9d bits \n", ix);
+ printf("\rTesting (not safe-prime): %9d bits ", ix);
fflush(stdout);
DO(mp_prime_rand(&a, 8, ix, (rand_int() & 1) ? 0 : MP_PRIME_2MSB_ON));
EXPECT(mp_count_bits(&a) == ix);
@@ -896,6 +896,238 @@ 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 */
+static const char *SPSP_2_100_LARGE[4] = {
+ "3L2x7YRmz7g4q+DwxESBacAClxrNiuspLCf8BUEphtky+5VNHLAb2ZZLLI0bu6cAOtNkUXenakBCCL"
+ "Vn7gqOpkcrQ/ptxZdk+4gnI99wFjgcfM512N71ZzbwvLe+5Pzat2k+nHIjE0w/WbQvzk4a2/syAY8S"
+ "i1B5XRjXYVAQOLyNWhsFpXeWXUgqiNzv7avfwBA3ZOXt", /* bases 2 - 100 */
+ "JOcSIwxGqGEjeQ2GsdlnFMwhc+xY7EtZo5Kf4BglOuakxTJaP8qrdZyduXaAZUdzyPgQLf7B8vqvVE"
+ "VLJwH7dLkLEiw19tfu3naT6DgQWzk+b5WuwWJzsTMdgWWH86M1h/Gjt2J/qABtTTH26C8bS4v/q9Fh"
+ "R8jqHNOiufUgHkDQdW9Z+BLlf6OVVh2VwPIOGVc7kFF", /* bases 2 - 107 */
+ "1ZCddPKHO7yeqI5ZeKG5ssTnzJeIDpWElJEZnHwejl4tsyly44JgwdiRmXgsi9FQfYhMzFZMgV6qWZZ"
+ "sIJl4RNgpD/PDb3nam++ECkzMBuNIXVpmZzw+Gj5xQmpKK+OX8pFSy2IQiKyKAOfSaivXEb2/dga2J/"
+ "Pc2d23lw+eP3WtBbfHc7TAQGgNI/6Xmcpl1G64eXCrJ", /* bases 2 - 103 */
+ "cCax282DurA+2Z54W3VLKSC2mwgpilQpGydCDHvXHNRKbJQRa5NtLLfa3sXvCmUWZ9okP2ZSsPDnw0X"
+ "dUQLzaz59vnw0rKbfsoA4nDBjMXR78Q889+KS4HFKfXkzxsiIKYo0kSfwPKYxFUi4Zj185kwwAPTAr2"
+ "IjegdWjQLeX1ZQM0HVUUF3WEVhHXcFzF0sMiJU5hl" /* bases 2 - 101 */
+};
+
+/* 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 digit 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"
+};
+
+const mp_digit prime_tab[] = {
+ 0x0002, 0x0003, 0x0005, 0x0007, 0x000B, 0x000D, 0x0011, 0x0013,
+ 0x0017, 0x001D, 0x001F, 0x0025, 0x0029, 0x002B, 0x002F, 0x0035,
+ 0x003B, 0x003D, 0x0043, 0x0047, 0x0049, 0x004F, 0x0053, 0x0059,
+ 0x0061, 0x0065, 0x0067, 0x006B, 0x006D, 0x0071, 0x007F, 0x0083,
+ 0x0089, 0x008B, 0x0095, 0x0097, 0x009D, 0x00A3, 0x00A7, 0x00AD,
+ 0x00B3, 0x00B5, 0x00BF, 0x00C1, 0x00C5, 0x00C7, 0x00D3, 0x00DF,
+ 0x00E3, 0x00E5, 0x00E9, 0x00EF, 0x00F1, 0x00FB, 0x0101, 0x0107,
+ 0x010D, 0x010F, 0x0115, 0x0119, 0x011B, 0x0125, 0x0133, 0x0137,
+
+ 0x0139, 0x013D, 0x014B, 0x0151, 0x015B, 0x015D, 0x0161, 0x0167,
+ 0x016F, 0x0175, 0x017B, 0x017F, 0x0185, 0x018D, 0x0191, 0x0199,
+ 0x01A3, 0x01A5, 0x01AF, 0x01B1, 0x01B7, 0x01BB, 0x01C1, 0x01C9,
+ 0x01CD, 0x01CF, 0x01D3, 0x01DF, 0x01E7, 0x01EB, 0x01F3, 0x01F7,
+ 0x01FD, 0x0209, 0x020B, 0x021D, 0x0223, 0x022D, 0x0233, 0x0239,
+ 0x023B, 0x0241, 0x024B, 0x0251, 0x0257, 0x0259, 0x025F, 0x0265,
+ 0x0269, 0x026B, 0x0277, 0x0281, 0x0283, 0x0287, 0x028D, 0x0293,
+ 0x0295, 0x02A1, 0x02A5, 0x02AB, 0x02B3, 0x02BD, 0x02C5, 0x02CF,
+
+ 0x02D7, 0x02DD, 0x02E3, 0x02E7, 0x02EF, 0x02F5, 0x02F9, 0x0301,
+ 0x0305, 0x0313, 0x031D, 0x0329, 0x032B, 0x0335, 0x0337, 0x033B,
+ 0x033D, 0x0347, 0x0355, 0x0359, 0x035B, 0x035F, 0x036D, 0x0371,
+ 0x0373, 0x0377, 0x038B, 0x038F, 0x0397, 0x03A1, 0x03A9, 0x03AD,
+ 0x03B3, 0x03B9, 0x03C7, 0x03CB, 0x03D1, 0x03D7, 0x03DF, 0x03E5,
+ 0x03F1, 0x03F5, 0x03FB, 0x03FD, 0x0407, 0x0409, 0x040F, 0x0419,
+ 0x041B, 0x0425, 0x0427, 0x042D, 0x043F, 0x0443, 0x0445, 0x0449,
+ 0x044F, 0x0455, 0x045D, 0x0463, 0x0469, 0x047F, 0x0481, 0x048B,
+
+ 0x0493, 0x049D, 0x04A3, 0x04A9, 0x04B1, 0x04BD, 0x04C1, 0x04C7,
+ 0x04CD, 0x04CF, 0x04D5, 0x04E1, 0x04EB, 0x04FD, 0x04FF, 0x0503,
+ 0x0509, 0x050B, 0x0511, 0x0515, 0x0517, 0x051B, 0x0527, 0x0529,
+ 0x052F, 0x0551, 0x0557, 0x055D, 0x0565, 0x0577, 0x0581, 0x058F,
+ 0x0593, 0x0595, 0x0599, 0x059F, 0x05A7, 0x05AB, 0x05AD, 0x05B3,
+ 0x05BF, 0x05C9, 0x05CB, 0x05CF, 0x05D1, 0x05D5, 0x05DB, 0x05E7,
+ 0x05F3, 0x05FB, 0x0607, 0x060D, 0x0611, 0x0617, 0x061F, 0x0623,
+ 0x062B, 0x062F, 0x063D, 0x0641, 0x0647, 0x0649, 0x064D, 0x0653
+};
+
+#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;
+ mp_digit j;
+ 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);
+ }
+
+ /* SPSP to bases 2 -- 100 */
+ mp_set(&b, 2u);
+ for (i = 0; i < 4; i++) {
+ DO(mp_read_radix(&a, SPSP_2_100_LARGE[i], 64));
+ for (j = 2u; j <= 100u; j++) {
+ result = false;
+ mp_set(&b, j);
+ DO(mp_prime_miller_rabin(&a, &b, &result));
+ EXPECT(result == true);
+ }
+ /* 107 is a prime that works */
+ mp_set(&b, 107u);
+ DO(mp_prime_miller_rabin(&a, &b, &result));
+ EXPECT(result == false);
+ }
+
+ /* SPSP to bases 2 -- 100, automatic */
+ mp_set(&b, 2u);
+ for (i = 0; i < 4; i++) {
+ DO(mp_read_radix(&a, SPSP_2_100_LARGE[i], 64));
+ for (j = 2u; j <= (mp_digit)mp_prime_rabin_miller_trials(mp_count_bits(&a)); j++) {
+ result = false;
+ mp_set(&b, (mp_digit)prime_tab[j]);
+ DO(mp_prime_miller_rabin(&a, &b, &result));
+ }
+ /* These numbers are not big enough for the heuristics to work */
+ 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;
@@ -905,6 +1137,17 @@ static int test_mp_prime_is_prime(void)
mp_int a, b;
DOR(mp_init_multi(&a, &b, NULL));
+ /* strong Miller-Rabin pseudoprimes to the first 100 primes (gernerated with Arnault's method) */
+ printf("Testing mp_prime_is_prime() with SPSPs to the first 100 primes\n");
+ for (ix = 0; ix < 4; ix++) {
+ DO(mp_read_radix(&a,SPSP_2_100_LARGE[ix],64));
+ DO(mp_prime_is_prime(&a, mp_prime_rabin_miller_trials(mp_count_bits(&a)), &cnt));
+ if (cnt) {
+ printf("SPSP_2_100_LARGE[%d] is not prime but mp_prime_is_prime says it is.\n", ix);
+ goto LBL_ERR;
+ }
+ }
+
/* strong Miller-Rabin pseudoprime to the first 200 primes (F. Arnault) */
printf("Testing mp_prime_is_prime() with Arnault's pseudoprime 803...901");
DO(mp_read_radix(&a,
@@ -965,7 +1208,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) {
@@ -2465,6 +2708,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 711ad0e58..1b64d3b08 100644
--- a/doc/bn.tex
+++ b/doc/bn.tex
@@ -2018,15 +2018,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}
@@ -2036,9 +2027,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
@@ -2220,9 +2208,8 @@ \subsection{Required Number of Tests}
the diagonal from $(512,2^{-80})$ downwards and to the right to gain a lower probability of getting
a composite declared a pseudoprime for the same amount of work or less.
-If this version of the library has the strong Lucas--Selfridge and/or the Frobenius--Underwood test
-implemented only one or two rounds of the Miller--Rabin test with a random base is necessary for
-numbers larger than or equal to $1024$ bits.
+If this version of the library has the extra strong Lucas test implemented only one or two rounds
+of the Miller--Rabin test with a random base is necessary for numbers larger than or equal to $1024$ bits.
This function is meant for RSA. The number of rounds for DSA is $\lceil -log_2(p)/2\rceil$ with $p$
the probability which is just the half of the absolute value of $p$ if given as a power of two.
@@ -2233,12 +2220,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.
@@ -2250,8 +2237,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 after 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
@@ -2272,11 +2258,11 @@ \section{Primality Testing}
mp_err mp_prime_is_prime(const mp_int *a, int t, bool *result)
\end{alltt}
This will perform a trial division followed by two rounds of Miller--Rabin with bases 2 and 3 and a
-Lucas--Selfridge test. The Frobenius--Underwood is available as a compile--time option with the
+extra strong Lucas test. The Frobenius--Underwood time is available as a compile--time option with the
preprocessor macro \texttt{LTM\_USE\_FROBENIUS\_TEST}. See file \texttt{bn\_mp\_prime\_is\_prime.c}
for the necessary details. It shall be noted that both functions are much slower than the
Miller--Rabin test and if speed is an essential issue, the macro \texttt{LTM\_USE\_ONLY\_MR}
-switches the Frobenius--Underwood test and the Lucas--Selfridge test off and their code will not
+switches the Frobenius--Underwood test and the Lucas test off and their code will not
even be compiled into the library.
If $t$ is set to a positive value $t$ additional rounds of the Miller--Rabin test with random bases
diff --git a/libtommath_VS2008.vcproj b/libtommath_VS2008.vcproj
index 13158a09d..bacc5249d 100644
--- a/libtommath_VS2008.vcproj
+++ b/libtommath_VS2008.vcproj
@@ -613,7 +613,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..79ff76ecc 100644
--- a/mp_prime_is_prime.c
+++ b/mp_prime_is_prime.c
@@ -16,9 +16,9 @@ static unsigned int s_floor_ilog2(int value)
mp_err mp_prime_is_prime(const mp_int *a, int t, bool *result)
{
mp_int b;
- int ix;
+ int ix, bits;
bool res;
- mp_err err;
+ mp_err err = MP_OKAY;
/* default to no */
*result = false;
@@ -59,11 +59,19 @@ mp_err mp_prime_is_prime(const mp_int *a, int t, bool *result)
if ((err = s_mp_prime_is_divisible(a, &res)) != MP_OKAY) {
return err;
}
-
/* return if it was trivially divisible */
if (res) {
return MP_OKAY;
}
+ /* floor(log_2(a)) */
+ bits = mp_count_bits(a) - 1;
+
+ /* If the whole prime table up to p = 1619 has been tested, than all
+ numbers below 1621^2 = 2,627,641 are prime now. log_2(1621^2) ~ 21.33 */
+ if (bits < 21) {
+ *result = true;
+ return MP_OKAY;
+ }
/*
Run the Miller-Rabin test with base 2 for the BPSW test.
@@ -78,6 +86,15 @@ mp_err mp_prime_is_prime(const mp_int *a, int t, bool *result)
if (!res) {
goto LBL_B;
}
+ /* If the whole prime table up to p = 1619 and the Miller-Rabin tests to base two
+ has been applied, than all numbers below 4,469,471 (~2^{22.1}) are prime now.
+ With 1659 SPSPs < 2^32 left */
+#if ((defined S_MP_PRIME_IS_DIVISIBLE_C) && (MP_PRIME_TAB_SIZE >= 256))
+ if (bits < 22) {
+ *result = true;
+ goto LBL_B;
+ }
+#endif
/*
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.
@@ -91,23 +108,36 @@ mp_err mp_prime_is_prime(const mp_int *a, int t, bool *result)
goto LBL_B;
}
+ /* If the whole prime table up to p = 1619 and the Miller-Rabin tests to bases
+ two and three have been applied, than all numbers below 11,541,307 (~2^{23.5}) are prime now.
+ With 89 SPSPs < 2^32 left */
+#if ((defined S_MP_PRIME_IS_DIVISIBLE_C) && (MP_PRIME_TAB_SIZE >= 256))
+ if (bits < 23) {
+ *result = true;
+ goto LBL_B;
+ }
+#endif
/*
- * Both, the Frobenius-Underwood test and the the Lucas-Selfridge test are quite
+ * Both, the Frobenius-Underwood test and the the extra strong Lucas 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)) {
+ if ((err = mp_prime_extra_strong_lucas(a, &res)) != MP_OKAY) {
goto LBL_B;
}
if (!res) {
goto LBL_B;
}
-#else
- if ((err = mp_prime_strong_lucas_selfridge(a, &res)) != MP_OKAY) {
+ /*
+ The Frobenius-Underwood pseudoprimes are sufficiently different from the
+ Extra Strong Lucas pseudoprimes with the parameters used in this library
+ to offer it as an additionally test (but it nearly doubles the runtime).
+ */
+#ifdef LTM_USE_FROBENIUS_TEST
+ err = mp_prime_frobenius_underwood(a, &res);
+ if ((err != MP_OKAY) && (err != MP_ITER)) {
goto LBL_B;
}
if (!res) {
@@ -117,8 +147,16 @@ mp_err mp_prime_is_prime(const mp_int *a, int t, bool *result)
}
#endif
- /* run at least one Miller-Rabin test with a random base */
+
if (t == 0) {
+#ifndef LTM_USE_ONLY_MR
+ /* The BPSW version as used here has no counter-example below 2^64 */
+ if (bits < 64) {
+ *result = true;
+ goto LBL_B;
+ }
+#endif
+ /* run at least one Miller-Rabin test with a random base if n > 2^64 */
t = 1;
}
@@ -127,35 +165,125 @@ mp_err mp_prime_is_prime(const mp_int *a, int t, bool *result)
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.
+ The caller has to check the maximum size.
*/
if (t < 0) {
int p_max = 0;
-
+#ifndef LTM_USE_ONLY_MR
+ if (bits < 64) {
+ /* Just complete the BPSW test */
+ if ((err = mp_prime_extra_strong_lucas(a, &res)) != MP_OKAY) {
+ goto LBL_B;
+ }
+ *result = res;
+ goto LBL_B;
+ }
+#else
/*
- Sorenson, Jonathan; Webster, Jonathan (2015).
- "Strong Pseudoprimes to Twelve Prime Bases".
+ Bases 2 and 3 are already done. Base 1459:
+ 1530787 = 2473 * 619 <- already out by trial division
+ 1518290707 = 19483 * 77929 <- trial division by 19483 (<2^15) or check for n == 1518290707 (~2^30.5)
+
+ This holds for a while. Next SPRPs < 2^35 to check for {2, 3, 1459}:
+
+ n factors log_2(n) has a factor < 2^28
+ 6770862367: 41143 164569 32.65669244751501848078 y
+ 15579919981: 88261 176521 33.85896877256133553111 y
+ 16149644101: 63541 254161 33.91078332064236217721 y
+ 17849326081: 50497 353473 34.05515055377487005769 y
+ 23510118061: 108421 216841 34.45256273267917105145 y
+ 24988416967: 79039 316153 34.54054045749318459634 y
+ 27031263841: 116257 232513 34.65390991522339129949 y
+ 28448982721: 97381 292141 34.72765801442453753728 y
+
*/
- /* 0x437ae92817f9fc85b7e5 = 318665857834031151167461 */
- if ((err = mp_read_radix(&b, "437ae92817f9fc85b7e5", 16)) != MP_OKAY) {
+#if ((defined S_MP_PRIME_IS_DIVISIBLE_C) && (MP_PRIME_TAB_SIZE >= 256))
+ if (bits < 32) {
+#if (MP_DIGIT_BIT < 31)
+ if ((err = mp_div_d(a, 19483u, NULL, &b)) != MP_OKAY) goto LBL_B;
+ if (mp_iszero(&b)) {
+ goto LBL_B;
+ }
+#else
+ if (mp_cmp_d(a, 1518290707u) == MP_EQ) {
+ goto LBL_B;
+ }
+#endif
+ mp_set(&b, 1459u);
+ if ((err = mp_prime_miller_rabin(a, &b, &res)) != MP_OKAY) {
+ goto LBL_B;
+ }
+ *result = res;
goto LBL_B;
}
+ }
+#else
+ /* 2, 7, 61 found by Gerhard Jaeschke 1993 */
+ mp_digit bases32 = {7u, 61u};
+#endif
+ /* 2, 325, 9375, 28178, 450775, 9780504, 1795265022 found by Jim Sinclair 2011 */
+ mp_word bases64 = {325ull, 9375ull, 28178ull, 450775ull, 9780504ull, 1795265022ull};
+ if (bits < 32) {
+ for (ix = 0; ix < 2; ix++) {
+ mp_set(&b, bases32[ix]);
+ if ((err = mp_prime_miller_rabin(a, &b, &res)) != MP_OKAY) {
+ goto LBL_B;
+ }
+ if (!res) {
+ goto LBL_B;
+ }
+ }
+ *result = true;
+ goto LBL_B;
+ } else if ((bits >= 32) && (bits < 64)) {
+ for (ix = 0; ix < 6; ix++) {
+ mp_set_u32(&b, bases64[ix]);
+ if ((err = mp_prime_miller_rabin(a, &b, &res)) != MP_OKAY) {
+ goto LBL_B;
+ }
+ if (!res) {
+ goto LBL_B;
+ }
+ }
+ *result = true;
+ goto LBL_B;
+ }
+#endif
+ /*
+ Sorenson, Jonathan; Webster, Jonathan, "Strong Pseudoprimes to Twelve Prime Bases". (2015) https://arxiv.org/abs/1509.00864
+ Z. Zhang, "Finding strong pseudoprimes to several bases," Math. Comp., 70:234 (2001) 863--87
+ Z. Zhang, "Finding C3-strong pseudoprimes," Math. Comp., 74:250 (2005) 1009--1024
+ Zhang, Zhenxiang, "Two kinds of strong pseudoprimes up to 10^36," Math. Comp., 76:260 (2007) 2095--2107
+ */
+ /*
+ Comparing bigints is not free, so give the magnitude of "n" a rough check
+ before spending computing time.
+ */
- if (mp_cmp(a, &b) == MP_LT) {
- p_max = 12;
- } else {
+ else if ((bits >= 64) && (bits <= 78)) {
+ /* 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 {
+ p_max = 13;
+ }
+ } else if ((bits >= 78) && (bits <= 81)) {
/* 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;
}
+ } else {
+ err = MP_VAL;
+ goto LBL_B;
}
/* we did bases 2 and 3 already, skip them */
@@ -173,10 +301,9 @@ mp_err mp_prime_is_prime(const mp_int *a, int t, bool *result)
Do "t" M-R tests with random bases between 3 and "a".
See Fips 186.4 p. 126ff
*/
- else if (t > 0) {
+ 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
@@ -185,16 +312,16 @@ mp_err mp_prime_is_prime(const mp_int *a, int t, bool *result)
* 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;
+ bits = mp_count_bits(a);
+ mask = (1u << s_floor_ilog2(bits)) - 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);
+ bits = (bits/10) * 7;
+ len = 2 * (bits * bits);
E.g.: a number of size 2^2048 would be reduced to the upper limit
@@ -224,7 +351,6 @@ mp_err mp_prime_is_prime(const mp_int *a, int t, bool *result)
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;
@@ -252,8 +378,8 @@ mp_err mp_prime_is_prime(const mp_int *a, int t, bool *result)
* smaller than "a"
*/
len = mp_count_bits(&b);
- if (len >= size_a) {
- len = (len - size_a) + 1;
+ if (len >= bits) {
+ len = (len - bits) + 1;
if ((err = mp_div_2d(&b, len, &b, NULL)) != MP_OKAY) {
goto LBL_B;
}
diff --git a/sources.cmake b/sources.cmake
index bbb2aeab6..bcb135241 100644
--- a/sources.cmake
+++ b/sources.cmake
@@ -77,7 +77,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 86f348727..42d84fd20 100644
--- a/tommath.def
+++ b/tommath.def
@@ -80,7 +80,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 7994b6e1a..0d33ea2eb 100644
--- a/tommath.h
+++ b/tommath.h
@@ -502,11 +502,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
*/
@@ -517,10 +512,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 e08bc5f3c..63410b18c 100644
--- a/tommath_class.h
+++ b/tommath_class.h
@@ -86,7 +86,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
@@ -605,12 +605,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)
@@ -642,8 +657,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