From 6e81627308c597bef5461ae39be830cd46846d73 Mon Sep 17 00:00:00 2001 From: czurnieden Date: Sun, 5 Mar 2023 08:16:59 +0100 Subject: [PATCH] Overhaul of the prime-tests - Removal of the Fermat test mp_prime_fermat - Replacement of the Strong Lucas-Selfridge test with the Extra Strong Lucas test with Robert Baillie's parameters P = 3 and Q = 1 - Finer grained early-outs - All determistic tests < 2^64 emprically verified - Additional tests to check the implementations of the Miller-Rabin and Extra Strong Lucas tests - Addition of tests of the LTM_USE_ONLY_MR to the CI - Documentation update --- .github/workflows/main.yml | 8 +- demo/test.c | 304 +++++++++++++++++++++++++++++++++- doc/bn.tex | 79 +++++---- libtommath_VS2008.vcproj | 2 +- makefile | 2 +- makefile.mingw | 2 +- makefile.msvc | 2 +- makefile.shared | 2 +- makefile.unix | 2 +- mp_prime_extra_strong_lucas.c | 137 +++++++++++++++ mp_prime_fermat.c | 41 ----- mp_prime_is_prime.c | 267 +++++++++++++++++++++++++---- sources.cmake | 2 +- tommath.def | 2 +- tommath.h | 10 +- tommath_class.h | 27 ++- 16 files changed, 748 insertions(+), 141 deletions(-) create mode 100644 mp_prime_extra_strong_lucas.c delete mode 100644 mp_prime_fermat.c diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml index 457216a7e..5e65b130a 100644 --- a/.github/workflows/main.yml +++ b/.github/workflows/main.yml @@ -44,13 +44,13 @@ jobs: # Run always with valgrind (no sanitizer, but debug info) - { BUILDOPTIONS: '--with-cc=gcc --with-m64 --with-valgrind', SANITIZER: '', COMPILE_DEBUG: '1', COMPILE_LTO: '0', CONV_WARNINGS: '', OTHERDEPS: '' } # Alternative big-int version of mp_log(_n) - - { BUILDOPTIONS: '--with-cc=gcc --with-m64 --cflags=-DS_MP_WORD_TOO_SMALL_C="" --with-valgrind', SANITIZER: '', COMPILE_DEBUG: '1', COMPILE_LTO: '0', CONV_WARNINGS: '', OTHERDEPS: '' } + - { BUILDOPTIONS: '--with-cc=gcc --with-m64 --cflags=-DS_MP_WORD_TOO_SMALL_C="" --cflags=-DLTM_USE_ONLY_MR --with-valgrind', SANITIZER: '', COMPILE_DEBUG: '1', COMPILE_LTO: '0', CONV_WARNINGS: '', OTHERDEPS: '' } # Shared library build - { BUILDOPTIONS: '--with-cc=gcc --make-option=-f --make-option=makefile.shared', SANITIZER: '', COMPILE_DEBUG: '0', COMPILE_LTO: '1', CONV_WARNINGS: '', OTHERDEPS: 'libtool-bin' } # GCC for the 32-bit architecture (no valgrind) - { BUILDOPTIONS: '--with-cc=gcc --with-m32', SANITIZER: '', COMPILE_DEBUG: '0', COMPILE_LTO: '0', CONV_WARNINGS: '', OTHERDEPS: 'libc6-dev-i386 gcc-multilib' } # Alternative big-int version of mp_log(_n) for the 32-bit architecture (no valgrind) - - { BUILDOPTIONS: '--with-cc=gcc --with-m32 --cflags=-DS_MP_WORD_TOO_SMALL_C="" ', SANITIZER: '', COMPILE_DEBUG: '0', COMPILE_LTO: '0', CONV_WARNINGS: '', OTHERDEPS: 'libc6-dev-i386 gcc-multilib' } + - { BUILDOPTIONS: '--with-cc=gcc --with-m32 --cflags=-DS_MP_WORD_TOO_SMALL_C="" --cflags=-DLTM_USE_ONLY_MR ', SANITIZER: '', COMPILE_DEBUG: '0', COMPILE_LTO: '0', CONV_WARNINGS: '', OTHERDEPS: 'libc6-dev-i386 gcc-multilib' } # clang for the 32-bit architecture (no valgrind) - { BUILDOPTIONS: '--with-cc=clang-10 --with-m32', SANITIZER: '', COMPILE_DEBUG: '0', COMPILE_LTO: '0', CONV_WARNINGS: '', OTHERDEPS: 'clang-10 llvm-10 gcc-multilib' } # RSA superclass with tests (no sanitizer, but debug info) @@ -108,8 +108,8 @@ jobs: - { BUILDOPTIONS: '--with-cc=gcc --cflags=-DMP_16BIT --limit-valgrind', SANITIZER: '1', COMPILE_DEBUG: '0', COMPILE_LTO: '0', CONV_WARNINGS: '', OTHERDEPS: '' } - { BUILDOPTIONS: '--with-cc=gcc --cflags=-DMP_32BIT --limit-valgrind', SANITIZER: '1', COMPILE_DEBUG: '0', COMPILE_LTO: '0', CONV_WARNINGS: '', OTHERDEPS: '' } # Alternative big-int version of mp_log(_n) - - { BUILDOPTIONS: '--with-cc=gcc --cflags=-DMP_16BIT --cflags=-DS_MP_WORD_TOO_SMALL_C="" --limit-valgrind', SANITIZER: '1', COMPILE_DEBUG: '0', COMPILE_LTO: '0', CONV_WARNINGS: '', OTHERDEPS: '' } - - { BUILDOPTIONS: '--with-cc=gcc --cflags=-DMP_32BIT --cflags=-DS_MP_WORD_TOO_SMALL_C="" --limit-valgrind', SANITIZER: '1', COMPILE_DEBUG: '0', COMPILE_LTO: '0', CONV_WARNINGS: '', OTHERDEPS: '' } + - { BUILDOPTIONS: '--with-cc=gcc --cflags=-DMP_16BIT --cflags=-DS_MP_WORD_TOO_SMALL_C="" --cflags=-DLTM_USE_ONLY_MR --limit-valgrind', SANITIZER: '1', COMPILE_DEBUG: '0', COMPILE_LTO: '0', CONV_WARNINGS: '', OTHERDEPS: '' } + - { BUILDOPTIONS: '--with-cc=gcc --cflags=-DMP_32BIT --cflags=-DS_MP_WORD_TOO_SMALL_C="" --cflags=-DLTM_USE_ONLY_MR --limit-valgrind', SANITIZER: '1', COMPILE_DEBUG: '0', COMPILE_LTO: '0', CONV_WARNINGS: '', OTHERDEPS: '' } # clang for the x86-64 architecture with restricted limb sizes - { BUILDOPTIONS: '--with-cc=clang --cflags=-DMP_16BIT --limit-valgrind', SANITIZER: '1', COMPILE_DEBUG: '0', COMPILE_LTO: '0', CONV_WARNINGS: '', OTHERDEPS: 'clang llvm' } diff --git a/demo/test.c b/demo/test.c index 939fcc0ed..27b992b02 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,15 +896,264 @@ 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 */ +}; + +#ifndef LTM_USE_ONLY_MR +/* 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 +}; +#endif + +/* 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; +} + +#ifndef LTM_USE_ONLY_MR +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; +} +#endif + static int test_mp_prime_is_prime(void) { int ix; mp_err e; - bool cnt, fu; + bool cnt; +#ifndef LTM_USE_ONLY_MR + bool fu; +#endif 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, @@ -948,6 +1197,7 @@ static int test_mp_prime_is_prime(void) DO(mp_prime_is_prime(&b, mp_prime_rabin_miller_trials(mp_count_bits(&b)), &cnt)); /* large problem */ EXPECT(cnt); +#ifndef LTM_USE_ONLY_MR DO(mp_prime_frobenius_underwood(&b, &fu)); EXPECT(fu); if ((e != MP_OKAY) || !cnt) { @@ -959,13 +1209,14 @@ static int test_mp_prime_is_prime(void) putchar('\n'); goto LBL_ERR; } - +#endif } +#ifndef LTM_USE_ONLY_MR /* Check regarding problem #143 */ 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) { @@ -974,6 +1225,47 @@ static int test_mp_prime_is_prime(void) putchar('\n'); goto LBL_ERR; } +#endif + /* Check deterministic tests */ +#ifdef LTM_USE_ONLY_MR +#if ((defined S_MP_PRIME_IS_DIVISIBLE_C) && (MP_PRIME_TAB_SIZE >= 256)) + /* 2-SPRP 4188889 = 431 * 9719 < 2^22 */ + DO(mp_read_radix(&a,"4188889",10)); + DO(mp_prime_is_prime(&a, 0, &cnt)); + EXPECT(cnt == false); + /* Last prime < 2^22 */ + DO(mp_read_radix(&a,"4194301",10)); + DO(mp_prime_is_prime(&a, 0, &cnt)); + EXPECT(cnt == true); + /* 2,3-SPRP 6787327 = 1303 * 5209 < 2^23 */ + DO(mp_read_radix(&a,"6787327",10)); + DO(mp_prime_is_prime(&a, 0, &cnt)); + EXPECT(cnt == false); + /* Last prime < 2^23 */ + DO(mp_read_radix(&a,"8388593",10)); + DO(mp_prime_is_prime(&a, 0, &cnt)); + EXPECT(cnt == true); + + /* 2,3,1459-SPRP < 2^32*/ + DO(mp_read_radix(&a,"1518290707",10)); + DO(mp_prime_is_prime(&a, -1, &cnt)); + EXPECT(cnt == false); +#endif + /* 2,3,7,61-SPRP < 2^43*/ + DO(mp_read_radix(&a,"7038007247701",10)); + DO(mp_prime_is_prime(&a, -1, &cnt)); + EXPECT(cnt == false); + + /* 2,325,9375,28178,450775,9780504-SPRP < 2^64 + which is also a + 2,3,325,9375,28178,450775,9780504-SPRP + */ + DO(mp_read_radix(&a,"18411296009130176041",10)); + DO(mp_prime_is_prime(&a, -1, &cnt)); + EXPECT(cnt == false); + +#endif + mp_clear_multi(&a, &b, NULL); return EXIT_SUCCESS; @@ -2465,6 +2757,10 @@ 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), +#ifndef LTM_USE_ONLY_MR + T1(mp_prime_extra_strong_lucas, MP_PRIME_EXTRA_STRONG_LUCAS), +#endif + 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..d5696112b 100644 --- a/doc/bn.tex +++ b/doc/bn.tex @@ -113,7 +113,7 @@ \section{Building LibTomMath} replace GCC for building the static and the shared library. Editing the makefiles is not needed, just set the shell variable \texttt{CC} as shown below. \begin{alltt} -CC=/home/czurnieden/intel/bin/icc make +CC=/home/$USER/intel/bin/icc make \end{alltt} ICC does not know all options available for GCC and LibTomMath uses two diagnostics @@ -1343,8 +1343,7 @@ \subsection{Multiplication by two} is implemented as a left--shift operation of $a$ by $b$ bits. It is also not very uncommon to need just the power of two $2^b$; for example as a start--value -for -the Newton method. +for the Newton method. \index{mp\_2expt} \begin{alltt} @@ -1440,7 +1439,7 @@ \section{Integer Division and Remainder} mp_err mp_div (const mp_int *a, const mp_int *b, mp_int *c, mp_int *d); \end{alltt} -This divides $a$ by $b$ and stores the quotient in $c$ and $d$. The signed quotient is computed +This divides $a$ by $b$ and stores the quotient in $c$ and the remainder in $d$. The signed quotient is computed such that $bc + d = a$. Note that either of $c$ or $d$ can be set to \texttt{NULL} if their value is not required. If $b$ is zero the function returns \texttt{MP\_VAL}. @@ -2004,7 +2003,7 @@ \subsection{Example} mp_error_to_string(e)); exit(EXIT_FAILURE); } - printf("%d\n",output); + printf("%d\textbackslash{}n",output); mp_clear(&x); exit(EXIT_SUCCESS); @@ -2018,15 +2017,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 +2026,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 +2207,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,13 +2219,13 @@ \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 -Rabin--Miller test with bases $2$ and $3$ resemble the BPSW test. The single internal use is a +Performs a extra strong Lucas test. The extra strong Lucas test together with the +Rabin--Miller test with base $2$ resembles 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 +2236,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 @@ -2271,13 +2256,8 @@ \section{Primality Testing} \begin{alltt} 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 -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 -even be compiled into the library. +This will perform a trial division followed by two rounds of Miller--Rabin to bases 2 and 3 and an +extra strong Lucas test. If $t$ is set to a positive value $t$ additional rounds of the Miller--Rabin test with random bases will be performed to allow for Fips 186.4 (vid.~p.~126ff) compliance. The function @@ -2289,16 +2269,33 @@ \section{Primality Testing} positive value this function will run $t + 1$ Miller--Rabin tests with random bases. If $t$ is set to a negative value the test will run the deterministic Miller--Rabin test for the -primes up to $3\,317\,044\,064\,679\,887\ 385\,961\,981$\footnote{The semiprime $1287836182261\cdot - 2575672364521$ with both factors smaller than $2^{64}$. An alternative with all factors smaller - than - $2^32$ is $4290067842\cdot 262853\cdot 1206721\cdot 2134439 + 3$}. That limit has to be checked -by -the caller. +primes up to $3\,317\,044\,064\,679\,887\ 385\,961\,981$ Input larger than the above limit +will return \texttt{MP\_VAL}. If $a$ passes all of the tests $result$ is set to \texttt{true}, otherwise it is set to \texttt{false}. +\subsection{Compile Time Switches} +There are several compile-time branches available. +\begin{description} +\item[\texttt{LTM\_USE\_ONLY\_MR}]\hfill \\ +Neither run the Frobenius--Underwood nor the Extra--Strong--Lucas test, only Miller-Rabin. +This does not change the deterministic tests with \texttt{t < 0} but the probabilistic test might +be a bit weaker. + +\item[\texttt{LTM\_USE\_FROBENIUS\_TEST}]\hfill \\ +This runs the Frobenius--Underwood test after the the Extra--Strong--Lucas test. This macro is not +compatible with the macro \texttt{LTM\_USE\_ONLY\_MR}. + +\item[\texttt{LTM\_USE\_ZHANG}]\hfill \\ +This macro allows the code to use larger bounds up to +$1\,543\,267\,864\,443\,420\,616\,877\,677\,640\,751\,301$ computed by Zhenxiang Zhang in +``Two Kinds of Strong Pseudoprimes up to $10^36$'' (Mathematics of computation, 76(260), 2095-2107) +but they have not been verified independently yet, use with caution. +\end{description} + + + \section{Next Prime} \index{mp\_prime\_next\_prime} \begin{alltt} @@ -2746,4 +2743,6 @@ \subsection{Shortcuts} \end{appendices} \input{bn.ind} + + \end{document} 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..7146d6513 100644 --- a/mp_prime_is_prime.c +++ b/mp_prime_is_prime.c @@ -16,9 +16,55 @@ 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; + +#ifdef LTM_USE_ONLY_MR + + /* 2, 7, 61 found by Gerhard Jaeschke 1993 (correctness confirmed by author) */ + /* + Because we already have 2 and 3 at this point all four bases are good up to + + 33717240721: 129841 259681 34.9728 y + + This is the single SPRP up to 51706044253 (~2^35.5896). Next are + + 163204128181: 285661 571321 37.2479 y + 501086407781: 288989 1733929 38.8663 y + ... 23588 others skipped ... + 18441334942415579101: 876577981 21037871521 63.9996 n + 18444384017352327673: 1920644893 9603224461 63.9998 n + + Adding base 5 brings us to + + 10087771603687: 1588063 6352249 43.1977 y + + The record for five bases is at 7999252175582851 for now, about 2^52.8288 bit large + but hase bases larger than 32 bit and none of the known 5-base records has bases + smaller than 28 bit. + + Adding base 63803 (already > 2^15) brings us to + + 849491953715047: 14573023 58292089 49.5936 y + + Record for 6 bases is at 585226005592931977 ~2^59.0217 but one base is over 60 bit + large. + + */ + mp_digit bases32[] = {5u, 7u, 61u}; + /* 2, 325, 9375, 28178, 450775, 9780504, 1795265022 found by Jim Sinclair 2011 + (correctness confirmed by author) + + An alternative would be 2, 3, 673, 325, 9375, 28178, 450775, 9780504 which has + 8 bases but all bases are smaller than 28 bit. + The first 2-SPRP > 2^64 is 18446744073709551617 which is not a 3-SPRP + */ + mp_word bases64[] = {325ull, 9375ull, 28178ull, 450775ull, 9780504ull, 1795265022ull}; +#if (MP_DIGIT_BIT < 31) + mp_digit rem; +#endif +#endif /* default to no */ *result = false; @@ -59,11 +105,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 +132,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 +154,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 SPRPs < 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 +193,19 @@ 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 is deterministic below 2^64 + (correctness confirmed by author) + */ + 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,36 +214,157 @@ 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 with {2, 3, 1459}: + + n factors log_2(n) has a factor < 2^28 + 6770862367: 41143 164569 32.6567 y + 15579919981: 88261 176521 33.859 y + 16149644101: 63541 254161 33.9108 y + 17849326081: 50497 353473 34.0556 y + 23510118061: 108421 216841 34.4526 y + 24988416967: 79039 316153 34.5405 y + 27031263841: 116257 232513 34.6539 y + 28448982721: 97381 292141 34.7276 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, &rem)) != MP_OKAY) goto LBL_B; + if (rem == 0u) { + 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; } +#endif + if (bits < 43) { + for (ix = 0; ix < 3; 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 >= 43) && (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 /* End of LTM_USE_ONLY_MR */ + /* + 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 - if (mp_cmp(a, &b) == MP_LT) { - p_max = 12; - } else { + Not implemented per default here but can be switched on via preprocessor directive "LTM_USE_ZHANG": + Zhang, Zhenxiang, "Two kinds of strong pseudoprimes up to 10^36," Math. Comp., 76:260 (2007) 2095--2107 + (10^36 > 2^119) + */ + 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 <= 82)) { /* 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; + } + /* The verified deterministic part ends here */ +#ifndef LTM_USE_ZHANG + else { + err = MP_VAL; + goto LBL_B; + } +#endif + } + /* The unverified deterministic part starts here */ +#ifdef LTM_USE_ZHANG + /* + USE WITH CAUTION! + + These bounds have not yet been independently verified. + But they are a good and relatively cheap method for sieving if an expensive full + deterministic primetest follows. (The Miller-Rabin test has no false negatives). + */ + /* 6003094289670105800312596501 = 54786377365501 * 109572754731001 ~2^92.2777 14 rounds */ + else if ((bits > 81) && (bits <= 92)) { + p_max = 14; + } + /* 59276361075595573263446330101 = 172157429516701 * 344314859033401 ~2^95.5814 15 rounds */ + else if ((bits > 92) && (bits <= 95)) { + p_max = 15; + } + /* 564132928021909221014087501701 = 531099297693901 * 1062198595387801 ~2^98.8319 16 rounds*/ + else if ((bits > 95) && (bits <= 98)) { + p_max = 16; + } + /* 1543267864443420616877677640751301 = 27778299663977101 * 55556599327954201 ~2^110.2496 18 rounds */ + else if ((bits > 98) && (bits <= 111)) { + if ((err = mp_read_radix(&b, "4c16c7697197146a6b8eb49518c5", 16)) != MP_OKAY) { + goto LBL_B; + } + if (mp_cmp(a, &b) == MP_LT) { + p_max = 18; } else { err = MP_VAL; goto LBL_B; } } +#endif + else { + err = MP_VAL; + goto LBL_B; + } /* we did bases 2 and 3 already, skip them */ for (ix = 2; ix < p_max; ix++) { @@ -173,10 +381,8 @@ 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 +391,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 +430,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 +457,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