diff --git a/src/core/include/mp-units/ext/prime.h b/src/core/include/mp-units/ext/prime.h index 51050f41..1206dc7e 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, @@ -183,28 +183,222 @@ struct NumberDecomposition { return false; } -[[nodiscard]] consteval bool is_prime_by_trial_division(std::uintmax_t n) +// 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(std::int64_t raw_a, std::uint64_t n) { - for (std::uintmax_t f = 2; f * f <= n; f += 1 + (f % 2)) { - if (n % f == 0) { - return false; - } + // Rule 1: n=1 case. + if (n == 1u) { + return 1; } - return true; + + // 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 std::uint64_t new_a = n % a; + n = a; + a = new_a; + } + + return 0; } -// 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) +[[nodiscard]] consteval bool is_perfect_square(std::uint64_t n) { - if (n % k == 0) { - return k; + if (n < 2u) { + return true; } - if (k * k > n) { - return n; + + std::uint64_t prev = n / 2u; + while (true) { + const std::uint64_t curr = (prev + n / prev) / 2u; + if (curr * curr == n) { + return true; + } + if (curr >= prev) { + return false; + } + prev = curr; } - return std::nullopt; +} + +// 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 { + std::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(std::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 { + 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, + std::uint64_t n, LucasDParameter D) +{ + 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 % 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); + + return {.u = mul_mod(u, v, n), .v = new_v}; +} + +[[nodiscard]] consteval LucasSequenceElement increment_strong_lucas_index(const LucasSequenceElement& element, + std::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 % 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); + + return {.u = new_u, .v = new_v}; +} + +[[nodiscard]] consteval LucasSequenceElement find_strong_lucas_element(std::uint64_t i, std::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(std::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; +} + +// 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(std::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); } template @@ -214,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) { @@ -276,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 bfd7aae5..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); @@ -122,4 +77,103 @@ 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 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)); + +// 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"); + +// 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