diff --git a/src/core/include/units/bits/prime.h b/src/core/include/units/bits/prime.h new file mode 100644 index 00000000..b9e06e72 --- /dev/null +++ b/src/core/include/units/bits/prime.h @@ -0,0 +1,165 @@ +// 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. + +#pragma once + +#include +#include +#include +#include +#include +#include + +namespace units::detail { + +constexpr bool is_prime_by_trial_division(std::size_t n) +{ + for (std::size_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. +constexpr std::optional first_factor_maybe(std::size_t n, std::size_t k) +{ + if (n % k == 0) { + return k; + } + if (k * k > n) { + return n; + } + return std::nullopt; +} + +template +constexpr std::array first_n_primes() +{ + std::array primes; + 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])) { + ++primes[i]; + } + } + return primes; +} + +template +constexpr void call_for_coprimes_up_to(std::size_t n, const std::array& basis, Callable&& call) +{ + for (std::size_t i = 0u; i < n; ++i) { + if (std::apply([&i](auto... primes) { return ((i % primes != 0) && ...); }, basis)) { + call(i); + } + } +} + +template +constexpr std::size_t num_coprimes_up_to(std::size_t n, const std::array& basis) +{ + std::size_t count = 0u; + call_for_coprimes_up_to(n, basis, [&count](auto) { ++count; }); + return count; +} + +template +constexpr 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::size_t cp) { coprimes[i++] = cp; }); + + return coprimes; +} + +template +constexpr std::size_t product(const std::array& values) +{ + std::size_t product = 1; + for (const auto& v : values) { + product *= v; + } + return product; +} + +// 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 std::size_t wheel_size = product(basis); + static constexpr auto coprimes_in_first_wheel = + coprimes_up_to(wheel_size, basis); + + static constexpr std::size_t find_first_factor(std::size_t n) + { + for (const auto& p : basis) { + if (const auto k = first_factor_maybe(n, p)) { + return *k; + } + } + + for (auto it = std::next(std::begin(coprimes_in_first_wheel)); it != std::end(coprimes_in_first_wheel); ++it) { + if (const auto k = first_factor_maybe(n, *it)) { + return *k; + } + } + + for (std::size_t wheel = wheel_size; wheel < n; wheel += wheel_size) { + for (const auto& p : coprimes_in_first_wheel) { + if (const auto k = first_factor_maybe(n, wheel + p)) { + return *k; + } + } + } + + return n; + } + + static constexpr bool is_prime(std::size_t n) { return (n > 1) && find_first_factor(n) == n; } +}; + +} // namespace units::detail diff --git a/src/core/include/units/magnitude.h b/src/core/include/units/magnitude.h index 05000e9f..94140397 100644 --- a/src/core/include/units/magnitude.h +++ b/src/core/include/units/magnitude.h @@ -23,12 +23,19 @@ #pragma once #include +#include #include +#include #include #include +#include #include namespace units { +namespace detail { +// Higher numbers use fewer trial divisions, at the price of more storage space. +using factorizer = wheel_factorizer<4>; +} // namespace detail /** * @brief Any type which can be used as a basis vector in a BasePower. @@ -244,17 +251,6 @@ constexpr auto pow(BasePower auto bp, ratio p) // A variety of implementation detail helpers. namespace detail { -// Find the smallest prime factor of `n`. -constexpr std::intmax_t find_first_factor(std::intmax_t n) -{ - for (std::intmax_t f = 2; f * f <= n; f += 1 + (f % 2)) { - if (n % f == 0) { - return f; - } - } - return n; -} - // The exponent of `factor` in the prime factorization of `n`. constexpr std::intmax_t multiplicity(std::intmax_t factor, std::intmax_t n) { @@ -278,7 +274,7 @@ constexpr std::intmax_t remove_power(std::intmax_t base, std::intmax_t pow, std: } // A way to check whether a number is prime at compile time. -constexpr bool is_prime(std::intmax_t n) { return (n > 1) && (find_first_factor(n) == n); } +constexpr bool is_prime(std::intmax_t n) { return (n >= 0) && factorizer::is_prime(static_cast(n)); } constexpr bool is_valid_base_power(const BasePower auto& bp) { @@ -287,6 +283,20 @@ constexpr bool is_valid_base_power(const BasePower auto& bp) } if constexpr (std::is_same_v) { + // Some prime numbers are so big, that we can't check their primality without exhausting limits on constexpr steps + // and/or iterations. We can still _perform_ the factorization for these by using the `known_first_factor` + // workaround. However, we can't _check_ that they are prime, because this workaround depends on the input being + // usable in a constexpr expression. This is true for `prime_factorization` (below), where the input `N` is a + // template parameter, but is not true for our case, where the input `bp.get_base()` is a function parameter. (See + // http://www.open-std.org/jtc1/sc22/wg21/docs/papers/2019/p1045r1.html for some background reading on this + // distinction.) + // + // In our case: we simply give up on excluding every possible ill-formed base power, and settle for catching the + // most likely and common mistakes. + if (const bool too_big_to_check = (bp.get_base() > 1'000'000'000)) { + return true; + } + return is_prime(bp.get_base()); } else { return bp.get_base() > 0; @@ -473,12 +483,30 @@ constexpr auto operator/(Magnitude auto l, Magnitude auto r) { return l * pow<-1 //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// // `as_magnitude()` 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! +template +inline constexpr std::optional known_first_factor = std::nullopt; + namespace detail { // Helper to perform prime factorization at compile time. template requires(N > 0) struct prime_factorization { - static constexpr std::intmax_t first_base = find_first_factor(N); + static constexpr std::intmax_t get_or_compute_first_factor() + { + if constexpr (known_first_factor.has_value()) { + return known_first_factor.value(); + } else { + return static_cast(factorizer::find_first_factor(N)); + } + } + + static constexpr std::intmax_t first_base = get_or_compute_first_factor(); static constexpr std::intmax_t first_power = multiplicity(first_base, N); static constexpr std::intmax_t remainder = remove_power(first_base, first_power, N); diff --git a/test/unit_test/runtime/magnitude_test.cpp b/test/unit_test/runtime/magnitude_test.cpp index be32bb73..9db7b696 100644 --- a/test/unit_test/runtime/magnitude_test.cpp +++ b/test/unit_test/runtime/magnitude_test.cpp @@ -25,7 +25,13 @@ #include #include -namespace units { +using namespace units; +using namespace units::detail; + +template<> +inline constexpr std::optional units::known_first_factor<9223372036854775783> = 9223372036854775783; + +namespace { // A set of non-standard bases for testing purposes. struct noninteger_base { @@ -149,6 +155,17 @@ TEST_CASE("make_ratio performs prime factorization correctly") // The failure was due to a prime factor which is larger than 2^31. as_magnitude(); } + + 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. + as_magnitude<9'223'372'036'854'775'783>(); + } } TEST_CASE("magnitude converts to numerical value") @@ -334,7 +351,8 @@ TEST_CASE("can distinguish integral, rational, and irrational magnitudes") } } -namespace detail { +//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// +// Detail function tests below. TEST_CASE("int_power computes integer powers") { @@ -358,16 +376,6 @@ TEST_CASE("int_power computes integer powers") TEST_CASE("Prime helper functions") { - SECTION("find_first_factor()") - { - CHECK(find_first_factor(1) == 1); - CHECK(find_first_factor(2) == 2); - CHECK(find_first_factor(4) == 2); - CHECK(find_first_factor(6) == 2); - CHECK(find_first_factor(15) == 3); - CHECK(find_first_factor(17) == 17); - } - SECTION("multiplicity") { CHECK(multiplicity(2, 8) == 3); @@ -514,6 +522,4 @@ TEST_CASE("strictly_increasing") } } -} // namespace detail - -} // namespace units +} // namespace diff --git a/test/unit_test/static/CMakeLists.txt b/test/unit_test/static/CMakeLists.txt index 14c9604b..2d974c9f 100644 --- a/test/unit_test/static/CMakeLists.txt +++ b/test/unit_test/static/CMakeLists.txt @@ -52,6 +52,7 @@ add_library(unit_tests_static kind_test.cpp math_test.cpp point_origin_test.cpp + prime_test.cpp ratio_test.cpp references_test.cpp si_test.cpp diff --git a/test/unit_test/static/prime_test.cpp b/test/unit_test/static/prime_test.cpp new file mode 100644 index 00000000..020a1fe7 --- /dev/null +++ b/test/unit_test/static/prime_test.cpp @@ -0,0 +1,76 @@ +// 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. + +#include +#include +#include + +using namespace units::detail; + +namespace { + +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(109561) == is_prime_by_trial_division(109561)); + +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)); + +} // namespace