mpusz/mp-units
1
0
Fork 1
mirror of https://github.com/mpusz/mp-units.git synced 2025-06-25 09:11:34 +02:00
Code Issues Packages Projects Releases Wiki Activity
Files
master
mp-units/test/static/prime_test.cpp

180 lines
9.3 KiB
C++
Raw Permalink Normal View History

Use wheel factorization for prime numbers Certain existing units in the library require very large prime numbers---so large, in fact, that our naive trial division hits the _iteration limit_ for `constexpr` loops. We don't want to force users to provide a compiler option override, so we'd better find another way. The solution is to use the "wheel factorization" algorithm: https://en.wikipedia.org/wiki/Wheel_factorization This lets us skip most composite numbers in our trial division. The implementation presented here is configurable in terms of the size of the "basis" of primes we use. Bigger bases let us skip more primes, but at the cost of storing more numbers. Fortunately, it turns out that N=3 was good enough for our purposes.
2022-03-10 23:29:02 +00:00
// The MIT License (MIT)
//
// Copyright (c) 2018 Mateusz Pusz
//
// Permission is hereby granted, free of charge, to any person obtaining a copy
// of this software and associated documentation files (the "Software"), to deal
// in the Software without restriction, including without limitation the rights
// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
// copies of the Software, and to permit persons to whom the Software is
// furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included in all
// copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
// SOFTWARE.
refactor: :boom: `framework.h` introduced + IWYU
2024-04-24 20:53:54 +02:00
#include <mp-units/ext/prime.h>
feat: `import std;` support added Resolves #595
2024-07-16 17:36:00 +02:00
#ifdef MP_UNITS_IMPORT_STD
import std;
#else
refactor: :boom: `framework.h` introduced + IWYU
2024-04-24 20:53:54 +02:00
#include <array>
#include <cstddef>
Use wheel factorization for prime numbers Certain existing units in the library require very large prime numbers---so large, in fact, that our naive trial division hits the _iteration limit_ for `constexpr` loops. We don't want to force users to provide a compiler option override, so we'd better find another way. The solution is to use the "wheel factorization" algorithm: https://en.wikipedia.org/wiki/Wheel_factorization This lets us skip most composite numbers in our trial division. The implementation presented here is configurable in terms of the size of the "basis" of primes we use. Bigger bases let us skip more primes, but at the cost of storing more numbers. Fortunately, it turns out that N=3 was good enough for our purposes.
2022-03-10 23:29:02 +00:00
#include <utility>
feat: `import std;` support added Resolves #595
2024-07-16 17:36:00 +02:00
#endif
Use wheel factorization for prime numbers Certain existing units in the library require very large prime numbers---so large, in fact, that our naive trial division hits the _iteration limit_ for `constexpr` loops. We don't want to force users to provide a compiler option override, so we'd better find another way. The solution is to use the "wheel factorization" algorithm: https://en.wikipedia.org/wiki/Wheel_factorization This lets us skip most composite numbers in our trial division. The implementation presented here is configurable in terms of the size of the "basis" of primes we use. Bigger bases let us skip more primes, but at the cost of storing more numbers. Fortunately, it turns out that N=3 was good enough for our purposes.
2022-03-10 23:29:02 +00:00
feat: preventing possible future `units` crash when `using namespace std;` Resolves #317
2022-12-29 20:18:48 +01:00
using namespace mp_units::detail;
Switch tests to use top-level, anonymous namespace
2022-03-19 15:54:29 +00:00
namespace {
Use wheel factorization for prime numbers Certain existing units in the library require very large prime numbers---so large, in fact, that our naive trial division hits the _iteration limit_ for `constexpr` loops. We don't want to force users to provide a compiler option override, so we'd better find another way. The solution is to use the "wheel factorization" algorithm: https://en.wikipedia.org/wiki/Wheel_factorization This lets us skip most composite numbers in our trial division. The implementation presented here is configurable in terms of the size of the "basis" of primes we use. Bigger bases let us skip more primes, but at the cost of storing more numbers. Fortunately, it turns out that N=3 was good enough for our purposes.
2022-03-10 23:29:02 +00:00
Add helpers for modular arithmetic The prime-testing techniques we will use (Miller-Rabin, Strong Lucas) all make heavy usage of modular arithmetic. Therefore, we lay those foundations here, adding utilities to perform the basic arithmetic operations robustly. Since these are internal-only helper functions, we don't bother checking the preconditions, although we state them clearly in the contract comment for each utility. After C++26, we could add contracts for these. Helps #509.
2024-11-11 12:55:41 -05:00
inline constexpr auto MAX_U64 = std::numeric_limits<std::uint64_t>::max();
// Modular arithmetic.
static_assert(add_mod(1u, 2u, 5u) == 3u);
static_assert(add_mod(4u, 4u, 5u) == 3u);
static_assert(add_mod(MAX_U64 - 1u, MAX_U64 - 2u, MAX_U64) == MAX_U64 - 3u);
static_assert(sub_mod(2u, 1u, 5u) == 1u);
static_assert(sub_mod(1u, 2u, 5u) == 4u);
static_assert(sub_mod(MAX_U64 - 2u, MAX_U64 - 1u, MAX_U64) == MAX_U64 - 1u);
static_assert(sub_mod(1u, MAX_U64 - 1u, MAX_U64) == 2u);
static_assert(mul_mod(6u, 7u, 10u) == 2u);
static_assert(mul_mod(13u, 11u, 50u) == 43u);
static_assert(mul_mod(MAX_U64 / 2u, 10u, MAX_U64) == MAX_U64 - 5u);
static_assert(half_mod_odd(0u, 11u) == 0u);
static_assert(half_mod_odd(10u, 11u) == 5u);
static_assert(half_mod_odd(1u, 11u) == 6u);
static_assert(half_mod_odd(9u, 11u) == 10u);
static_assert(half_mod_odd(MAX_U64 - 1u, MAX_U64) == (MAX_U64 - 1u) / 2u);
static_assert(half_mod_odd(MAX_U64 - 2u, MAX_U64) == MAX_U64 - 1u);
static_assert(pow_mod(5u, 8u, 9u) == ((5u * 5u * 5u * 5u) * (5u * 5u * 5u * 5u)) % 9u);
static_assert(pow_mod(2u, 64u, MAX_U64) == 1u);
Add Miller-Rabin probable prime test This can mark a number as either "probably prime", or "definitely composite". The first parameter is the base, and the second is the number to test. Future PRs will build up the Strong Lucas test which complements this, and then form the Baillie-PSW test by composing the two. Helps #506.
2024-11-13 19:56:05 -05:00
// Miller-Rabin primality testing.
static_assert(miller_rabin_probable_prime(2u, 5u));
static_assert(miller_rabin_probable_prime(2u, 7u));
static_assert(!miller_rabin_probable_prime(2u, 9u));
static_assert(miller_rabin_probable_prime(2u, 11u));
static_assert(miller_rabin_probable_prime(2u, 2047u), "Known base 2 pseudoprime");
static_assert(miller_rabin_probable_prime(2u, 3277u), "Known base 2 pseudoprime");
static_assert(miller_rabin_probable_prime(3u, 121u), "Known base 3 pseudoprime");
static_assert(miller_rabin_probable_prime(3u, 703u), "Known base 3 pseudoprime");
static_assert(miller_rabin_probable_prime(2u, 225'653'407'801u), "Large known prime");
static_assert(miller_rabin_probable_prime(2u, 334'524'384'739u), "Large known prime");
static_assert(miller_rabin_probable_prime(2u, 9'007'199'254'740'881u), "Large known prime");
static_assert(miller_rabin_probable_prime(2u, 18'446'744'073'709'551'557u), "Largest 64-bit prime");
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.
2024-11-14 18:50:21 -05:00
// 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));
Add some more
2024-11-15 09:39:39 -05:00
constexpr std::uint64_t BIG_SQUARE = [](auto x) { return x * x; }((std::uint64_t{1u} << 32) - 1u);
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.
2024-11-14 18:50:21 -05:00
static_assert(!is_perfect_square(BIG_SQUARE - 1u));
static_assert(is_perfect_square(BIG_SQUARE));
static_assert(!is_perfect_square(BIG_SQUARE + 1u));
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.
2024-11-14 20:54:31 -05:00
// 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");
Implement Baillie-PSW
2024-11-14 21:27:09 -05:00
// 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");
Switch tests to use top-level, anonymous namespace
2022-03-19 15:54:29 +00:00
} // namespace
Reference in New Issue Copy Permalink