From d56ffc08b8fbaa5049a19e2152a89c703efcbbc9 Mon Sep 17 00:00:00 2001 From: Chip Hogg Date: Thu, 14 Nov 2024 18:50:21 -0500 Subject: [PATCH 01/10] Add utilities to make Strong Lucas tests easier The Strong Lucas test coming in the next PR will already be complicated enough. It'll be convenient, and less distracting, if we already have functions for certain operations we'll need. One thing we'll need to do is detect inputs that are perfect squares. Fortunately, this is pretty easy to do robustly and quickly, with Newton's method. We _don't_ want to use `std::sqrt`, because that takes us into the floating point domain for no good reason, which could give us wrong answers for larger integers. The other thing we need is Jacobi symbols. These are a lot more obscure, but thankfully, still resonably straightforward to compute. The Wikipedia page (https://en.wikipedia.org/wiki/Jacobi_symbol) has a good explanation, and in particular, good instructions for computing values. With these utilities in place, the Strong Lucas code should be easier to review. --- src/core/include/mp-units/ext/prime.h | 77 +++++++++++++++++++++++++++ test/static/prime_test.cpp | 43 +++++++++++++++ 2 files changed, 120 insertions(+) diff --git a/src/core/include/mp-units/ext/prime.h b/src/core/include/mp-units/ext/prime.h index 51050f41..3f127eb4 100644 --- a/src/core/include/mp-units/ext/prime.h +++ b/src/core/include/mp-units/ext/prime.h @@ -183,6 +183,83 @@ struct NumberDecomposition { return false; } +// The Jacobi symbol, notated as `(a/n)`, is defined for odd positive `n` and any integer `a`, taking values +// in the set `{-1, 0, 1}`. Besides being a completely multiplicative function (so that, for example, both +// (a*b/n) = (a/n) * (b/n), and (a/n*m) = (a/n) * (a/m)), it obeys the following symmetry rules, which enable +// its calculation: +// +// 1. (a/1) = 1, and (1/n) = 1, for all a and n. +// +// 2. (a/n) = 0 whenever a and n have a nontrivial common factor. +// +// 3. (a/n) = (b/n) whenever (a % n) = (b % n). +// +// 4. (2a/n) = (a/n) if n % 8 = 1 or 7, and -(a/n) if n % 8 = 3 or 5. +// +// 5. (a/n) = (n/a) * x if a and n are both odd, positive, and coprime. Here, x is 1 if either (a % 4) = 1 +// or (n % 4) = 1, and -1 otherwise. +// +// 6. (-1/n) = 1 if n % 4 = 1, and -1 if n % 4 = 3. +[[nodiscard]] consteval int jacobi_symbol(int64_t raw_a, uint64_t n) +{ + // Rule 1: n=1 case. + if (n == 1u) { + return 1; + } + + // Starting conditions: transform `a` to strictly non-negative values, setting `result` to the sign that we + // pick up (if any) from following these rules (i.e., rules 3 and 6). + int result = ((raw_a >= 0) || (n % 4u == 1u)) ? 1 : -1; + auto a = static_cast(raw_a < 0 ? -raw_a : raw_a) % n; + + while (a != 0u) { + // Rule 4. + const int sign_for_even = (n % 8u == 1u || n % 8u == 7u) ? 1 : -1; + while (a % 2u == 0u) { + a /= 2u; + result *= sign_for_even; + } + + // Rule 1: a=1 case. + if (a == 1u) { + return result; + } + + // Rule 2. + if (std::gcd(a, n) != 1u) { + return 0; + } + + // Note that at this point, we know that `a` and `n` are coprime, and are both odd and positive. + // Therefore, we meet the preconditions for rule 5 (the "flip-and-reduce" rule). + result *= (n % 4u == 1u || a % 4u == 1u) ? 1 : -1; + const uint64_t new_a = n % a; + n = a; + a = new_a; + } + + return 0; +} + +[[nodiscard]] consteval bool is_perfect_square(uint64_t n) +{ + if (n < 2u) { + return true; + } + + uint64_t prev = n / 2u; + while (true) { + const uint64_t curr = (prev + n / prev) / 2u; + if (curr * curr == n) { + return true; + } + if (curr >= prev) { + return false; + } + prev = curr; + } +} + [[nodiscard]] consteval bool is_prime_by_trial_division(std::uintmax_t n) { for (std::uintmax_t f = 2; f * f <= n; f += 1 + (f % 2)) { diff --git a/test/static/prime_test.cpp b/test/static/prime_test.cpp index bfd7aae5..e669ab4e 100644 --- a/test/static/prime_test.cpp +++ b/test/static/prime_test.cpp @@ -122,4 +122,47 @@ static_assert(miller_rabin_probable_prime(2u, 9'007'199'254'740'881u), "Large kn static_assert(miller_rabin_probable_prime(2u, 18'446'744'073'709'551'557u), "Largest 64-bit prime"); +// Jacobi symbols --- a building block for the Strong Lucas probable prime test, needed for Baillie-PSW. +static_assert(jacobi_symbol(1, 1u) == 1, "Jacobi symbol always 1 when 'numerator' is 1"); +static_assert(jacobi_symbol(1, 3u) == 1, "Jacobi symbol always 1 when 'numerator' is 1"); +static_assert(jacobi_symbol(1, 5u) == 1, "Jacobi symbol always 1 when 'numerator' is 1"); +static_assert(jacobi_symbol(1, 987654321u) == 1, "Jacobi symbol always 1 when 'numerator' is 1"); + +static_assert(jacobi_symbol(3, 1u) == 1, "Jacobi symbol always 1 when 'denominator' is 1"); +static_assert(jacobi_symbol(5, 1u) == 1, "Jacobi symbol always 1 when 'denominator' is 1"); +static_assert(jacobi_symbol(-1234567890, 1u) == 1, "Jacobi symbol always 1 when 'denominator' is 1"); + +static_assert(jacobi_symbol(10, 5u) == 0, "Jacobi symbol always 0 when there's a common factor"); +static_assert(jacobi_symbol(25, 15u) == 0, "Jacobi symbol always 0 when there's a common factor"); +static_assert(jacobi_symbol(-24, 9u) == 0, "Jacobi symbol always 0 when there's a common factor"); + +static_assert(jacobi_symbol(14, 9u) == +jacobi_symbol(7, 9u), + "Divide numerator by 2: positive when (denom % 8) in {1, 7}"); +static_assert(jacobi_symbol(14, 15u) == +jacobi_symbol(7, 15u), + "Divide numerator by 2: positive when (denom % 8) in {1, 7}"); +static_assert(jacobi_symbol(14, 11u) == -jacobi_symbol(7, 11u), + "Divide numerator by 2: negative when (denom % 8) in {3, 5}"); +static_assert(jacobi_symbol(14, 13u) == -jacobi_symbol(7, 13u), + "Divide numerator by 2: negative when (denom % 8) in {3, 5}"); + +static_assert(jacobi_symbol(19, 9u) == +jacobi_symbol(9, 19u), "Flip is identity when (n % 4) = 1"); +static_assert(jacobi_symbol(17, 7u) == +jacobi_symbol(7, 17u), "Flip is identity when (a % 4) = 1"); +static_assert(jacobi_symbol(19, 7u) == -jacobi_symbol(9, 7u), "Flip changes sign when (n % 4) = 3 and (a % 4) = 3"); + +static_assert(jacobi_symbol(1001, 9907u) == -1, "Example from Wikipedia page"); +static_assert(jacobi_symbol(19, 45u) == 1, "Example from Wikipedia page"); +static_assert(jacobi_symbol(8, 21u) == -1, "Example from Wikipedia page"); +static_assert(jacobi_symbol(5, 21u) == 1, "Example from Wikipedia page"); + +// Tests for perfect square finder +static_assert(is_perfect_square(0u)); +static_assert(is_perfect_square(1u)); +static_assert(!is_perfect_square(2u)); +static_assert(is_perfect_square(4u)); + +constexpr uint64_t BIG_SQUARE = [](auto x) { return x * x; }((uint64_t{1u} << 32) - 1u); +static_assert(!is_perfect_square(BIG_SQUARE - 1u)); +static_assert(is_perfect_square(BIG_SQUARE)); +static_assert(!is_perfect_square(BIG_SQUARE + 1u)); + } // namespace From b49a4b711b1c8ecd2d06ed055668969e90e47a7c Mon Sep 17 00:00:00 2001 From: Chip Hogg Date: Fri, 15 Nov 2024 09:17:40 -0500 Subject: [PATCH 02/10] Add missing `std::` prefixes --- src/core/include/mp-units/ext/prime.h | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/core/include/mp-units/ext/prime.h b/src/core/include/mp-units/ext/prime.h index 3f127eb4..7d6f7c5b 100644 --- a/src/core/include/mp-units/ext/prime.h +++ b/src/core/include/mp-units/ext/prime.h @@ -200,7 +200,7 @@ struct NumberDecomposition { // or (n % 4) = 1, and -1 otherwise. // // 6. (-1/n) = 1 if n % 4 = 1, and -1 if n % 4 = 3. -[[nodiscard]] consteval int jacobi_symbol(int64_t raw_a, uint64_t n) +[[nodiscard]] consteval int jacobi_symbol(std::int64_t raw_a, std::uint64_t n) { // Rule 1: n=1 case. if (n == 1u) { @@ -210,7 +210,7 @@ struct NumberDecomposition { // Starting conditions: transform `a` to strictly non-negative values, setting `result` to the sign that we // pick up (if any) from following these rules (i.e., rules 3 and 6). int result = ((raw_a >= 0) || (n % 4u == 1u)) ? 1 : -1; - auto a = static_cast(raw_a < 0 ? -raw_a : raw_a) % n; + auto a = static_cast(raw_a < 0 ? -raw_a : raw_a) % n; while (a != 0u) { // Rule 4. @@ -233,7 +233,7 @@ struct NumberDecomposition { // Note that at this point, we know that `a` and `n` are coprime, and are both odd and positive. // Therefore, we meet the preconditions for rule 5 (the "flip-and-reduce" rule). result *= (n % 4u == 1u || a % 4u == 1u) ? 1 : -1; - const uint64_t new_a = n % a; + const std::uint64_t new_a = n % a; n = a; a = new_a; } @@ -241,15 +241,15 @@ struct NumberDecomposition { return 0; } -[[nodiscard]] consteval bool is_perfect_square(uint64_t n) +[[nodiscard]] consteval bool is_perfect_square(std::uint64_t n) { if (n < 2u) { return true; } - uint64_t prev = n / 2u; + std::uint64_t prev = n / 2u; while (true) { - const uint64_t curr = (prev + n / prev) / 2u; + const std::uint64_t curr = (prev + n / prev) / 2u; if (curr * curr == n) { return true; } From 26639dab343a66c7d3dd168900681c7f6e0c8d41 Mon Sep 17 00:00:00 2001 From: Chip Hogg Date: Fri, 15 Nov 2024 09:39:39 -0500 Subject: [PATCH 03/10] Add some more --- test/static/prime_test.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/static/prime_test.cpp b/test/static/prime_test.cpp index e669ab4e..5b279527 100644 --- a/test/static/prime_test.cpp +++ b/test/static/prime_test.cpp @@ -160,7 +160,7 @@ static_assert(is_perfect_square(1u)); static_assert(!is_perfect_square(2u)); static_assert(is_perfect_square(4u)); -constexpr uint64_t BIG_SQUARE = [](auto x) { return x * x; }((uint64_t{1u} << 32) - 1u); +constexpr std::uint64_t BIG_SQUARE = [](auto x) { return x * x; }((std::uint64_t{1u} << 32) - 1u); static_assert(!is_perfect_square(BIG_SQUARE - 1u)); static_assert(is_perfect_square(BIG_SQUARE)); static_assert(!is_perfect_square(BIG_SQUARE + 1u)); From 9e8dfec265344e8eb4975fd6a8836f9ccf4b1ba0 Mon Sep 17 00:00:00 2001 From: Chip Hogg Date: Thu, 14 Nov 2024 20:54:31 -0500 Subject: [PATCH 04/10] Add the Strong Lucas Probable Prime test This is more involved than the Miller-Rabin test, but we can tame the complexity by breaking it down into helper functions, performing tasks such as: - Increment the index of the (U_k, V_k) sequence elements by one. - Double the index of the (U_k, V_k) sequence elements. - Find an appropriate D parameter. etc. With these helpers, the algorithm becomes straightforward (see, for instance, https://en.wikipedia.org/wiki/Lucas_pseudoprime#Strong_Lucas_pseudoprimes). We start by ruling out perfect squares (because if we don't, then the search for `D` will never terminate). Then we find our `D`, and decompose `n + 1` into `s` and `d` parameters (exactly as we did for Miller-Rabin, except there we used `n - 1`). At this point, the strong test is easy: check whether `U_d` is 0, then check `V_d`, as well as `V` for all successive doublings of the index less than `s`. A similar testing strategy as for the Miller Rabin gives us sufficient confidence. 1. Test that we get small primes right. 2. Test that we get known pseudoprimes "correctly wrong". 3. Test some really big primes. (Remember, a probable prime test must mark every actual prime as "probably prime".) Helps #509. --- src/core/include/mp-units/ext/prime.h | 118 ++++++++++++++++++++++++++ test/static/prime_test.cpp | 32 +++++++ 2 files changed, 150 insertions(+) diff --git a/src/core/include/mp-units/ext/prime.h b/src/core/include/mp-units/ext/prime.h index 7d6f7c5b..4adc88e8 100644 --- a/src/core/include/mp-units/ext/prime.h +++ b/src/core/include/mp-units/ext/prime.h @@ -260,6 +260,124 @@ struct NumberDecomposition { } } +// The "D" parameter in the Strong Lucas Probable Prime test. +// +// Default construction produces the first value to try, according to Selfridge's parameter selection. +// Calling `successor()` on this repeatedly will produce the sequence of values to try. +struct LucasDParameter { + uint64_t mag = 5u; + bool pos = true; + + friend constexpr int as_int(LucasDParameter p) + { + int D = static_cast(p.mag); + return p.pos ? D : -D; + } + friend constexpr LucasDParameter successor(LucasDParameter p) { return {.mag = p.mag + 2u, .pos = !p.pos}; } +}; + +// The first `D` in the infinite sequence {5, -7, 9, -11, ...} whose Jacobi symbol is -1. This is the D we +// want to use for the Strong Lucas Probable Prime test. +// +// Precondition: `n` is not a perfect square. +[[nodiscard]] consteval LucasDParameter find_first_D_with_jacobi_symbol_neg_one(uint64_t n) +{ + LucasDParameter D{}; + while (jacobi_symbol(as_int(D), n) != -1) { + D = successor(D); + } + return D; +} + +// An element of the Lucas sequence, with parameters tuned for the Strong Lucas test. +// +// The default values give the first element (i.e., k=1) of the sequence. +struct LucasSequenceElement { + uint64_t u = 1u; + uint64_t v = 1u; +}; + +// Produce the Lucas element whose index is twice the input element's index. +[[nodiscard]] consteval LucasSequenceElement double_strong_lucas_index(const LucasSequenceElement& element, uint64_t n, + LucasDParameter D) +{ + const auto& [u, v] = element; + + uint64_t v_squared = mul_mod(v, v, n); + uint64_t D_u_squared = mul_mod(D.mag, mul_mod(u, u, n), n); + uint64_t new_v = D.pos ? add_mod(v_squared, D_u_squared, n) : sub_mod(v_squared, D_u_squared, n); + new_v = half_mod_odd(new_v, n); + + return {.u = mul_mod(u, v, n), .v = new_v}; +} + +[[nodiscard]] consteval LucasSequenceElement increment_strong_lucas_index(const LucasSequenceElement& element, + uint64_t n, LucasDParameter D) +{ + const auto& [u, v] = element; + + const auto new_u = half_mod_odd(add_mod(u, v, n), n); + + const auto D_u = mul_mod(D.mag, u, n); + auto new_v = D.pos ? add_mod(v, D_u, n) : sub_mod(v, D_u, n); + new_v = half_mod_odd(new_v, n); + + return {.u = new_u, .v = new_v}; +} + +[[nodiscard]] consteval LucasSequenceElement find_strong_lucas_element(uint64_t i, uint64_t n, LucasDParameter D) +{ + LucasSequenceElement element{}; + + bool bits[64] = {}; + std::size_t n_bits = 0u; + while (i > 1u) { + bits[n_bits++] = (i & 1u); + i >>= 1u; + } + + for (auto j = n_bits; j > 0u; --j) { + element = double_strong_lucas_index(element, n, D); + if (bits[j - 1u]) { + element = increment_strong_lucas_index(element, n, D); + } + } + + return element; +} + +// Perform a Strong Lucas Probable Prime test on `n`. +// +// Precondition: (n >= 2). +// Precondition: (n is odd). +[[nodiscard]] consteval bool strong_lucas_probable_prime(uint64_t n) +{ + MP_UNITS_EXPECTS_DEBUG(n >= 2u); + MP_UNITS_EXPECTS_DEBUG(n % 2u == 1u); + + if (is_perfect_square(n)) { + return false; + } + + const auto D = find_first_D_with_jacobi_symbol_neg_one(n); + + const auto [s, d] = decompose(n + 1u); + + auto element = find_strong_lucas_element(d, n, D); + if (element.u == 0u) { + return true; + } + + for (auto i = 0u; i < s; ++i) { + if (element.v == 0u) { + return true; + } + element = double_strong_lucas_index(element, n, D); + } + + return false; +} + [[nodiscard]] consteval bool is_prime_by_trial_division(std::uintmax_t n) { for (std::uintmax_t f = 2; f * f <= n; f += 1 + (f % 2)) { diff --git a/test/static/prime_test.cpp b/test/static/prime_test.cpp index 5b279527..5378e211 100644 --- a/test/static/prime_test.cpp +++ b/test/static/prime_test.cpp @@ -165,4 +165,36 @@ static_assert(!is_perfect_square(BIG_SQUARE - 1u)); static_assert(is_perfect_square(BIG_SQUARE)); static_assert(!is_perfect_square(BIG_SQUARE + 1u)); +// Tests for the Strong Lucas Probable Prime test. +static_assert(as_int(LucasDParameter{.mag = 5, .pos = true}) == 5); +static_assert(as_int(LucasDParameter{.mag = 7, .pos = false}) == -7); + +static_assert(as_int(LucasDParameter{}) == 5, "First D parameter in the sequence is 5"); +static_assert(as_int(successor(LucasDParameter{})) == -7, "Incrementing adds 2 to the mag, and flips the sign"); +static_assert(as_int(successor(successor(LucasDParameter{}))) == 9); +static_assert(as_int(successor(successor(successor(LucasDParameter{})))) == -11); + +static_assert(strong_lucas_probable_prime(3u), "Known small prime"); +static_assert(strong_lucas_probable_prime(5u), "Known small prime"); +static_assert(strong_lucas_probable_prime(7u), "Known small prime"); +static_assert(!strong_lucas_probable_prime(9u), "Known small composite"); + +// Test some Miller-Rabin pseudoprimes (https://oeis.org/A001262), which should NOT be marked prime. +static_assert(!strong_lucas_probable_prime(2047u), "Miller-Rabin pseudoprime"); +static_assert(!strong_lucas_probable_prime(3277u), "Miller-Rabin pseudoprime"); +static_assert(!strong_lucas_probable_prime(486737u), "Miller-Rabin pseudoprime"); + +// Test some Strong Lucas pseudoprimes (https://oeis.org/A217255). +static_assert(strong_lucas_probable_prime(5459u), "Strong Lucas pseudoprime"); +static_assert(strong_lucas_probable_prime(5777u), "Strong Lucas pseudoprime"); +static_assert(strong_lucas_probable_prime(10877u), "Strong Lucas pseudoprime"); +static_assert(strong_lucas_probable_prime(324899u), "Strong Lucas pseudoprime"); + +// Test some actual primes +static_assert(strong_lucas_probable_prime(225'653'407'801u), "Large known prime"); +static_assert(strong_lucas_probable_prime(334'524'384'739u), "Large known prime"); +static_assert(strong_lucas_probable_prime(9'007'199'254'740'881u), "Large known prime"); + +static_assert(strong_lucas_probable_prime(18'446'744'073'709'551'557u), "Largest 64-bit prime"); + } // namespace From 346a0012832868533632530e13c82c06b6a7b04b Mon Sep 17 00:00:00 2001 From: Chip Hogg Date: Fri, 15 Nov 2024 10:56:53 -0500 Subject: [PATCH 05/10] Add more missing `std::` prefixes --- src/core/include/mp-units/ext/prime.h | 25 +++++++++++++------------ 1 file changed, 13 insertions(+), 12 deletions(-) diff --git a/src/core/include/mp-units/ext/prime.h b/src/core/include/mp-units/ext/prime.h index 4adc88e8..56b83a44 100644 --- a/src/core/include/mp-units/ext/prime.h +++ b/src/core/include/mp-units/ext/prime.h @@ -265,7 +265,7 @@ struct NumberDecomposition { // Default construction produces the first value to try, according to Selfridge's parameter selection. // Calling `successor()` on this repeatedly will produce the sequence of values to try. struct LucasDParameter { - uint64_t mag = 5u; + std::uint64_t mag = 5u; bool pos = true; friend constexpr int as_int(LucasDParameter p) @@ -280,7 +280,7 @@ struct LucasDParameter { // want to use for the Strong Lucas Probable Prime test. // // Precondition: `n` is not a perfect square. -[[nodiscard]] consteval LucasDParameter find_first_D_with_jacobi_symbol_neg_one(uint64_t n) +[[nodiscard]] consteval LucasDParameter find_first_D_with_jacobi_symbol_neg_one(std::uint64_t n) { LucasDParameter D{}; while (jacobi_symbol(as_int(D), n) != -1) { @@ -293,26 +293,26 @@ struct LucasDParameter { // // The default values give the first element (i.e., k=1) of the sequence. struct LucasSequenceElement { - uint64_t u = 1u; - uint64_t v = 1u; + std::uint64_t u = 1u; + std::uint64_t v = 1u; }; // Produce the Lucas element whose index is twice the input element's index. -[[nodiscard]] consteval LucasSequenceElement double_strong_lucas_index(const LucasSequenceElement& element, uint64_t n, - LucasDParameter D) +[[nodiscard]] consteval LucasSequenceElement double_strong_lucas_index(const LucasSequenceElement& element, + std::uint64_t n, LucasDParameter D) { const auto& [u, v] = element; - uint64_t v_squared = mul_mod(v, v, n); - uint64_t D_u_squared = mul_mod(D.mag, mul_mod(u, u, n), n); - uint64_t new_v = D.pos ? add_mod(v_squared, D_u_squared, n) : sub_mod(v_squared, D_u_squared, n); + std::uint64_t v_squared = mul_mod(v, v, n); + std::uint64_t D_u_squared = mul_mod(D.mag, mul_mod(u, u, n), n); + std::uint64_t new_v = D.pos ? add_mod(v_squared, D_u_squared, n) : sub_mod(v_squared, D_u_squared, n); new_v = half_mod_odd(new_v, n); return {.u = mul_mod(u, v, n), .v = new_v}; } [[nodiscard]] consteval LucasSequenceElement increment_strong_lucas_index(const LucasSequenceElement& element, - uint64_t n, LucasDParameter D) + std::uint64_t n, LucasDParameter D) { const auto& [u, v] = element; @@ -325,7 +325,8 @@ struct LucasSequenceElement { return {.u = new_u, .v = new_v}; } -[[nodiscard]] consteval LucasSequenceElement find_strong_lucas_element(uint64_t i, uint64_t n, LucasDParameter D) +[[nodiscard]] consteval LucasSequenceElement find_strong_lucas_element(std::uint64_t i, std::uint64_t n, + LucasDParameter D) { LucasSequenceElement element{}; @@ -350,7 +351,7 @@ struct LucasSequenceElement { // // Precondition: (n >= 2). // Precondition: (n is odd). -[[nodiscard]] consteval bool strong_lucas_probable_prime(uint64_t n) +[[nodiscard]] consteval bool strong_lucas_probable_prime(std::uint64_t n) { MP_UNITS_EXPECTS_DEBUG(n >= 2u); MP_UNITS_EXPECTS_DEBUG(n % 2u == 1u); From d963ce8f936aa18ecc2a10795c86cc25c92d9f49 Mon Sep 17 00:00:00 2001 From: Chip Hogg Date: Fri, 15 Nov 2024 11:58:49 -0500 Subject: [PATCH 06/10] Apply mod to D Looks like we were violating a precondition. --- src/core/include/mp-units/ext/prime.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/core/include/mp-units/ext/prime.h b/src/core/include/mp-units/ext/prime.h index 56b83a44..918d1967 100644 --- a/src/core/include/mp-units/ext/prime.h +++ b/src/core/include/mp-units/ext/prime.h @@ -304,7 +304,7 @@ struct LucasSequenceElement { const auto& [u, v] = element; std::uint64_t v_squared = mul_mod(v, v, n); - std::uint64_t D_u_squared = mul_mod(D.mag, mul_mod(u, u, n), n); + std::uint64_t D_u_squared = mul_mod(D.mag % n, mul_mod(u, u, n), n); std::uint64_t new_v = D.pos ? add_mod(v_squared, D_u_squared, n) : sub_mod(v_squared, D_u_squared, n); new_v = half_mod_odd(new_v, n); @@ -318,7 +318,7 @@ struct LucasSequenceElement { const auto new_u = half_mod_odd(add_mod(u, v, n), n); - const auto D_u = mul_mod(D.mag, u, n); + const auto D_u = mul_mod(D.mag % n, u, n); auto new_v = D.pos ? add_mod(v, D_u, n) : sub_mod(v, D_u, n); new_v = half_mod_odd(new_v, n); From 82040c1788236e72674f54a8c74acdde71fd605c Mon Sep 17 00:00:00 2001 From: Chip Hogg Date: Fri, 15 Nov 2024 12:58:22 -0500 Subject: [PATCH 07/10] Handle edge case Subtracting from `n` doesn't work in literally every case: it fails for `0` input. --- src/core/include/mp-units/ext/prime.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/core/include/mp-units/ext/prime.h b/src/core/include/mp-units/ext/prime.h index 918d1967..3c3dc22d 100644 --- a/src/core/include/mp-units/ext/prime.h +++ b/src/core/include/mp-units/ext/prime.h @@ -100,7 +100,7 @@ namespace mp_units::detail { return add_mod( // Transform into "negative space" to make the first parameter as small as possible; // then, transform back. - n - mul_mod(n % a, num_batches, n), + (n - mul_mod(n % a, num_batches, n)) % n, // Handle the leftover product (which is guaranteed to fit in the integer type). (a * (b % batch_size)) % n, From 3b586a685f52cfacc2f35c964da42acb045e35d9 Mon Sep 17 00:00:00 2001 From: Chip Hogg Date: Thu, 14 Nov 2024 21:27:09 -0500 Subject: [PATCH 08/10] Implement Baillie-PSW --- src/core/include/mp-units/ext/prime.h | 22 ++++++++++++++++++++++ test/static/prime_test.cpp | 24 ++++++++++++++++++++++++ 2 files changed, 46 insertions(+) diff --git a/src/core/include/mp-units/ext/prime.h b/src/core/include/mp-units/ext/prime.h index 3c3dc22d..2b03a5ec 100644 --- a/src/core/include/mp-units/ext/prime.h +++ b/src/core/include/mp-units/ext/prime.h @@ -379,6 +379,28 @@ struct LucasSequenceElement { return false; } +// The Baillie-PSW test is technically a "probable prime" test. However, it is known to be correct for all +// 64-bit integers, and no counterexample of any size has ever been found. Thus, for this library's purposes, +// we can treat it as deterministic... probably. +[[nodiscard]] consteval bool baillie_psw_probable_prime(uint64_t n) +{ + if (n < 2u) { + return false; + } + if (n < 4u) { + return true; + } + if (n % 2u == 0u) { + return false; + } + + if (!miller_rabin_probable_prime(2u, n)) { + return false; + } + + return strong_lucas_probable_prime(n); +} + [[nodiscard]] consteval bool is_prime_by_trial_division(std::uintmax_t n) { for (std::uintmax_t f = 2; f * f <= n; f += 1 + (f % 2)) { diff --git a/test/static/prime_test.cpp b/test/static/prime_test.cpp index 5378e211..c7c317fa 100644 --- a/test/static/prime_test.cpp +++ b/test/static/prime_test.cpp @@ -197,4 +197,28 @@ static_assert(strong_lucas_probable_prime(9'007'199'254'740'881u), "Large known static_assert(strong_lucas_probable_prime(18'446'744'073'709'551'557u), "Largest 64-bit prime"); +// Tests for Baillie-PSW, which is known to be correct for all 64-bit integers. +static_assert(baillie_psw_probable_prime(3u), "Known small prime"); +static_assert(baillie_psw_probable_prime(5u), "Known small prime"); +static_assert(baillie_psw_probable_prime(7u), "Known small prime"); +static_assert(!baillie_psw_probable_prime(9u), "Known small composite"); + +// Test some Miller-Rabin pseudoprimes (https://oeis.org/A001262), which should NOT be marked prime. +static_assert(!baillie_psw_probable_prime(2047u), "Miller-Rabin pseudoprime"); +static_assert(!baillie_psw_probable_prime(3277u), "Miller-Rabin pseudoprime"); +static_assert(!baillie_psw_probable_prime(486737u), "Miller-Rabin pseudoprime"); + +// Test some Strong Lucas pseudoprimes (https://oeis.org/A217255), which should NOT be marked prime. +static_assert(!baillie_psw_probable_prime(5459u), "Strong Lucas pseudoprime"); +static_assert(!baillie_psw_probable_prime(5777u), "Strong Lucas pseudoprime"); +static_assert(!baillie_psw_probable_prime(10877u), "Strong Lucas pseudoprime"); +static_assert(!baillie_psw_probable_prime(324899u), "Strong Lucas pseudoprime"); + +// Test some actual primes +static_assert(baillie_psw_probable_prime(225'653'407'801u), "Large known prime"); +static_assert(baillie_psw_probable_prime(334'524'384'739u), "Large known prime"); +static_assert(baillie_psw_probable_prime(9'007'199'254'740'881u), "Large known prime"); + +static_assert(baillie_psw_probable_prime(18'446'744'073'709'551'557u), "Largest 64-bit prime"); + } // namespace From de68b3a9f1ff40304eef3c3545f35aca8e605409 Mon Sep 17 00:00:00 2001 From: Chip Hogg Date: Fri, 15 Nov 2024 06:30:16 -0500 Subject: [PATCH 09/10] Replace old factoring with Baillie-PSW --- src/core/include/mp-units/ext/prime.h | 117 ++++-------------- .../include/mp-units/framework/magnitude.h | 24 +--- src/systems/include/mp-units/systems/hep.h | 3 - test/static/magnitude_test.cpp | 17 --- test/static/prime_test.cpp | 45 ------- 5 files changed, 29 insertions(+), 177 deletions(-) diff --git a/src/core/include/mp-units/ext/prime.h b/src/core/include/mp-units/ext/prime.h index 2b03a5ec..61346ce9 100644 --- a/src/core/include/mp-units/ext/prime.h +++ b/src/core/include/mp-units/ext/prime.h @@ -401,30 +401,6 @@ struct LucasSequenceElement { return strong_lucas_probable_prime(n); } -[[nodiscard]] consteval bool is_prime_by_trial_division(std::uintmax_t n) -{ - for (std::uintmax_t f = 2; f * f <= n; f += 1 + (f % 2)) { - if (n % f == 0) { - return false; - } - } - return true; -} - -// Return the first factor of n, as long as it is either k or n. -// -// Precondition: no integer smaller than k evenly divides n. -[[nodiscard]] constexpr std::optional first_factor_maybe(std::uintmax_t n, std::uintmax_t k) -{ - if (n % k == 0) { - return k; - } - if (k * k > n) { - return n; - } - return std::nullopt; -} - template [[nodiscard]] consteval std::array first_n_primes() { @@ -432,42 +408,13 @@ template primes[0] = 2; for (std::size_t i = 1; i < N; ++i) { primes[i] = primes[i - 1] + 1; - while (!is_prime_by_trial_division(primes[i])) { + while (!baillie_psw_probable_prime(primes[i])) { ++primes[i]; } } return primes; } -template -consteval void call_for_coprimes_up_to(std::uintmax_t n, const std::array& basis, Callable call) -{ - for (std::uintmax_t i = 0u; i < n; ++i) { - if (std::apply([&i](auto... primes) { return ((i % primes != 0) && ...); }, basis)) { - call(i); - } - } -} - -template -[[nodiscard]] consteval std::size_t num_coprimes_up_to(std::uintmax_t n, const std::array& basis) -{ - std::size_t count = 0u; - call_for_coprimes_up_to(n, basis, [&count](auto) { ++count; }); - return count; -} - -template -[[nodiscard]] consteval auto coprimes_up_to(std::size_t n, const std::array& basis) -{ - std::array coprimes{}; - std::size_t i = 0u; - - call_for_coprimes_up_to(n, basis, [&coprimes, &i](std::uintmax_t cp) { coprimes[i++] = cp; }); - - return coprimes; -} - template [[nodiscard]] consteval std::uintmax_t product(const std::array& values) { @@ -494,48 +441,34 @@ constexpr auto get_first_of(const Rng& rng, UnaryFunction f) return get_first_of(begin(rng), end(rng), f); } -// A configurable instantiation of the "wheel factorization" algorithm [1] for prime numbers. -// -// Instantiate with N to use a "basis" of the first N prime numbers. Higher values of N use fewer trial divisions, at -// the cost of additional space. The amount of space consumed is roughly the total number of numbers that are a) less -// than the _product_ of the basis elements (first N primes), and b) coprime with every element of the basis. This -// means it grows rapidly with N. Consider this approximate chart: -// -// N | Num coprimes | Trial divisions needed -// --+--------------+----------------------- -// 1 | 1 | 50.0 % -// 2 | 2 | 33.3 % -// 3 | 8 | 26.7 % -// 4 | 48 | 22.9 % -// 5 | 480 | 20.8 % -// -// Note the diminishing returns, and the rapidly escalating costs. Consider this behaviour when choosing the value of N -// most appropriate for your needs. -// -// [1] https://en.wikipedia.org/wiki/Wheel_factorization -template -struct wheel_factorizer { - static constexpr auto basis = first_n_primes(); - static constexpr auto wheel_size = product(basis); - static constexpr auto coprimes_in_first_wheel = - coprimes_up_to(wheel_size, basis); +[[nodiscard]] consteval std::uintmax_t find_first_factor(std::uintmax_t n) +{ + constexpr auto first_100_primes = first_n_primes<100>(); - [[nodiscard]] static consteval std::uintmax_t find_first_factor(std::uintmax_t n) - { - if (const auto k = detail::get_first_of(basis, [&](auto p) { return first_factor_maybe(n, p); })) return *k; + for (const auto& p : first_100_primes) { + if (n % p == 0u) { + return p; + } + if (p * p > n) { + return n; + } + } - if (const auto k = detail::get_first_of(std::next(begin(coprimes_in_first_wheel)), end(coprimes_in_first_wheel), - [&](auto p) { return first_factor_maybe(n, p); })) - return *k; - - for (std::size_t wheel = wheel_size; wheel < n; wheel += wheel_size) - if (const auto k = - detail::get_first_of(coprimes_in_first_wheel, [&](auto p) { return first_factor_maybe(n, wheel + p); })) - return *k; + // If we got this far, and we haven't found a factor, and haven't terminated... maybe it's prime? Do a fast check. + if (baillie_psw_probable_prime(n)) { return n; } - [[nodiscard]] static consteval bool is_prime(std::size_t n) { return (n > 1) && find_first_factor(n) == n; } -}; + // If we're here, we know `n` is composite, so continue with trial division for all odd numbers. + std::uintmax_t factor = first_100_primes.back() + 2u; + while (factor * factor <= n) { + if (n % factor == 0u) { + return factor; + } + factor += 2u; + } + + return n; // Technically unreachable. +} } // namespace mp_units::detail diff --git a/src/core/include/mp-units/framework/magnitude.h b/src/core/include/mp-units/framework/magnitude.h index 6ab5a1d5..cdb796f9 100644 --- a/src/core/include/mp-units/framework/magnitude.h +++ b/src/core/include/mp-units/framework/magnitude.h @@ -51,13 +51,6 @@ import std; namespace mp_units { -namespace detail { - -// Higher numbers use fewer trial divisions, at the price of more storage space. -using factorizer = wheel_factorizer<4>; - -} // namespace detail - #if defined MP_UNITS_COMP_CLANG || MP_UNITS_COMP_CLANG < 18 MP_UNITS_EXPORT template @@ -629,14 +622,10 @@ using common_magnitude_type = decltype(_common_magnitude_type_impl(M)); //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// // `mag()` implementation. -// Sometimes we need to give the compiler a "shortcut" when factorizing large numbers (specifically, numbers whose -// _first factor_ is very large). If we don't, we can run into limits on the number of constexpr steps or iterations. -// -// To provide the first factor for a given number, specialize this variable template. -// -// WARNING: The program behaviour will be undefined if you provide a wrong answer, so check your math! MP_UNITS_EXPORT template -constexpr std::optional known_first_factor = std::nullopt; +[[deprecated("`known_first_factor` is no longer necessary and can simply be removed")]] +constexpr std::optional + known_first_factor = std::nullopt; namespace detail { @@ -646,12 +635,7 @@ template struct prime_factorization { [[nodiscard]] static consteval std::intmax_t get_or_compute_first_factor() { - constexpr auto opt = known_first_factor; - if constexpr (opt.has_value()) { - return opt.value(); // NOLINT(bugprone-unchecked-optional-access) - } else { - return static_cast(factorizer::find_first_factor(N)); - } + return static_cast(find_first_factor(N)); } static constexpr std::intmax_t first_base = get_or_compute_first_factor(); diff --git a/src/systems/include/mp-units/systems/hep.h b/src/systems/include/mp-units/systems/hep.h index c64844f8..b4bb997e 100644 --- a/src/systems/include/mp-units/systems/hep.h +++ b/src/systems/include/mp-units/systems/hep.h @@ -35,9 +35,6 @@ MP_UNITS_EXPORT namespace mp_units { -template<> -constexpr std::optional known_first_factor<334'524'384'739> = 334'524'384'739; - namespace hep { // energy diff --git a/test/static/magnitude_test.cpp b/test/static/magnitude_test.cpp index 146e709a..72bc8af3 100644 --- a/test/static/magnitude_test.cpp +++ b/test/static/magnitude_test.cpp @@ -31,9 +31,6 @@ import mp_units; using namespace units; using namespace units::detail; -template<> -constexpr std::optional units::known_first_factor<9223372036854775783> = 9223372036854775783; - namespace { // A set of non-standard bases for testing purposes. @@ -182,20 +179,6 @@ static_assert(std::is_same_v{})), mag_2_> // mag_ratio<16'605'390'666'050, 10'000'000'000'000>(); // } -// SECTION("Can bypass computing primes by providing known_first_factor") -// { -// // Sometimes, even wheel factorization isn't enough to handle the compilers' limits on constexpr steps and/or -// // iterations. To work around these cases, we can explicitly provide the correct answer directly to the -// compiler. -// // -// // In this case, we test that we can represent the largest prime that fits in a signed 64-bit int. The reason -// this -// // test can pass is that we have provided the answer, by specializing the `known_first_factor` variable template -// // above in this file. -// mag<9'223'372'036'854'775'783>(); -// } -// } - // TEST_CASE("magnitude converts to numerical value") // { // SECTION("Positive integer powers of integer bases give integer values") diff --git a/test/static/prime_test.cpp b/test/static/prime_test.cpp index c7c317fa..d62e49f7 100644 --- a/test/static/prime_test.cpp +++ b/test/static/prime_test.cpp @@ -35,51 +35,6 @@ namespace { inline constexpr auto MAX_U64 = std::numeric_limits::max(); -template -constexpr bool check_primes(std::index_sequence) -{ - return ((Is < 2 || wheel_factorizer::is_prime(Is) == is_prime_by_trial_division(Is)) && ...); -} - -static_assert(check_primes<2>(std::make_index_sequence<122>{})); - -// This is the smallest number that can catch the bug where we use only _prime_ numbers in the first wheel, rather than -// numbers which are _coprime to the basis_. -// -// The basis for N = 4 is {2, 3, 5, 7}, so the wheel size is 210. 11 * 11 = 121 is within the first wheel. It is -// coprime with every element of the basis, but it is _not_ prime. If we keep only prime numbers, then we will neglect -// using numbers of the form (210 * n + 121) as trial divisors, which is a problem if any are prime. For n = 1, we have -// a divisor of (210 + 121 = 331), which happens to be prime but will not be used. Thus, (331 * 331 = 109561) is a -// composite number which could wrongly appear prime if we skip over 331. -static_assert(wheel_factorizer<4>::is_prime(109'561) == is_prime_by_trial_division(109'561)); - -static_assert(wheel_factorizer<1>::coprimes_in_first_wheel.size() == 1); -static_assert(wheel_factorizer<2>::coprimes_in_first_wheel.size() == 2); -static_assert(wheel_factorizer<3>::coprimes_in_first_wheel.size() == 8); -static_assert(wheel_factorizer<4>::coprimes_in_first_wheel.size() == 48); -static_assert(wheel_factorizer<5>::coprimes_in_first_wheel.size() == 480); - -static_assert(wheel_factorizer<3>::coprimes_in_first_wheel[0] == 1); -static_assert(wheel_factorizer<3>::coprimes_in_first_wheel[1] == 7); -static_assert(wheel_factorizer<3>::coprimes_in_first_wheel[2] == 11); -static_assert(wheel_factorizer<3>::coprimes_in_first_wheel[3] == 13); -static_assert(wheel_factorizer<3>::coprimes_in_first_wheel[4] == 17); -static_assert(wheel_factorizer<3>::coprimes_in_first_wheel[5] == 19); -static_assert(wheel_factorizer<3>::coprimes_in_first_wheel[6] == 23); -static_assert(wheel_factorizer<3>::coprimes_in_first_wheel[7] == 29); - -static_assert(!wheel_factorizer<1>::is_prime(0)); -static_assert(!wheel_factorizer<1>::is_prime(1)); -static_assert(wheel_factorizer<1>::is_prime(2)); - -static_assert(!wheel_factorizer<2>::is_prime(0)); -static_assert(!wheel_factorizer<2>::is_prime(1)); -static_assert(wheel_factorizer<2>::is_prime(2)); - -static_assert(!wheel_factorizer<3>::is_prime(0)); -static_assert(!wheel_factorizer<3>::is_prime(1)); -static_assert(wheel_factorizer<3>::is_prime(2)); - // Modular arithmetic. static_assert(add_mod(1u, 2u, 5u) == 3u); static_assert(add_mod(4u, 4u, 5u) == 3u); From a93a5d6d2ef3b0f34792d2e72b9a6519de1458f0 Mon Sep 17 00:00:00 2001 From: Chip Hogg Date: Fri, 15 Nov 2024 13:38:48 -0500 Subject: [PATCH 10/10] Add `std::` prefix --- src/core/include/mp-units/ext/prime.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/core/include/mp-units/ext/prime.h b/src/core/include/mp-units/ext/prime.h index 61346ce9..1206dc7e 100644 --- a/src/core/include/mp-units/ext/prime.h +++ b/src/core/include/mp-units/ext/prime.h @@ -382,7 +382,7 @@ struct LucasSequenceElement { // The Baillie-PSW test is technically a "probable prime" test. However, it is known to be correct for all // 64-bit integers, and no counterexample of any size has ever been found. Thus, for this library's purposes, // we can treat it as deterministic... probably. -[[nodiscard]] consteval bool baillie_psw_probable_prime(uint64_t n) +[[nodiscard]] consteval bool baillie_psw_probable_prime(std::uint64_t n) { if (n < 2u) { return false;