diff --git a/src/core/include/units/bits/prime.h b/src/core/include/units/bits/prime.h new file mode 100644 index 00000000..0871b17b --- /dev/null +++ b/src/core/include/units/bits/prime.h @@ -0,0 +1,141 @@ +// 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 + +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; +} + +constexpr std::size_t num_primes_between(std::size_t start, std::size_t end) { + std::size_t n = 0; + + for (auto k = start; k < end; ++k) { + if (is_prime_by_trial_division(k)) { ++n; } + } + + return n; +} + +template +constexpr auto primes_between() { + std::array results; + auto iter = std::begin(results); + + for (auto k = Start; k < End; ++k) { + if (is_prime_by_trial_division(k)) { + *iter = k; + ++iter; + } + } + + assert(iter == std::end(results)); + return results; +} + +// 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 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 prime numbers that are +// less than the _product_ of the first N prime numbers. This means it grows rapidly with N. Consider this approximate +// chart: +// +// N | Num primes stored | Trial divisions needed ~= (stored + 1 - N) / product(basis) +// --+-------------------+----------------------- +// 1 | 1 | 50.0 % +// 2 | 3 | 33.3 % +// 3 | 10 | 26.7 % +// 4 | 46 | 20.5 % +// 5 | 343 | 14.7 % +// +// Consider this behaviour when choosing the value of N most appropriate for your needs, and watch out for diminishing +// returns. +// +// [1] https://en.wikipedia.org/wiki/Wheel_factorization +template +struct WheelFactorizer { + static constexpr auto basis = first_n_primes(); + static constexpr std::size_t wheel_size = product(basis); + static constexpr auto primes_in_first_wheel = primes_between(); + + 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 (const auto& p : primes_in_first_wheel) { + if (const auto k = first_factor_maybe(n, p)) { return *k; } + } + + for (std::size_t wheel = wheel_size; wheel < n; wheel += wheel_size) { + if (const auto k = first_factor_maybe(n, wheel + 1)) { return *k; } + + for (const auto &p : primes_in_first_wheel) { + if (const auto k = first_factor_maybe(n, wheel + p)) { return *k; } + } + } + } + + 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 3860f326..cb672c88 100644 --- a/src/core/include/units/magnitude.h +++ b/src/core/include/units/magnitude.h @@ -23,11 +23,17 @@ #pragma once #include +#include #include #include #include namespace units { +namespace detail +{ +// Higher numbers use fewer trial divisions, at the price of more storage space. +using Factorizer = WheelFactorizer<3>; +} // namespace detail /** * @brief Any type which can be used as a basis vector in a BasePower. @@ -229,16 +235,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) { @@ -261,7 +257,9 @@ 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) { if (bp.power == 0) { return false; } @@ -442,7 +440,7 @@ namespace detail { template requires (N > 0) struct prime_factorization { - static constexpr std::intmax_t first_base = find_first_factor(N); + static constexpr std::intmax_t first_base = Factorizer::find_first_factor(N); 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 3def404b..667adfc3 100644 --- a/test/unit_test/runtime/magnitude_test.cpp +++ b/test/unit_test/runtime/magnitude_test.cpp @@ -140,6 +140,15 @@ 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 handle prime number which would exceed GCC iteration limit") + { + // GCC 10 has a constexpr loop iteration limit of 262144. A naive algorithm, which performs trial division on 2 and + // all odd numbers up to sqrt(N), will exceed this limit for the following prime. Thus, for this test to pass, we + // need to be using a more efficient algorithm. (We could increase the limit, but we don't want users to have to + // mess with compiler flags just to compile the code.) + as_magnitude<334524384739>(); + } } TEST_CASE("magnitude converts to numerical value") @@ -344,15 +353,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); CHECK(multiplicity(2, 1024) == 10); 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..650f5b64 --- /dev/null +++ b/test/unit_test/static/prime_test.cpp @@ -0,0 +1,57 @@ +// 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 + +namespace units::detail +{ + +template +constexpr bool check_primes(std::index_sequence) { + return ((Is < 2 || WheelFactorizer::is_prime(Is) == is_prime_by_trial_division(Is)) && ...); +} + +static_assert(check_primes<2>(std::make_index_sequence<122>{})); + +constexpr auto some_primes = primes_between<7, 19>(); // 7, 11, 13, 17 +static_assert(std::is_same_v, std::array>); +static_assert(some_primes[0] == 7); +static_assert(some_primes[1] == 11); +static_assert(some_primes[2] == 13); +static_assert(some_primes[3] == 17); + +static_assert(!WheelFactorizer<1>::is_prime(0)); +static_assert(!WheelFactorizer<1>::is_prime(1)); +static_assert(WheelFactorizer<1>::is_prime(2)); + +static_assert(!WheelFactorizer<2>::is_prime(0)); +static_assert(!WheelFactorizer<2>::is_prime(1)); +static_assert(WheelFactorizer<2>::is_prime(2)); + +static_assert(!WheelFactorizer<3>::is_prime(0)); +static_assert(!WheelFactorizer<3>::is_prime(1)); +static_assert(WheelFactorizer<3>::is_prime(2)); + +} // namespace units::detail