From 04b80f0827b56baa43fcd50b69ffcfe763a88c92 Mon Sep 17 00:00:00 2001 From: Chip Hogg Date: Thu, 10 Mar 2022 23:29:02 +0000 Subject: [PATCH 01/20] 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. --- src/core/include/units/bits/prime.h | 141 ++++++++++++++++++++++ src/core/include/units/magnitude.h | 22 ++-- test/unit_test/runtime/magnitude_test.cpp | 18 +-- test/unit_test/static/CMakeLists.txt | 1 + test/unit_test/static/prime_test.cpp | 57 +++++++++ 5 files changed, 218 insertions(+), 21 deletions(-) create mode 100644 src/core/include/units/bits/prime.h create mode 100644 test/unit_test/static/prime_test.cpp 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 From 24b284fcbbd0467935494441da775fc458395b35 Mon Sep 17 00:00:00 2001 From: Chip Hogg Date: Thu, 10 Mar 2022 23:41:36 +0000 Subject: [PATCH 02/20] Add tests to support claims in comment --- test/unit_test/static/prime_test.cpp | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/test/unit_test/static/prime_test.cpp b/test/unit_test/static/prime_test.cpp index 650f5b64..8360a60e 100644 --- a/test/unit_test/static/prime_test.cpp +++ b/test/unit_test/static/prime_test.cpp @@ -42,6 +42,11 @@ static_assert(some_primes[1] == 11); static_assert(some_primes[2] == 13); static_assert(some_primes[3] == 17); +static_assert(WheelFactorizer<1>::primes_in_first_wheel.size() == 0); +static_assert(WheelFactorizer<2>::primes_in_first_wheel.size() == 1); +static_assert(WheelFactorizer<3>::primes_in_first_wheel.size() == 7); +static_assert(WheelFactorizer<4>::primes_in_first_wheel.size() == 42); + static_assert(!WheelFactorizer<1>::is_prime(0)); static_assert(!WheelFactorizer<1>::is_prime(1)); static_assert(WheelFactorizer<1>::is_prime(2)); From c8a44adee2e934a98e9762c425a0650c5377cb15 Mon Sep 17 00:00:00 2001 From: Chip Hogg Date: Thu, 10 Mar 2022 23:41:47 +0000 Subject: [PATCH 03/20] Add missing header for `std::integral` --- src/core/include/units/magnitude.h | 1 + 1 file changed, 1 insertion(+) diff --git a/src/core/include/units/magnitude.h b/src/core/include/units/magnitude.h index cb672c88..53630029 100644 --- a/src/core/include/units/magnitude.h +++ b/src/core/include/units/magnitude.h @@ -24,6 +24,7 @@ #include #include +#include #include #include #include From 6c73947fe00cfddd86392f444e8d8c7df46bfc47 Mon Sep 17 00:00:00 2001 From: Chip Hogg Date: Thu, 10 Mar 2022 23:50:33 +0000 Subject: [PATCH 04/20] Satisfy complaint This line shouldn't _actually_ ever be reachable, but I can't fault the compiler for not figuring that out. --- src/core/include/units/bits/prime.h | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/core/include/units/bits/prime.h b/src/core/include/units/bits/prime.h index 0871b17b..4c830b75 100644 --- a/src/core/include/units/bits/prime.h +++ b/src/core/include/units/bits/prime.h @@ -131,6 +131,8 @@ struct WheelFactorizer { if (const auto k = first_factor_maybe(n, wheel + p)) { return *k; } } } + + return n; } static constexpr bool is_prime(std::size_t n) { From 73a56115a197d605ccd979c566232369091ae3b1 Mon Sep 17 00:00:00 2001 From: Chip Hogg Date: Fri, 11 Mar 2022 03:42:19 +0000 Subject: [PATCH 05/20] Fix wheel factorization algorithm --- src/core/include/units/bits/prime.h | 88 +++++++++++++++------------- test/unit_test/static/prime_test.cpp | 33 +++++++---- 2 files changed, 69 insertions(+), 52 deletions(-) diff --git a/src/core/include/units/bits/prime.h b/src/core/include/units/bits/prime.h index 4c830b75..4baed1c4 100644 --- a/src/core/include/units/bits/prime.h +++ b/src/core/include/units/bits/prime.h @@ -26,6 +26,7 @@ #include #include #include +#include namespace units::detail { @@ -37,32 +38,6 @@ constexpr bool is_prime_by_trial_division(std::size_t n) { 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. @@ -83,6 +58,34 @@ constexpr std::array first_n_primes() { 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; }); + + assert (i == N); + + return coprimes; +} + template constexpr std::size_t product(const std::array &values) { std::size_t product = 1; @@ -93,41 +96,42 @@ constexpr std::size_t product(const std::array &values) { // 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: +// 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 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 % +// 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 % // -// Consider this behaviour when choosing the value of N most appropriate for your needs, and watch out for diminishing -// returns. +// 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 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 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 (const auto& p : primes_in_first_wheel) { - 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) { if (const auto k = first_factor_maybe(n, wheel + 1)) { return *k; } - for (const auto &p : primes_in_first_wheel) { + for (const auto &p : coprimes_in_first_wheel) { if (const auto k = first_factor_maybe(n, wheel + p)) { return *k; } } } diff --git a/test/unit_test/static/prime_test.cpp b/test/unit_test/static/prime_test.cpp index 8360a60e..ab852974 100644 --- a/test/unit_test/static/prime_test.cpp +++ b/test/unit_test/static/prime_test.cpp @@ -35,17 +35,30 @@ constexpr bool check_primes(std::index_sequence) { 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); +// 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(WheelFactorizer<4>::is_prime(109561) == is_prime_by_trial_division(109561)); -static_assert(WheelFactorizer<1>::primes_in_first_wheel.size() == 0); -static_assert(WheelFactorizer<2>::primes_in_first_wheel.size() == 1); -static_assert(WheelFactorizer<3>::primes_in_first_wheel.size() == 7); -static_assert(WheelFactorizer<4>::primes_in_first_wheel.size() == 42); +static_assert(WheelFactorizer<1>::coprimes_in_first_wheel.size() == 1); +static_assert(WheelFactorizer<2>::coprimes_in_first_wheel.size() == 2); +static_assert(WheelFactorizer<3>::coprimes_in_first_wheel.size() == 8); +static_assert(WheelFactorizer<4>::coprimes_in_first_wheel.size() == 48); +static_assert(WheelFactorizer<5>::coprimes_in_first_wheel.size() == 480); + +static_assert(WheelFactorizer<3>::coprimes_in_first_wheel[0] == 1); +static_assert(WheelFactorizer<3>::coprimes_in_first_wheel[1] == 7); +static_assert(WheelFactorizer<3>::coprimes_in_first_wheel[2] == 11); +static_assert(WheelFactorizer<3>::coprimes_in_first_wheel[3] == 13); +static_assert(WheelFactorizer<3>::coprimes_in_first_wheel[4] == 17); +static_assert(WheelFactorizer<3>::coprimes_in_first_wheel[5] == 19); +static_assert(WheelFactorizer<3>::coprimes_in_first_wheel[6] == 23); +static_assert(WheelFactorizer<3>::coprimes_in_first_wheel[7] == 29); static_assert(!WheelFactorizer<1>::is_prime(0)); static_assert(!WheelFactorizer<1>::is_prime(1)); From bfa8db613903fd50ce45c037b0d8ad505e85682f Mon Sep 17 00:00:00 2001 From: Chip Hogg Date: Fri, 11 Mar 2022 03:48:57 +0000 Subject: [PATCH 06/20] Use std::accumulate --- src/core/include/units/bits/prime.h | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/src/core/include/units/bits/prime.h b/src/core/include/units/bits/prime.h index 4baed1c4..47eb2dbe 100644 --- a/src/core/include/units/bits/prime.h +++ b/src/core/include/units/bits/prime.h @@ -25,6 +25,7 @@ #include #include #include +#include #include #include @@ -88,9 +89,7 @@ constexpr auto coprimes_up_to(std::size_t n, const std::array &b template constexpr std::size_t product(const std::array &values) { - std::size_t product = 1; - for (const auto &v : values) { product *= v; } - return product; + return std::accumulate(std::begin(values), std::end(values), std::size_t{1u}, std::multiplies{}); } // A configurable instantiation of the "wheel factorization" algorithm [1] for prime numbers. From 59d9cd14079a0acaf0283635e1057eabbe329901 Mon Sep 17 00:00:00 2001 From: Chip Hogg Date: Fri, 11 Mar 2022 03:50:08 +0000 Subject: [PATCH 07/20] static_cast for first factor --- src/core/include/units/magnitude.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/core/include/units/magnitude.h b/src/core/include/units/magnitude.h index 53630029..dc3ac304 100644 --- a/src/core/include/units/magnitude.h +++ b/src/core/include/units/magnitude.h @@ -441,7 +441,7 @@ namespace detail { template requires (N > 0) struct prime_factorization { - static constexpr std::intmax_t first_base = Factorizer::find_first_factor(N); + static constexpr std::intmax_t first_base = static_cast(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); From 70640a1017c11e9bc8af831d632a745a177f9d73 Mon Sep 17 00:00:00 2001 From: Chip Hogg Date: Fri, 11 Mar 2022 03:56:57 +0000 Subject: [PATCH 08/20] Remove constexpr-incompatible assert() --- src/core/include/units/bits/prime.h | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/core/include/units/bits/prime.h b/src/core/include/units/bits/prime.h index 47eb2dbe..5633dc4e 100644 --- a/src/core/include/units/bits/prime.h +++ b/src/core/include/units/bits/prime.h @@ -82,8 +82,6 @@ constexpr auto coprimes_up_to(std::size_t n, const std::array &b call_for_coprimes_up_to(n, basis, [&coprimes, &i](std::size_t cp) { coprimes[i++] = cp; }); - assert (i == N); - return coprimes; } From a719a8b91233fc03ed5e0940d5bb36d15cc524f1 Mon Sep 17 00:00:00 2001 From: Chip Hogg Date: Fri, 11 Mar 2022 03:57:21 +0000 Subject: [PATCH 09/20] Try upping the basis size --- src/core/include/units/magnitude.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/core/include/units/magnitude.h b/src/core/include/units/magnitude.h index dc3ac304..5f01a0fa 100644 --- a/src/core/include/units/magnitude.h +++ b/src/core/include/units/magnitude.h @@ -33,7 +33,7 @@ namespace units { namespace detail { // Higher numbers use fewer trial divisions, at the price of more storage space. -using Factorizer = WheelFactorizer<3>; +using Factorizer = WheelFactorizer<4>; } // namespace detail /** From 87073856a7f18b365c0670d87b2380942705dd9d Mon Sep 17 00:00:00 2001 From: Chip Hogg Date: Sat, 12 Mar 2022 18:50:01 +0000 Subject: [PATCH 10/20] Try upping the basis size further We are well into a regime of diminishing returns, but we'd better start by seeing if the easy thing works. Besides, setting this to 7 trips the step limit in _generating_ the algorithm! --- src/core/include/units/magnitude.h | 2 +- test/unit_test/static/prime_test.cpp | 1 + 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/src/core/include/units/magnitude.h b/src/core/include/units/magnitude.h index 5f01a0fa..fbe93316 100644 --- a/src/core/include/units/magnitude.h +++ b/src/core/include/units/magnitude.h @@ -33,7 +33,7 @@ namespace units { namespace detail { // Higher numbers use fewer trial divisions, at the price of more storage space. -using Factorizer = WheelFactorizer<4>; +using Factorizer = WheelFactorizer<6>; } // namespace detail /** diff --git a/test/unit_test/static/prime_test.cpp b/test/unit_test/static/prime_test.cpp index ab852974..cfb57273 100644 --- a/test/unit_test/static/prime_test.cpp +++ b/test/unit_test/static/prime_test.cpp @@ -50,6 +50,7 @@ static_assert(WheelFactorizer<2>::coprimes_in_first_wheel.size() == 2); static_assert(WheelFactorizer<3>::coprimes_in_first_wheel.size() == 8); static_assert(WheelFactorizer<4>::coprimes_in_first_wheel.size() == 48); static_assert(WheelFactorizer<5>::coprimes_in_first_wheel.size() == 480); +static_assert(WheelFactorizer<6>::coprimes_in_first_wheel.size() == 5760); static_assert(WheelFactorizer<3>::coprimes_in_first_wheel[0] == 1); static_assert(WheelFactorizer<3>::coprimes_in_first_wheel[1] == 7); From 1e8460d401ada0004eb8599986ebca83e01d05e4 Mon Sep 17 00:00:00 2001 From: Chip Hogg Date: Sat, 19 Mar 2022 13:21:47 +0000 Subject: [PATCH 11/20] Revert "Try upping the basis size further" This reverts commit 87073856a7f18b365c0670d87b2380942705dd9d. It didn't fix the problem, and it caused some new ones. We need a different approach. --- src/core/include/units/magnitude.h | 2 +- test/unit_test/static/prime_test.cpp | 1 - 2 files changed, 1 insertion(+), 2 deletions(-) diff --git a/src/core/include/units/magnitude.h b/src/core/include/units/magnitude.h index fbe93316..5f01a0fa 100644 --- a/src/core/include/units/magnitude.h +++ b/src/core/include/units/magnitude.h @@ -33,7 +33,7 @@ namespace units { namespace detail { // Higher numbers use fewer trial divisions, at the price of more storage space. -using Factorizer = WheelFactorizer<6>; +using Factorizer = WheelFactorizer<4>; } // namespace detail /** diff --git a/test/unit_test/static/prime_test.cpp b/test/unit_test/static/prime_test.cpp index cfb57273..ab852974 100644 --- a/test/unit_test/static/prime_test.cpp +++ b/test/unit_test/static/prime_test.cpp @@ -50,7 +50,6 @@ static_assert(WheelFactorizer<2>::coprimes_in_first_wheel.size() == 2); static_assert(WheelFactorizer<3>::coprimes_in_first_wheel.size() == 8); static_assert(WheelFactorizer<4>::coprimes_in_first_wheel.size() == 48); static_assert(WheelFactorizer<5>::coprimes_in_first_wheel.size() == 480); -static_assert(WheelFactorizer<6>::coprimes_in_first_wheel.size() == 5760); static_assert(WheelFactorizer<3>::coprimes_in_first_wheel[0] == 1); static_assert(WheelFactorizer<3>::coprimes_in_first_wheel[1] == 7); From 28c4fe3c08ad92ecd1f39c0d32c687603774517a Mon Sep 17 00:00:00 2001 From: Chip Hogg Date: Sat, 19 Mar 2022 14:29:11 +0000 Subject: [PATCH 12/20] Run clang-format-15 on changed files --- src/core/include/units/bits/prime.h | 67 +++++--- src/core/include/units/magnitude.h | 176 +++++++++++-------- test/unit_test/runtime/magnitude_test.cpp | 201 +++++++++------------- test/unit_test/static/prime_test.cpp | 9 +- 4 files changed, 238 insertions(+), 215 deletions(-) diff --git a/src/core/include/units/bits/prime.h b/src/core/include/units/bits/prime.h index 5633dc4e..aaa20e4f 100644 --- a/src/core/include/units/bits/prime.h +++ b/src/core/include/units/bits/prime.h @@ -29,12 +29,14 @@ #include #include -namespace units::detail -{ +namespace units::detail { -constexpr bool is_prime_by_trial_division(std::size_t n) { +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; } + if (n % f == 0) { + return false; + } } return true; } @@ -42,25 +44,34 @@ constexpr bool is_prime_by_trial_division(std::size_t n) { // 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; } +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() { +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]; } + 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) { +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); @@ -69,14 +80,16 @@ constexpr void call_for_coprimes_up_to(std::size_t n, const std::array -constexpr std::size_t num_coprimes_up_to(std::size_t n, const std::array &basis) { +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) { +constexpr auto coprimes_up_to(std::size_t n, const std::array& basis) +{ std::array coprimes; std::size_t i = 0u; @@ -86,7 +99,8 @@ constexpr auto coprimes_up_to(std::size_t n, const std::array &b } template -constexpr std::size_t product(const std::array &values) { +constexpr std::size_t product(const std::array& values) +{ return std::accumulate(std::begin(values), std::end(values), std::size_t{1u}, std::multiplies{}); } @@ -116,29 +130,36 @@ struct WheelFactorizer { 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) { + 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; } + 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; } + if (const auto k = first_factor_maybe(n, *it)) { + 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; } + if (const auto k = first_factor_maybe(n, wheel + 1)) { + return *k; + } - for (const auto &p : coprimes_in_first_wheel) { - if (const auto k = first_factor_maybe(n, wheel + p)) { return *k; } + 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; - } + static constexpr bool is_prime(std::size_t n) { return (n > 1) && find_first_factor(n) == n; } }; -} // namespace units::detail +} // namespace units::detail diff --git a/src/core/include/units/magnitude.h b/src/core/include/units/magnitude.h index 5f01a0fa..3cec3b98 100644 --- a/src/core/include/units/magnitude.h +++ b/src/core/include/units/magnitude.h @@ -22,19 +22,19 @@ #pragma once -#include +#include #include +#include #include #include #include #include namespace units { -namespace detail -{ +namespace detail { // Higher numbers use fewer trial divisions, at the price of more storage space. using Factorizer = WheelFactorizer<4>; -} // namespace detail +} // namespace detail /** * @brief Any type which can be used as a basis vector in a BasePower. @@ -52,7 +52,8 @@ using Factorizer = WheelFactorizer<4>; * GCC 10) which don't yet permit floating point NTTPs. */ template -concept BaseRep = std::is_same_v || std::is_same_v, long double>; +concept BaseRep = std::is_same_v || std::is_same_v < std::remove_cvref_t, +long double > ; /** * @brief A basis vector in our magnitude representation, raised to some rational power. @@ -114,7 +115,7 @@ template inline constexpr bool is_base_power = false; template inline constexpr bool is_base_power> = true; -} // namespace detail +} // namespace detail /** * @brief Concept to detect whether a _type_ is a valid base power. @@ -125,34 +126,40 @@ inline constexpr bool is_base_power> = true; template concept BasePower = detail::is_base_power; -namespace detail +namespace detail { + +constexpr auto inverse(BasePower auto bp) { -constexpr auto inverse(BasePower auto bp) { bp.power.num *= -1; return bp; } // `widen_t` gives the widest arithmetic type in the same category, for intermediate computations. -template requires std::is_arithmetic_v -using widen_t = std::conditional_t< - std::is_floating_point_v, - long double, - std::conditional_t, std::intmax_t, std::uintmax_t>>; +template + requires std::is_arithmetic_v +using widen_t = std::conditional_t, long double, + std::conditional_t, std::intmax_t, std::uintmax_t>>; // Raise an arbitrary arithmetic type to a positive integer power at compile time. -template requires std::is_arithmetic_v -constexpr T int_power(T base, std::integral auto exp){ +template + requires std::is_arithmetic_v +constexpr T int_power(T base, std::integral auto exp) +{ // As this function should only be called at compile time, the exceptions herein function as // "parameter-compatible static_asserts", and should not result in exceptions at runtime. - if (exp < 0) { throw std::invalid_argument{"int_power only supports positive integer powers"}; } + if (exp < 0) { + throw std::invalid_argument{"int_power only supports positive integer powers"}; + } - constexpr auto checked_multiply = [] (auto a, auto b) { + constexpr auto checked_multiply = [](auto a, auto b) { const auto result = a * b; - if (result / a != b) { throw std::overflow_error{"Wraparound detected"}; } + if (result / a != b) { + throw std::overflow_error{"Wraparound detected"}; + } return result; }; - constexpr auto checked_square = [checked_multiply] (auto a) { return checked_multiply(a, a); }; + constexpr auto checked_square = [checked_multiply](auto a) { return checked_multiply(a, a); }; // TODO(chogg): Unify this implementation with the one in pow.h. That one takes its exponent as a // template parameter, rather than a function parameter. @@ -169,7 +176,8 @@ constexpr T int_power(T base, std::integral auto exp){ } -template requires std::is_arithmetic_v +template + requires std::is_arithmetic_v constexpr widen_t compute_base_power(BasePower auto bp) { // This utility can only handle integer powers. To compute rational powers at compile time, we'll @@ -177,8 +185,12 @@ constexpr widen_t compute_base_power(BasePower auto bp) // // Note that since this function should only be called at compile time, the point of these // exceptions is to act as "static_assert substitutes", not to throw actual exceptions at runtime. - if (bp.power.den != 1) { throw std::invalid_argument{"Rational powers not yet supported"}; } - if (bp.power.exp < 0) { throw std::invalid_argument{"Unsupported exp value"}; } + if (bp.power.den != 1) { + throw std::invalid_argument{"Rational powers not yet supported"}; + } + if (bp.power.exp < 0) { + throw std::invalid_argument{"Unsupported exp value"}; + } if (bp.power.num < 0) { if constexpr (std::is_integral_v) { @@ -197,10 +209,10 @@ constexpr widen_t compute_base_power(BasePower auto bp) // The input is the desired result, but in a (wider) intermediate type. The point of this function // is to cast to the desired type, but avoid overflow in doing so. template - requires std::is_arithmetic_v - && std::is_arithmetic_v - && (std::is_integral_v == std::is_integral_v) -constexpr To checked_static_cast(From x) { + requires std::is_arithmetic_v && std::is_arithmetic_v && + (std::is_integral_v == std::is_integral_v) +constexpr To checked_static_cast(From x) +{ // This function should only ever be called at compile time. The purpose of these exceptions is // to produce compiler errors, because we cannot `static_assert` on function arguments. if constexpr (std::is_integral_v) { @@ -215,20 +227,22 @@ constexpr To checked_static_cast(From x) { return static_cast(x); } -} // namespace detail +} // namespace detail /** * @brief Equality detection for two base powers. */ template -constexpr bool operator==(T t, U u) { +constexpr bool operator==(T t, U u) +{ return std::is_same_v && (t.get_base() == u.get_base()) && (t.power == u.power); } /** * @brief A BasePower, raised to a rational power E. */ -constexpr auto pow(BasePower auto bp, ratio p) { +constexpr auto pow(BasePower auto bp, ratio p) +{ bp.power = bp.power * p; return bp; } @@ -240,8 +254,7 @@ namespace detail { constexpr std::intmax_t multiplicity(std::intmax_t factor, std::intmax_t n) { std::intmax_t m = 0; - while (n % factor == 0) - { + while (n % factor == 0) { n /= factor; ++m; } @@ -253,20 +266,26 @@ constexpr std::intmax_t multiplicity(std::intmax_t factor, std::intmax_t n) // Undefined unless base > 1, pow >= 0, and (base ^ pow) evenly divides n. constexpr std::intmax_t remove_power(std::intmax_t base, std::intmax_t pow, std::intmax_t n) { - while (pow-- > 0) { n /= base; } + while (pow-- > 0) { + n /= base; + } return n; } // A way to check whether a number is prime at compile time. -constexpr bool is_prime(std::intmax_t n) { - return (n >= 0) && Factorizer::is_prime(static_cast(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; } +constexpr bool is_valid_base_power(const BasePower auto& bp) +{ + if (bp.power == 0) { + return false; + } - if constexpr (std::is_same_v) { return is_prime(bp.get_base()); } - else { return bp.get_base() > 0; } + if constexpr (std::is_same_v) { + return is_prime(bp.get_base()); + } else { + return bp.get_base() > 0; + } } // A function object to apply a predicate to all consecutive pairs of values in a sequence. @@ -275,16 +294,19 @@ struct pairwise_all { Predicate predicate; template - constexpr bool operator()(Ts&&... ts) const { + constexpr bool operator()(Ts&&... ts) const + { // Carefully handle different sizes, avoiding unsigned integer underflow. constexpr auto num_comparisons = [](auto num_elements) { return (num_elements > 1) ? (num_elements - 1) : 0; }(sizeof...(Ts)); // Compare zero or more pairs of neighbours as needed. - return [this](std::tuple &&t, std::index_sequence) { + return [this](std::tuple && t, std::index_sequence) + { return (predicate(std::get(t), std::get(t)) && ...); - }(std::make_tuple(std::forward(ts)...), std::make_index_sequence()); + } + (std::make_tuple(std::forward(ts)...), std::make_index_sequence()); } }; @@ -294,8 +316,9 @@ pairwise_all(T) -> pairwise_all; // Check whether a sequence of (possibly heterogeneously typed) values are strictly increasing. template - requires (std::is_signed_v && ...) -constexpr bool strictly_increasing(Ts&&... ts) { + requires(std::is_signed_v && ...) +constexpr bool strictly_increasing(Ts&&... ts) +{ return pairwise_all{std::less{}}(std::forward(ts)...); } @@ -308,14 +331,13 @@ inline constexpr bool all_bases_in_order = strictly_increasing(BPs.get_base()... template inline constexpr bool is_base_power_pack_valid = all_base_powers_valid && all_bases_in_order; -constexpr bool is_rational(BasePower auto bp) { +constexpr bool is_rational(BasePower auto bp) +{ return std::is_integral_v && (bp.power.den == 1) && (bp.power.exp >= 0); } -constexpr bool is_integral(BasePower auto bp) { - return is_rational(bp) && bp.power.num > 0; -} -} // namespace detail +constexpr bool is_integral(BasePower auto bp) { return is_rational(bp) && bp.power.num > 0; } +} // namespace detail /** * @brief A representation for positive real numbers which optimizes taking products and rational powers. @@ -339,7 +361,7 @@ template inline constexpr bool is_magnitude = false; template inline constexpr bool is_magnitude> = true; -} // namespace detail +} // namespace detail /** * @brief Concept to detect whether T is a valid Magnitude. @@ -352,7 +374,8 @@ concept Magnitude = detail::is_magnitude; */ template requires std::is_floating_point_v || (std::is_integral_v && is_integral(magnitude{})) -constexpr T get_value(const magnitude &) { +constexpr T get_value(const magnitude&) +{ // Force the expression to be evaluated in a constexpr context, to catch, e.g., overflow. constexpr auto result = detail::checked_static_cast((detail::compute_base_power(BPs) * ...)); @@ -370,18 +393,26 @@ struct pi_base { // Magnitude equality implementation. template -constexpr bool operator==(magnitude, magnitude) { - if constexpr (sizeof...(LeftBPs) == sizeof...(RightBPs)) { return ((LeftBPs == RightBPs) && ...); } - else { return false; } +constexpr bool operator==(magnitude, magnitude) +{ + if constexpr (sizeof...(LeftBPs) == sizeof...(RightBPs)) { + return ((LeftBPs == RightBPs) && ...); + } else { + return false; + } } //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// // Magnitude rational powers implementation. template -constexpr auto pow(magnitude) { - if constexpr (E == 0) { return magnitude<>{}; } - else { return magnitude{}; } +constexpr auto pow(magnitude) +{ + if constexpr (E == 0) { + return magnitude<>{}; + } else { + return magnitude{}; + } } //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// @@ -394,9 +425,10 @@ constexpr auto operator*(Magnitude auto m, magnitude<>) { return m; } // Recursive case for the product of any two non-identity Magnitudes. template -constexpr auto operator*(magnitude, magnitude) { +constexpr auto operator*(magnitude, magnitude) +{ // Case for when H1 has the smaller base. - if constexpr(H1.get_base() < H2.get_base()){ + if constexpr (H1.get_base() < H2.get_base()) { if constexpr (sizeof...(T1) == 0) { // Shortcut for the "pure prepend" case, which makes it easier to implement some of the other cases. return magnitude{}; @@ -406,7 +438,7 @@ constexpr auto operator*(magnitude, magnitude) { } // Case for when H2 has the smaller base. - if constexpr(H1.get_base() > H2.get_base()){ + if constexpr (H1.get_base() > H2.get_base()) { return magnitude

{} * (magnitude{} * magnitude{}); } @@ -439,23 +471,25 @@ constexpr auto operator/(Magnitude auto l, Magnitude auto r) { return l * pow<-1 namespace detail { // Helper to perform prime factorization at compile time. template - requires (N > 0) + requires(N > 0) struct prime_factorization { static constexpr std::intmax_t first_base = static_cast(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); - static constexpr auto value = magnitude{} - * prime_factorization::value; + static constexpr auto value = + magnitude{} * prime_factorization::value; }; // Specialization for the prime factorization of 1 (base case). template<> -struct prime_factorization<1> { static constexpr magnitude<> value{}; }; +struct prime_factorization<1> { + static constexpr magnitude<> value{}; +}; template inline constexpr auto prime_factorization_v = prime_factorization::value; -} // namespace detail +} // namespace detail /** * @brief Convert any positive integer to a Magnitude. @@ -464,11 +498,11 @@ inline constexpr auto prime_factorization_v = prime_factorization::value; * manually adding base powers. */ template - requires (R.num > 0) -constexpr Magnitude auto as_magnitude() { - return pow(detail::prime_factorization_v<10>) - * detail::prime_factorization_v - / detail::prime_factorization_v; + requires(R.num > 0) +constexpr Magnitude auto as_magnitude() +{ + return pow(detail::prime_factorization_v<10>) * detail::prime_factorization_v / + detail::prime_factorization_v; } -} // namespace units +} // namespace units diff --git a/test/unit_test/runtime/magnitude_test.cpp b/test/unit_test/runtime/magnitude_test.cpp index 667adfc3..e9b359f2 100644 --- a/test/unit_test/runtime/magnitude_test.cpp +++ b/test/unit_test/runtime/magnitude_test.cpp @@ -20,54 +20,64 @@ // OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE // SOFTWARE. +#include #include #include -#include #include namespace units { // A set of non-standard bases for testing purposes. -struct noninteger_base { static constexpr long double value = 1.234L; }; -struct noncanonical_two_base { static constexpr long double value = 2.0L; }; -struct other_noncanonical_two_base { static constexpr long double value = 2.0L; }; -struct invalid_zero_base { static constexpr long double value = 0.0L; }; -struct invalid_negative_base { static constexpr long double value = -1.234L; }; +struct noninteger_base { + static constexpr long double value = 1.234L; +}; +struct noncanonical_two_base { + static constexpr long double value = 2.0L; +}; +struct other_noncanonical_two_base { + static constexpr long double value = 2.0L; +}; +struct invalid_zero_base { + static constexpr long double value = 0.0L; +}; +struct invalid_negative_base { + static constexpr long double value = -1.234L; +}; template -constexpr auto pi_to_the() { return magnitude{Power}>{}; } +constexpr auto pi_to_the() +{ + return magnitude{Power}>{}; +} template -void check_same_type_and_value(T actual, U expected) { +void check_same_type_and_value(T actual, U expected) +{ CHECK(std::is_same_v); CHECK(actual == expected); } TEST_CASE("base_power") { - SECTION("base rep deducible for integral base") - { + SECTION ("base rep deducible for integral base") { CHECK(base_power{2} == base_power{2, ratio{1}}); CHECK(base_power{2, 3} == base_power{2, ratio{3}}); CHECK(base_power{2, ratio{3, 4}} == base_power{2, ratio{3, 4}}); } - SECTION("get_base retrieves base for integral base") - { + SECTION ("get_base retrieves base for integral base") { CHECK(base_power{2}.get_base() == 2); CHECK(base_power{3, 5}.get_base() == 3); CHECK(base_power{5, ratio{1, 3}}.get_base() == 5); } - SECTION("get_base retrieves member value for non-integer base") - { + SECTION ("get_base retrieves member value for non-integer base") { CHECK(base_power{}.get_base() == 1.234L); CHECK(base_power{2}.get_base() == 1.234L); CHECK(base_power{ratio{5, 8}}.get_base() == 1.234L); } - SECTION("same-base values not equal if types are different") - { + SECTION ("same-base values not equal if types are different") { const auto a = base_power{}; const auto b = base_power{2}; const auto c = base_power{}; @@ -79,31 +89,25 @@ TEST_CASE("base_power") CHECK(a != c); } - SECTION("same-type values not equal if bases are different") - { + SECTION ("same-type values not equal if bases are different") { CHECK(base_power{2} != base_power{3}); CHECK(base_power{2, ratio{5, 4}} != base_power{3, ratio{5, 4}}); } - SECTION("same-type, same-base values not equal if powers are different") - { + SECTION ("same-type, same-base values not equal if powers are different") { CHECK(base_power{2} != base_power{2, 2}); CHECK(base_power{} != base_power{ratio{1, 3}}); } - SECTION("product with inverse equals identity") - { - auto check_product_with_inverse_is_identity = [] (auto x) { - CHECK(x * pow<-1>(x) == as_magnitude<1>()); - }; + SECTION ("product with inverse equals identity") { + auto check_product_with_inverse_is_identity = [](auto x) { CHECK(x * pow<-1>(x) == as_magnitude<1>()); }; check_product_with_inverse_is_identity(as_magnitude<3>()); check_product_with_inverse_is_identity(as_magnitude()); check_product_with_inverse_is_identity(pi_to_the()); } - SECTION("pow() multiplies exponent") - { + SECTION ("pow() multiplies exponent") { CHECK(pow(base_power{2}, 0) == base_power{2, 0}); CHECK(pow(base_power{2, 3}, ratio{-1, 2}) == base_power{2, ratio{-3, 2}}); CHECK(pow(base_power{ratio{3, 2}}, ratio{1, 3}) == base_power{ratio{1, 2}}); @@ -112,8 +116,7 @@ TEST_CASE("base_power") TEST_CASE("make_ratio performs prime factorization correctly") { - SECTION("Performs prime factorization when denominator is 1") - { + SECTION ("Performs prime factorization when denominator is 1") { CHECK(as_magnitude<1>() == magnitude<>{}); CHECK(as_magnitude<2>() == magnitude{}); CHECK(as_magnitude<3>() == magnitude{}); @@ -122,27 +125,23 @@ TEST_CASE("make_ratio performs prime factorization correctly") CHECK(as_magnitude<792>() == magnitude{}); } - SECTION("Supports fractions") - { + SECTION ("Supports fractions") { CHECK(as_magnitude() == magnitude{}); } - SECTION("Supports nonzero exp") - { + SECTION ("Supports nonzero exp") { constexpr ratio r{3, 1, 2}; REQUIRE(r.exp == 2); CHECK(as_magnitude() == as_magnitude<300>()); } - SECTION("Can handle prime factor which would be large enough to overflow int") - { + SECTION ("Can handle prime factor which would be large enough to overflow int") { // This was taken from a case which failed when we used `int` for our base to store prime numbers. // 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") - { + 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 @@ -153,8 +152,7 @@ TEST_CASE("make_ratio performs prime factorization correctly") TEST_CASE("magnitude converts to numerical value") { - SECTION("Positive integer powers of integer bases give integer values") - { + SECTION ("Positive integer powers of integer bases give integer values") { constexpr auto mag_412 = as_magnitude<412>(); check_same_type_and_value(get_value(mag_412), 412); check_same_type_and_value(get_value(mag_412), std::size_t{412}); @@ -162,23 +160,20 @@ TEST_CASE("magnitude converts to numerical value") check_same_type_and_value(get_value(mag_412), 412.0); } - SECTION("Negative integer powers of integer bases compute correct values") - { + SECTION ("Negative integer powers of integer bases compute correct values") { constexpr auto mag_0p125 = as_magnitude(); check_same_type_and_value(get_value(mag_0p125), 0.125f); check_same_type_and_value(get_value(mag_0p125), 0.125); } - SECTION("pi to the 1 supplies correct values") - { + SECTION ("pi to the 1 supplies correct values") { constexpr auto pi = pi_to_the<1>(); check_same_type_and_value(get_value(pi), std::numbers::pi_v); check_same_type_and_value(get_value(pi), std::numbers::pi_v); check_same_type_and_value(get_value(pi), std::numbers::pi_v); } - SECTION("pi to arbitrary power performs computations in most accurate type at compile time") - { + SECTION ("pi to arbitrary power performs computations in most accurate type at compile time") { if constexpr (sizeof(float) < sizeof(long double)) { constexpr auto pi_cubed = pi_to_the<3>(); @@ -192,8 +187,7 @@ TEST_CASE("magnitude converts to numerical value") } } - SECTION("Impossible requests are prevented at compile time") - { + SECTION ("Impossible requests are prevented at compile time") { // Naturally, we cannot actually write a test to verify a compiler error. But any of these can // be uncommented if desired to verify that it breaks the build. @@ -205,7 +199,7 @@ TEST_CASE("magnitude converts to numerical value") // Would work for pow<63>: // get_value(pow<64>(as_magnitude<2>())); - get_value(pow<308>(as_magnitude<10>())); // Compiles, correctly. + get_value(pow<308>(as_magnitude<10>())); // Compiles, correctly. // get_value(pow<309>(as_magnitude<10>())); // get_value(pow<3099>(as_magnitude<10>())); // get_value(pow<3099999>(as_magnitude<10>())); @@ -218,21 +212,18 @@ TEST_CASE("magnitude converts to numerical value") TEST_CASE("Equality works for magnitudes") { - SECTION("Equivalent ratios are equal") - { + SECTION ("Equivalent ratios are equal") { CHECK(as_magnitude<1>() == as_magnitude<1>()); CHECK(as_magnitude<3>() == as_magnitude<3>()); CHECK(as_magnitude() == as_magnitude()); } - SECTION("Different ratios are unequal") - { + SECTION ("Different ratios are unequal") { CHECK(as_magnitude<3>() != as_magnitude<5>()); CHECK(as_magnitude<3>() != as_magnitude()); } - SECTION("Supports constexpr") - { + SECTION ("Supports constexpr") { constexpr auto eq = (as_magnitude() == as_magnitude()); CHECK(!eq); } @@ -240,25 +231,20 @@ TEST_CASE("Equality works for magnitudes") TEST_CASE("Multiplication works for magnitudes") { - SECTION("Reciprocals reduce to null magnitude") - { + SECTION ("Reciprocals reduce to null magnitude") { CHECK(as_magnitude() * as_magnitude() == as_magnitude<1>()); } - SECTION("Products work as expected") - { + SECTION ("Products work as expected") { CHECK(as_magnitude() * as_magnitude() == as_magnitude()); } - SECTION("Products handle pi correctly") - { - CHECK( - pi_to_the<1>() * as_magnitude() * pi_to_the() == - magnitude{ratio{1, 2}}>{}); + SECTION ("Products handle pi correctly") { + CHECK(pi_to_the<1>() * as_magnitude() * pi_to_the() == + magnitude{ratio{1, 2}}>{}); } - SECTION("Supports constexpr") - { + SECTION ("Supports constexpr") { constexpr auto p = as_magnitude() * as_magnitude(); CHECK(p == as_magnitude()); } @@ -266,19 +252,16 @@ TEST_CASE("Multiplication works for magnitudes") TEST_CASE("Division works for magnitudes") { - SECTION("Dividing anything by itself reduces to null magnitude") - { + SECTION ("Dividing anything by itself reduces to null magnitude") { CHECK(as_magnitude() / as_magnitude() == as_magnitude<1>()); CHECK(as_magnitude<15>() / as_magnitude<15>() == as_magnitude<1>()); } - SECTION("Quotients work as expected") - { + SECTION ("Quotients work as expected") { CHECK(as_magnitude() / as_magnitude() == as_magnitude()); } - SECTION("Supports constexpr") - { + SECTION ("Supports constexpr") { constexpr auto q = as_magnitude() / as_magnitude(); CHECK(q == as_magnitude()); } @@ -286,29 +269,28 @@ TEST_CASE("Division works for magnitudes") TEST_CASE("Can raise Magnitudes to rational powers") { - SECTION("Anything to the 0 is 1") { + SECTION ("Anything to the 0 is 1") { CHECK(pow<0>(as_magnitude<1>()) == as_magnitude<1>()); CHECK(pow<0>(as_magnitude<123>()) == as_magnitude<1>()); CHECK(pow<0>(as_magnitude()) == as_magnitude<1>()); CHECK(pow<0>(pi_to_the()) == as_magnitude<1>()); } - SECTION("Anything to the 1 is itself") { + SECTION ("Anything to the 1 is itself") { CHECK(pow<1>(as_magnitude<1>()) == as_magnitude<1>()); CHECK(pow<1>(as_magnitude<123>()) == as_magnitude<123>()); CHECK(pow<1>(as_magnitude()) == as_magnitude()); CHECK(pow<1>(pi_to_the()) == pi_to_the()); } - SECTION("Can raise to arbitrary rational power") { + SECTION ("Can raise to arbitrary rational power") { CHECK(pow(pi_to_the()) == pi_to_the()); } } TEST_CASE("can distinguish integral, rational, and irrational magnitudes") { - SECTION("Integer magnitudes are integral and rational") - { + SECTION ("Integer magnitudes are integral and rational") { auto check_rational_and_integral = [](Magnitude auto m) { CHECK(is_integral(m)); CHECK(is_rational(m)); @@ -321,8 +303,7 @@ TEST_CASE("can distinguish integral, rational, and irrational magnitudes") check_rational_and_integral(as_magnitude()); } - SECTION("Fractional magnitudes are rational, but not integral") - { + SECTION ("Fractional magnitudes are rational, but not integral") { auto check_rational_but_not_integral = [](Magnitude auto m) { CHECK(!is_integral(m)); CHECK(is_rational(m)); @@ -334,8 +315,9 @@ TEST_CASE("can distinguish integral, rational, and irrational magnitudes") namespace detail { -TEST_CASE("int_power computes integer powers") { - SECTION("handles floating point") { +TEST_CASE("int_power computes integer powers") +{ + SECTION ("handles floating point") { check_same_type_and_value(int_power(0.123L, 0), 1.0L); check_same_type_and_value(int_power(0.246f, 1), 0.246f); check_same_type_and_value(int_power(0.5f, 3), 0.125f); @@ -344,7 +326,7 @@ TEST_CASE("int_power computes integer powers") { CHECK(std::is_same_v(base_power{10, 20}))>); } - SECTION("handles integral") { + SECTION ("handles integral") { check_same_type_and_value(int_power(8, 0), 1); check_same_type_and_value(int_power(9L, 1), 9L); check_same_type_and_value(int_power(2, 10), 1024); @@ -353,13 +335,13 @@ TEST_CASE("int_power computes integer powers") { TEST_CASE("Prime helper functions") { - SECTION("multiplicity") { + SECTION ("multiplicity") { CHECK(multiplicity(2, 8) == 3); CHECK(multiplicity(2, 1024) == 10); CHECK(multiplicity(11, 6655) == 3); } - SECTION("remove_power()") { + SECTION ("remove_power()") { CHECK(remove_power(17, 0, 5) == 5); CHECK(remove_power(2, 3, 24) == 3); CHECK(remove_power(11, 3, 6655) == 5); @@ -368,13 +350,11 @@ TEST_CASE("Prime helper functions") TEST_CASE("Prime factorization") { - SECTION("1 factors into the null magnitude") - { + SECTION ("1 factors into the null magnitude") { CHECK(prime_factorization_v<1> == magnitude<>{}); } - SECTION("Prime numbers factor into themselves") - { + SECTION ("Prime numbers factor into themselves") { CHECK(prime_factorization_v<2> == magnitude{}); CHECK(prime_factorization_v<3> == magnitude{}); CHECK(prime_factorization_v<5> == magnitude{}); @@ -384,29 +364,24 @@ TEST_CASE("Prime factorization") CHECK(prime_factorization_v<41> == magnitude{}); } - SECTION("Prime factorization finds factors and multiplicities") - { - CHECK(prime_factorization_v<792> == - magnitude{}); + SECTION ("Prime factorization finds factors and multiplicities") { + CHECK(prime_factorization_v<792> == magnitude{}); } } TEST_CASE("is_prime detects primes") { - SECTION("Non-positive numbers are not prime") - { + SECTION ("Non-positive numbers are not prime") { CHECK(!is_prime(-1328)); CHECK(!is_prime(-1)); CHECK(!is_prime(0)); } - SECTION("1 is not prime") - { + SECTION ("1 is not prime") { CHECK(!is_prime(1)); } - SECTION("Discriminates between primes and non-primes") - { + SECTION ("Discriminates between primes and non-primes") { CHECK(is_prime(2)); CHECK(is_prime(3)); CHECK(!is_prime(4)); @@ -422,7 +397,7 @@ TEST_CASE("is_prime detects primes") TEST_CASE("is_valid_base_power") { - SECTION("0 power is invalid") { + SECTION ("0 power is invalid") { REQUIRE(is_valid_base_power(base_power{2})); CHECK(!is_valid_base_power(base_power{2, 0})); @@ -433,7 +408,7 @@ TEST_CASE("is_valid_base_power") CHECK(!is_valid_base_power(base_power{0})); } - SECTION("non-prime integers are invalid") { + SECTION ("non-prime integers are invalid") { CHECK(!is_valid_base_power(base_power{-8})); CHECK(!is_valid_base_power(base_power{0})); CHECK(!is_valid_base_power(base_power{1})); @@ -444,7 +419,7 @@ TEST_CASE("is_valid_base_power") CHECK(!is_valid_base_power(base_power{4})); } - SECTION("non-positive floating point bases are invalid") { + SECTION ("non-positive floating point bases are invalid") { CHECK(!is_valid_base_power(base_power{})); CHECK(!is_valid_base_power(base_power{})); } @@ -452,25 +427,22 @@ TEST_CASE("is_valid_base_power") TEST_CASE("pairwise_all evaluates all pairs") { - const auto all_pairs_return_true = pairwise_all{[](auto, auto){ return true; }}; - const auto all_pairs_return_false = pairwise_all{[](auto, auto){ return false; }}; + const auto all_pairs_return_true = pairwise_all{[](auto, auto) { return true; }}; + const auto all_pairs_return_false = pairwise_all{[](auto, auto) { return false; }}; const auto all_increasing = pairwise_all{std::less{}}; - SECTION("always true for empty tuples") - { + SECTION ("always true for empty tuples") { CHECK(all_pairs_return_true()); CHECK(all_pairs_return_false()); } - SECTION("always true for single-element tuples") - { + SECTION ("always true for single-element tuples") { CHECK(all_pairs_return_true(1)); CHECK(all_pairs_return_false(3.14)); CHECK(all_pairs_return_true('x')); } - SECTION("true for longer tuples iff true for all neighbouring pairs") - { + SECTION ("true for longer tuples iff true for all neighbouring pairs") { CHECK(all_increasing(1, 1.5)); CHECK(all_increasing(1, 1.5, 2)); @@ -484,26 +456,23 @@ TEST_CASE("pairwise_all evaluates all pairs") TEST_CASE("strictly_increasing") { - SECTION("Empty input is sorted") - { + SECTION ("Empty input is sorted") { CHECK(strictly_increasing()); } - SECTION("Single-element input is sorted") - { + SECTION ("Single-element input is sorted") { CHECK(strictly_increasing(3)); CHECK(strictly_increasing(15.42)); CHECK(strictly_increasing('c')); } - SECTION("Multi-value inputs compare correctly") - { + SECTION ("Multi-value inputs compare correctly") { CHECK(strictly_increasing(3, 3.14)); CHECK(!strictly_increasing(3, 3.0)); CHECK(!strictly_increasing(4, 3.0)); } } -} // namespace detail +} // namespace detail -} // namespace units +} // namespace units diff --git a/test/unit_test/static/prime_test.cpp b/test/unit_test/static/prime_test.cpp index ab852974..26e9478d 100644 --- a/test/unit_test/static/prime_test.cpp +++ b/test/unit_test/static/prime_test.cpp @@ -21,15 +21,14 @@ // SOFTWARE. #include - #include #include -namespace units::detail -{ +namespace units::detail { template -constexpr bool check_primes(std::index_sequence) { +constexpr bool check_primes(std::index_sequence) +{ return ((Is < 2 || WheelFactorizer::is_prime(Is) == is_prime_by_trial_division(Is)) && ...); } @@ -72,4 +71,4 @@ static_assert(!WheelFactorizer<3>::is_prime(0)); static_assert(!WheelFactorizer<3>::is_prime(1)); static_assert(WheelFactorizer<3>::is_prime(2)); -} // namespace units::detail +} // namespace units::detail From a99e5f9032e2d1b98224fffccfad3b6deef6ae2e Mon Sep 17 00:00:00 2001 From: Chip Hogg Date: Sat, 19 Mar 2022 15:54:29 +0000 Subject: [PATCH 13/20] Switch tests to use top-level, anonymous namespace --- test/unit_test/runtime/magnitude_test.cpp | 12 +++++++----- test/unit_test/static/prime_test.cpp | 6 ++++-- 2 files changed, 11 insertions(+), 7 deletions(-) diff --git a/test/unit_test/runtime/magnitude_test.cpp b/test/unit_test/runtime/magnitude_test.cpp index e9b359f2..40c8c1b1 100644 --- a/test/unit_test/runtime/magnitude_test.cpp +++ b/test/unit_test/runtime/magnitude_test.cpp @@ -25,7 +25,10 @@ #include #include -namespace units { +using namespace units; +using namespace units::detail; + +namespace { // A set of non-standard bases for testing purposes. struct noninteger_base { @@ -313,7 +316,8 @@ TEST_CASE("can distinguish integral, rational, and irrational magnitudes") } } -namespace detail { +//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// +// Detail function tests below. TEST_CASE("int_power computes integer powers") { @@ -473,6 +477,4 @@ TEST_CASE("strictly_increasing") } } -} // namespace detail - -} // namespace units +} // namespace diff --git a/test/unit_test/static/prime_test.cpp b/test/unit_test/static/prime_test.cpp index 26e9478d..9cb017b6 100644 --- a/test/unit_test/static/prime_test.cpp +++ b/test/unit_test/static/prime_test.cpp @@ -24,7 +24,9 @@ #include #include -namespace units::detail { +using namespace units::detail; + +namespace { template constexpr bool check_primes(std::index_sequence) @@ -71,4 +73,4 @@ static_assert(!WheelFactorizer<3>::is_prime(0)); static_assert(!WheelFactorizer<3>::is_prime(1)); static_assert(WheelFactorizer<3>::is_prime(2)); -} // namespace units::detail +} // namespace From 166dd1e94424f6bf0d9c201ed03249c1fb16dba5 Mon Sep 17 00:00:00 2001 From: Chip Hogg Date: Sat, 19 Mar 2022 15:55:34 +0000 Subject: [PATCH 14/20] Work around numbers with very large first factors We introduce the `known_first_factor` variable template. --- src/core/include/units/magnitude.h | 35 ++++++++++++++++++++++- test/unit_test/runtime/magnitude_test.cpp | 15 +++++++++- 2 files changed, 48 insertions(+), 2 deletions(-) diff --git a/src/core/include/units/magnitude.h b/src/core/include/units/magnitude.h index 3cec3b98..ffd86e7e 100644 --- a/src/core/include/units/magnitude.h +++ b/src/core/include/units/magnitude.h @@ -28,6 +28,7 @@ #include #include #include +#include #include namespace units { @@ -282,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; @@ -468,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 = static_cast(Factorizer::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 40c8c1b1..6df348d4 100644 --- a/test/unit_test/runtime/magnitude_test.cpp +++ b/test/unit_test/runtime/magnitude_test.cpp @@ -28,6 +28,9 @@ 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. @@ -149,7 +152,17 @@ TEST_CASE("make_ratio performs prime factorization correctly") // 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>(); + as_magnitude<334'524'384'739>(); + } + + 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>(); } } From c3393838736fac7f7371873ec6ea31fb2ef430e1 Mon Sep 17 00:00:00 2001 From: Chip Hogg Date: Sat, 19 Mar 2022 16:07:39 +0000 Subject: [PATCH 15/20] Convert names to standard_case --- src/core/include/units/bits/prime.h | 2 +- src/core/include/units/magnitude.h | 6 ++-- test/unit_test/static/prime_test.cpp | 48 ++++++++++++++-------------- 3 files changed, 28 insertions(+), 28 deletions(-) diff --git a/src/core/include/units/bits/prime.h b/src/core/include/units/bits/prime.h index aaa20e4f..698f05bf 100644 --- a/src/core/include/units/bits/prime.h +++ b/src/core/include/units/bits/prime.h @@ -124,7 +124,7 @@ constexpr std::size_t product(const std::array& values) // // [1] https://en.wikipedia.org/wiki/Wheel_factorization template -struct WheelFactorizer { +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 = diff --git a/src/core/include/units/magnitude.h b/src/core/include/units/magnitude.h index ffd86e7e..94140397 100644 --- a/src/core/include/units/magnitude.h +++ b/src/core/include/units/magnitude.h @@ -34,7 +34,7 @@ namespace units { namespace detail { // Higher numbers use fewer trial divisions, at the price of more storage space. -using Factorizer = WheelFactorizer<4>; +using factorizer = wheel_factorizer<4>; } // namespace detail /** @@ -274,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 >= 0) && Factorizer::is_prime(static_cast(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) { @@ -502,7 +502,7 @@ struct prime_factorization { if constexpr (known_first_factor.has_value()) { return known_first_factor.value(); } else { - return static_cast(Factorizer::find_first_factor(N)); + return static_cast(factorizer::find_first_factor(N)); } } diff --git a/test/unit_test/static/prime_test.cpp b/test/unit_test/static/prime_test.cpp index 9cb017b6..020a1fe7 100644 --- a/test/unit_test/static/prime_test.cpp +++ b/test/unit_test/static/prime_test.cpp @@ -31,7 +31,7 @@ namespace { template constexpr bool check_primes(std::index_sequence) { - return ((Is < 2 || WheelFactorizer::is_prime(Is) == is_prime_by_trial_division(Is)) && ...); + return ((Is < 2 || wheel_factorizer::is_prime(Is) == is_prime_by_trial_division(Is)) && ...); } static_assert(check_primes<2>(std::make_index_sequence<122>{})); @@ -44,33 +44,33 @@ static_assert(check_primes<2>(std::make_index_sequence<122>{})); // 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(WheelFactorizer<4>::is_prime(109561) == is_prime_by_trial_division(109561)); +static_assert(wheel_factorizer<4>::is_prime(109561) == is_prime_by_trial_division(109561)); -static_assert(WheelFactorizer<1>::coprimes_in_first_wheel.size() == 1); -static_assert(WheelFactorizer<2>::coprimes_in_first_wheel.size() == 2); -static_assert(WheelFactorizer<3>::coprimes_in_first_wheel.size() == 8); -static_assert(WheelFactorizer<4>::coprimes_in_first_wheel.size() == 48); -static_assert(WheelFactorizer<5>::coprimes_in_first_wheel.size() == 480); +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(WheelFactorizer<3>::coprimes_in_first_wheel[0] == 1); -static_assert(WheelFactorizer<3>::coprimes_in_first_wheel[1] == 7); -static_assert(WheelFactorizer<3>::coprimes_in_first_wheel[2] == 11); -static_assert(WheelFactorizer<3>::coprimes_in_first_wheel[3] == 13); -static_assert(WheelFactorizer<3>::coprimes_in_first_wheel[4] == 17); -static_assert(WheelFactorizer<3>::coprimes_in_first_wheel[5] == 19); -static_assert(WheelFactorizer<3>::coprimes_in_first_wheel[6] == 23); -static_assert(WheelFactorizer<3>::coprimes_in_first_wheel[7] == 29); +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(!WheelFactorizer<1>::is_prime(0)); -static_assert(!WheelFactorizer<1>::is_prime(1)); -static_assert(WheelFactorizer<1>::is_prime(2)); +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(!WheelFactorizer<2>::is_prime(0)); -static_assert(!WheelFactorizer<2>::is_prime(1)); -static_assert(WheelFactorizer<2>::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(!WheelFactorizer<3>::is_prime(0)); -static_assert(!WheelFactorizer<3>::is_prime(1)); -static_assert(WheelFactorizer<3>::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 From 0f80c1010468b4a112293b62c684b9f70baf6159 Mon Sep 17 00:00:00 2001 From: Chip Hogg Date: Sat, 19 Mar 2022 16:58:43 +0000 Subject: [PATCH 16/20] Try "gentler" test case I verified that we hit GCC 10's constexpr limit with `wheel_factorizer<1>`, but pass with `wheel_factorizer<4>`. I hope this number is enough smaller than the square root of the previous value that the other compilers will be able to handle it. If not: we'll go lower. --- test/unit_test/runtime/magnitude_test.cpp | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/test/unit_test/runtime/magnitude_test.cpp b/test/unit_test/runtime/magnitude_test.cpp index c3138fb8..9958c715 100644 --- a/test/unit_test/runtime/magnitude_test.cpp +++ b/test/unit_test/runtime/magnitude_test.cpp @@ -156,15 +156,18 @@ TEST_CASE("make_ratio performs prime factorization correctly") as_magnitude(); } - SECTION ("Can handle prime number which would exceed GCC iteration limit") { + 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<334'524'384'739>(); + constexpr std::intmax_t big_prime = 414'131; + as_magnitude(); } - SECTION ("Can bypass computing primes by providing known_first_factor") { + 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. // @@ -383,7 +386,7 @@ TEST_CASE("int_power computes integer powers") TEST_CASE("Prime helper functions") { - SECTION ("multiplicity") + SECTION("multiplicity") { CHECK(multiplicity(2, 8) == 3); CHECK(multiplicity(2, 1024) == 10); From 438feb30014b9c2955f5a59ef985724311e767f7 Mon Sep 17 00:00:00 2001 From: Chip Hogg Date: Sat, 19 Mar 2022 17:05:55 +0000 Subject: [PATCH 17/20] Remove offending unit test Apparently, the constexpr depth which clang and MSVC can handle is too shallow for me to write a unit test that works on all supported compilers. --- test/unit_test/runtime/magnitude_test.cpp | 10 ---------- 1 file changed, 10 deletions(-) diff --git a/test/unit_test/runtime/magnitude_test.cpp b/test/unit_test/runtime/magnitude_test.cpp index 9958c715..9db7b696 100644 --- a/test/unit_test/runtime/magnitude_test.cpp +++ b/test/unit_test/runtime/magnitude_test.cpp @@ -156,16 +156,6 @@ TEST_CASE("make_ratio performs prime factorization correctly") 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.) - constexpr std::intmax_t big_prime = 414'131; - 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 From f495ad9e756ffea6c267da84ce51d78a1b717953 Mon Sep 17 00:00:00 2001 From: Chip Hogg Date: Sat, 19 Mar 2022 17:27:38 +0000 Subject: [PATCH 18/20] Replace `accumulate` with `reduce` Perhaps this will also satisfy Apple's Clang 13? Since `reduce` is newer, it may be more likely to be `constexpr` compatible. --- src/core/include/units/bits/prime.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/core/include/units/bits/prime.h b/src/core/include/units/bits/prime.h index 698f05bf..648dcce8 100644 --- a/src/core/include/units/bits/prime.h +++ b/src/core/include/units/bits/prime.h @@ -101,7 +101,7 @@ constexpr auto coprimes_up_to(std::size_t n, const std::array& b template constexpr std::size_t product(const std::array& values) { - return std::accumulate(std::begin(values), std::end(values), std::size_t{1u}, std::multiplies{}); + return std::reduce(values.begin(), values.end(), std::size_t{1}, std::multiplies{}); } // A configurable instantiation of the "wheel factorization" algorithm [1] for prime numbers. From 6872117bae6bf63de8ab86e89c3875f53675fa85 Mon Sep 17 00:00:00 2001 From: Chip Hogg Date: Sat, 19 Mar 2022 17:35:20 +0000 Subject: [PATCH 19/20] Replace `reduce` with bespoke implementation If _this_ isn't `constexpr` compatible, I'm going to propose removing support for the MacOS clang build. --- src/core/include/units/bits/prime.h | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/src/core/include/units/bits/prime.h b/src/core/include/units/bits/prime.h index 648dcce8..585a3b50 100644 --- a/src/core/include/units/bits/prime.h +++ b/src/core/include/units/bits/prime.h @@ -101,7 +101,11 @@ constexpr auto coprimes_up_to(std::size_t n, const std::array& b template constexpr std::size_t product(const std::array& values) { - return std::reduce(values.begin(), values.end(), std::size_t{1}, std::multiplies{}); + 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. From b589ba8d86e90871e4cf0d5f318f0cc44fb8ec80 Mon Sep 17 00:00:00 2001 From: Chip Hogg Date: Sat, 19 Mar 2022 22:09:51 +0000 Subject: [PATCH 20/20] Omit redundant computation This stems from an earlier mistake where I was using primes in the first wheel, rather than coprimes-to-the-basis. 1 is not prime, so we used to need to handle it separately (in an implementation which was, to be clear, wrong). It _is_ coprime, so now we get it for free! --- src/core/include/units/bits/prime.h | 4 ---- 1 file changed, 4 deletions(-) diff --git a/src/core/include/units/bits/prime.h b/src/core/include/units/bits/prime.h index 585a3b50..b9e06e72 100644 --- a/src/core/include/units/bits/prime.h +++ b/src/core/include/units/bits/prime.h @@ -149,10 +149,6 @@ struct wheel_factorizer { } 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 : coprimes_in_first_wheel) { if (const auto k = first_factor_maybe(n, wheel + p)) { return *k;