From 23dd21465dbf458843b44478facdf2934d49878d Mon Sep 17 00:00:00 2001 From: dkavolis <12998363+dkavolis@users.noreply.github.com> Date: Sat, 17 Oct 2020 20:50:58 +0100 Subject: [PATCH] rational powers --- src/include/units/bits/constexpr_math.h | 239 ++++++++++++++++++++++++ src/include/units/bits/dimension_op.h | 66 ++----- src/include/units/bits/pow.h | 16 ++ src/include/units/bits/root.h | 94 ++++++++++ src/include/units/math.h | 51 +++-- src/include/units/ratio.h | 94 +++++----- test/unit_test/runtime/math_test.cpp | 81 ++++++++ test/unit_test/static/math_test.cpp | 11 ++ test/unit_test/static/ratio_test.cpp | 11 +- 9 files changed, 555 insertions(+), 108 deletions(-) create mode 100644 src/include/units/bits/constexpr_math.h create mode 100644 src/include/units/bits/root.h diff --git a/src/include/units/bits/constexpr_math.h b/src/include/units/bits/constexpr_math.h new file mode 100644 index 00000000..9f8c2849 --- /dev/null +++ b/src/include/units/bits/constexpr_math.h @@ -0,0 +1,239 @@ +// 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 + +namespace units::detail { + +struct decimal_fp { + double significand; + std::intmax_t mantissa; +}; + +[[nodiscard]] constexpr decimal_fp to_decimal(double v) noexcept +{ + if (v == 0) { + return {.significand = 0.0, .mantissa = 0}; + } + + double significand = abs(v); + std::intmax_t mantissa = 0; + + while (significand < 1) { + significand *= 10.0; + --mantissa; + } + + while (significand >= 10) { + significand /= 10.0; + ++mantissa; + } + + if (v < 0) { + significand = -significand; + } + + return {.significand = significand, .mantissa = mantissa}; +} + +/* approximate natural log as https://math.stackexchange.com/a/977836 + far slower than std::log but works at compile time with similar accuracy + */ +[[nodiscard]] constexpr double constexpr_log(double v) noexcept +{ + Expects(v > 0); + + // lookup table to speed up convergence for all significand values + // significand values of 7 and greater benefit mostly as they now converge in 5 terms compared to O(10)-O(100) + // required without the table + // + // using python: + // >>> import math + // >>> for i in range(1, 100): + // ... print(f"/* log({i:>2d}) = */ {math.log(i):.16f},") + constexpr std::array log_table{ + /* log( 1) = */ 0.0000000000000000, + /* log( 2) = */ 0.6931471805599453, + /* log( 3) = */ 1.0986122886681098, + /* log( 4) = */ 1.3862943611198906, + /* log( 5) = */ 1.6094379124341003, + /* log( 6) = */ 1.7917594692280550, + /* log( 7) = */ 1.9459101490553132, + /* log( 8) = */ 2.0794415416798357, + /* log( 9) = */ 2.1972245773362196, + /* log(10) = */ 2.3025850929940459, + /* log(11) = */ 2.3978952727983707, + /* log(12) = */ 2.4849066497880004, + /* log(13) = */ 2.5649493574615367, + /* log(14) = */ 2.6390573296152584, + /* log(15) = */ 2.7080502011022101, + /* log(16) = */ 2.7725887222397811, + /* log(17) = */ 2.8332133440562162, + /* log(18) = */ 2.8903717578961645, + /* log(19) = */ 2.9444389791664403, + /* log(20) = */ 2.9957322735539909, + /* log(21) = */ 3.0445224377234230, + /* log(22) = */ 3.0910424533583161, + /* log(23) = */ 3.1354942159291497, + /* log(24) = */ 3.1780538303479458, + /* log(25) = */ 3.2188758248682006, + /* log(26) = */ 3.2580965380214821, + /* log(27) = */ 3.2958368660043291, + /* log(28) = */ 3.3322045101752038, + /* log(29) = */ 3.3672958299864741, + /* log(30) = */ 3.4011973816621555, + /* log(31) = */ 3.4339872044851463, + /* log(32) = */ 3.4657359027997265, + /* log(33) = */ 3.4965075614664802, + /* log(34) = */ 3.5263605246161616, + /* log(35) = */ 3.5553480614894135, + /* log(36) = */ 3.5835189384561099, + /* log(37) = */ 3.6109179126442243, + /* log(38) = */ 3.6375861597263857, + /* log(39) = */ 3.6635616461296463, + /* log(40) = */ 3.6888794541139363, + /* log(41) = */ 3.7135720667043080, + /* log(42) = */ 3.7376696182833684, + /* log(43) = */ 3.7612001156935624, + /* log(44) = */ 3.7841896339182610, + /* log(45) = */ 3.8066624897703196, + /* log(46) = */ 3.8286413964890951, + /* log(47) = */ 3.8501476017100584, + /* log(48) = */ 3.8712010109078911, + /* log(49) = */ 3.8918202981106265, + /* log(50) = */ 3.9120230054281460, + /* log(51) = */ 3.9318256327243257, + /* log(52) = */ 3.9512437185814275, + /* log(53) = */ 3.9702919135521220, + /* log(54) = */ 3.9889840465642745, + /* log(55) = */ 4.0073331852324712, + /* log(56) = */ 4.0253516907351496, + /* log(57) = */ 4.0430512678345503, + /* log(58) = */ 4.0604430105464191, + /* log(59) = */ 4.0775374439057197, + /* log(60) = */ 4.0943445622221004, + /* log(61) = */ 4.1108738641733114, + /* log(62) = */ 4.1271343850450917, + /* log(63) = */ 4.1431347263915326, + /* log(64) = */ 4.1588830833596715, + /* log(65) = */ 4.1743872698956368, + /* log(66) = */ 4.1896547420264252, + /* log(67) = */ 4.2046926193909657, + /* log(68) = */ 4.2195077051761070, + /* log(69) = */ 4.2341065045972597, + /* log(70) = */ 4.2484952420493594, + /* log(71) = */ 4.2626798770413155, + /* log(72) = */ 4.2766661190160553, + /* log(73) = */ 4.2904594411483910, + /* log(74) = */ 4.3040650932041702, + /* log(75) = */ 4.3174881135363101, + /* log(76) = */ 4.3307333402863311, + /* log(77) = */ 4.3438054218536841, + /* log(78) = */ 4.3567088266895917, + /* log(79) = */ 4.3694478524670215, + /* log(80) = */ 4.3820266346738812, + /* log(81) = */ 4.3944491546724391, + /* log(82) = */ 4.4067192472642533, + /* log(83) = */ 4.4188406077965983, + /* log(84) = */ 4.4308167988433134, + /* log(85) = */ 4.4426512564903167, + /* log(86) = */ 4.4543472962535073, + /* log(87) = */ 4.4659081186545837, + /* log(88) = */ 4.4773368144782069, + /* log(89) = */ 4.4886363697321396, + /* log(90) = */ 4.4998096703302650, + /* log(91) = */ 4.5108595065168497, + /* log(92) = */ 4.5217885770490405, + /* log(93) = */ 4.5325994931532563, + /* log(94) = */ 4.5432947822700038, + /* log(95) = */ 4.5538768916005408, + /* log(96) = */ 4.5643481914678361, + /* log(97) = */ 4.5747109785033828, + /* log(98) = */ 4.5849674786705723, + /* log(99) = */ 4.5951198501345898, + }; + decimal_fp x = to_decimal(v); + + // dividing the significand by nearest lower value in [1.0, 1.1, 1.2, ..., 9.9] will greatly improve convergence + x.significand *= 10; + const auto isignificand = static_cast(x.significand); + x.significand /= static_cast(isignificand); + const double result = static_cast(x.mantissa - 1) * log_table[9] + log_table[isignificand - 1]; + + // 1.0 <= significand < 1.1 converges rapidly + const double y = (x.significand - 1) / (x.significand + 1); + const double y_squared = y * y; + double sum = 0; + // 5 terms are needed for convergence to machine precision in the worst case scenario + for (int k = 4; k > 0; --k) { + sum = y_squared * (1 / (2 * static_cast(k) + 1) + sum); + } + sum = 2 * y * (1 + sum); // k = 0 term + return result + sum; +} + +/* approximate e^x as Taylor series e^x = 1 + x/1! + x^2/2! + x^3/3! +... where N is the order of the Taylor series + use https://math.stackexchange.com/a/1988927 to improve convergence for large values + + larger Factor values improve convergence for all values but reduce the precision +*/ +template +[[nodiscard]] constexpr double constexpr_exp(double v) noexcept requires requires { Factor > 0; } +{ + if constexpr (N == 0) { + return 1.0; + } else { + constexpr auto coefficients = []() { + std::array coeffs; + std::size_t factorial = 1; + for (std::size_t i = 0; i < N; ++i) { + factorial *= i + 1; + coeffs[i] = 1.0 / static_cast(factorial); + } + return coeffs; + }(); + + const double x = v / static_cast(Factor); + double result = 0; + for (auto i = static_cast(N - 1); i >= 0; --i) { + result = x * (coefficients[static_cast(i)] + result); + } + + // for factors of power of 2 this should be replaced by log2(Factor) multiplications by the compiler + return pow_impl(1 + result); + } +} + +// default template arguments provide reasonable precision even for fairly large exponents +// see constexpr_exp for template arguments +template +[[nodiscard]] constexpr double constexpr_pow(double v, double exponent) noexcept +{ + const double x = exponent * constexpr_log(v); + return constexpr_exp(x); +} + +} // namespace units::detail diff --git a/src/include/units/bits/dimension_op.h b/src/include/units/bits/dimension_op.h index 608d00f1..e12c7383 100644 --- a/src/include/units/bits/dimension_op.h +++ b/src/include/units/bits/dimension_op.h @@ -159,71 +159,43 @@ using dimension_multiply = TYPENAME detail::dimension_multiply_impl::typ template using dimension_divide = TYPENAME detail::dimension_multiply_impl>::type; -// dimension_sqrt -namespace detail { - -template -struct dimension_sqrt_impl; - -template -struct dimension_sqrt_impl { - using type = downcast_dimension>>; -}; - -template -struct dimension_sqrt_impl>> { - using type = D; -}; - -template -struct dimension_sqrt_impl { - using type = TYPENAME dimension_sqrt_impl::type; -}; - -template -struct dimension_sqrt_impl> { - using type = downcast_dimension...>>; -}; - -} // namespace detail - -template -using dimension_sqrt = TYPENAME detail::dimension_sqrt_impl::type; - // dimension_pow namespace detail { -template +template struct dimension_pow_impl; -template -struct dimension_pow_impl { - using type = downcast_dimension>>; +template +struct dimension_pow_impl { + using type = downcast_dimension>>; }; template -struct dimension_pow_impl { +struct dimension_pow_impl { using type = D; }; -template -struct dimension_pow_impl>, N> { +template +struct dimension_pow_impl>, Num, Den> { using type = D; }; -template -struct dimension_pow_impl { - using type = TYPENAME dimension_pow_impl, N>::type; +template +struct dimension_pow_impl { + using type = TYPENAME dimension_pow_impl, Num, Den>::type; }; -template -struct dimension_pow_impl, N> { - using type = downcast_dimension...>>; -}; +template +struct dimension_pow_impl, Num, Den> { + using type = downcast_dimension...>>; +}; } // namespace detail -template -using dimension_pow = TYPENAME detail::dimension_pow_impl::type; +template +using dimension_pow = TYPENAME detail::dimension_pow_impl::type; + +template +using dimension_sqrt = TYPENAME detail::dimension_pow_impl::type; } // namespace units diff --git a/src/include/units/bits/pow.h b/src/include/units/bits/pow.h index c6693103..28b2ab8f 100644 --- a/src/include/units/bits/pow.h +++ b/src/include/units/bits/pow.h @@ -58,4 +58,20 @@ constexpr Rep fpow10(std::intmax_t exp) return result; } +template +constexpr T pow_impl(const T& v) noexcept +{ + if constexpr (N == 0) { + return T(1); + } else if constexpr (N == 1) { + return v; + } else if constexpr (N < 0) { + return 1 / pow_impl<-N>(v); + } else if constexpr (N % 2 == 0) { // even + return pow_impl(v * v); + } else { // odd + return v * pow_impl<(N - 1) / 2>(v * v); + } +} + } // namespace units::detail diff --git a/src/include/units/bits/root.h b/src/include/units/bits/root.h new file mode 100644 index 00000000..c8c18691 --- /dev/null +++ b/src/include/units/bits/root.h @@ -0,0 +1,94 @@ +// The MIT License (MIT) +// +// Copyright (c) 2018 Mateusz Pusz +// +// Permission is hereby granted, free of charge, to any person obtaining a copy +// of this software and associated documentation files (the "Software"), to deal +// in the Software without restriction, including without limitation the rights +// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +// copies of the Software, and to permit persons to whom the Software is +// furnished to do so, subject to the following conditions: +// +// The above copyright notice and this permission notice shall be included in all +// copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +// SOFTWARE. + +#pragma once + +#include +#include +#include +#include +#include + +namespace units::detail { + +template +[[nodiscard]] constexpr std::intmax_t iroot_impl(std::intmax_t v, F const& pow_function) noexcept requires requires { N > 0; } +{ + if constexpr (N == 1) { + return v; + } else { + Expects(v >= 0); + if (v == 0) { + return 0; + } + + constexpr double exponent = 1.0 / static_cast(N); + const auto root = static_cast(pow_function(static_cast(v), exponent)); + + // integer roots may be truncated down by 1 or in even rarer cases up by 1 due to finite precision of pow and + // exponent, check both cases + if (v == pow_impl(root + 1)) { + return root + 1; + } + if (v < pow_impl(root)) { + return root - 1; + } + return root; + } +} + +// maximum v is std::numeric_limits::max() which is the worst case for exp convergence +// ExpOrder = 12 and Factor = 64 give a precision of about O(1e-15) for a wide range of 1 / N exponents +// Factor = 32 needs quite a few more terms to converge +// https://godbolt.org/z/odWq1o +template +[[nodiscard]] constexpr std::intmax_t iroot_compile(std::intmax_t v) noexcept requires requires { N > 0; } +{ + return iroot_impl(v, [](double x, double exponent) { return constexpr_pow(x, exponent); }); +} + +template +[[nodiscard]] std::intmax_t iroot_runtime(std::intmax_t v) noexcept requires requires { N > 0; } +{ + return iroot_impl(v, [](double x, double exponent) { + if constexpr (N == 2) { + return std::sqrt(x); + } else if constexpr (N == 3) { + return std::cbrt(x); + } else { + return std::pow(x, exponent); + } + }); +} + +template +[[nodiscard]] constexpr std::intmax_t iroot(std::intmax_t v) noexcept requires requires { N > 0; } +{ + // compile time version is much slower, use faster version at runtime + if (std::is_constant_evaluated()) { + return iroot_compile(v); + } + + return iroot_runtime(v); +} + +} // namespace units::detail diff --git a/src/include/units/math.h b/src/include/units/math.h index 534c0c0a..3d5bc059 100644 --- a/src/include/units/math.h +++ b/src/include/units/math.h @@ -31,46 +31,65 @@ namespace units { /** * @brief Computes the value of a quantity raised to the power `N` - * + * * Both the quantity value and its dimension are the base of the operation. - * - * @tparam N Exponent + * + * @tparam Num Exponent numerator + * @tparam Den Exponent denominator * @param q Quantity being the base of the operation - * @return Quantity The result of computation + * @return Quantity The result of computation */ -template +template [[nodiscard]] inline auto pow(const Q& q) noexcept - requires requires { std::pow(q.count(), N); } + requires requires { std::pow(q.count(), 1.0); Den != 0; } { using rep = TYPENAME Q::rep; - if constexpr(N == 0) { + if constexpr (Num == 0) { return rep(1); - } - else { - using dim = dimension_pow; - using unit = downcast_unit(Q::unit::ratio)>; - return quantity(static_cast(std::pow(q.count(), N))); + } else { + using dim = dimension_pow; + using unit = downcast_unit(Q::unit::ratio)>; + return quantity( + static_cast(std::pow(q.count(), static_cast(Num) / static_cast(Den)))); } } /** * @brief Computes the square root of a quantity - * + * * Both the quantity value and its dimension are the base of the operation. - * + * * @param q Quantity being the base of the operation - * @return Quantity The result of computation + * @return Quantity The result of computation */ template [[nodiscard]] inline Quantity auto sqrt(const Q& q) noexcept requires requires { std::sqrt(q.count()); } { - using dim = dimension_sqrt; + using dim = dimension_pow; using unit = downcast_unit; using rep = TYPENAME Q::rep; return quantity(static_cast(std::sqrt(q.count()))); } +/** + * @brief Computes the cubic root of a quantity + * + * Both the quantity value and its dimension are the base of the operation. + * + * @param q Quantity being the base of the operation + * @return Quantity The result of computation + */ +template +[[nodiscard]] inline Quantity auto cbrt(const Q& q) noexcept + requires requires { std::cbrt(q.count()); } +{ + using dim = dimension_pow; + using unit = downcast_unit; + using rep = TYPENAME Q::rep; + return quantity(static_cast(std::cbrt(q.count()))); +} + /** * @brief Computes Euler's raised to the given power * diff --git a/src/include/units/ratio.h b/src/include/units/ratio.h index 9848b142..71ac79dc 100644 --- a/src/include/units/ratio.h +++ b/src/include/units/ratio.h @@ -24,6 +24,8 @@ #include #include +#include +#include #include #include #include @@ -37,7 +39,7 @@ constexpr ratio inverse(const ratio& r); /** * @brief Provides compile-time rational arithmetic support. - * + * * This class is really similar to @c std::ratio but gets an additional `Exp` * template parameter that defines the exponent of the ratio. Another important * difference is the fact that the objects of that class are used as class NTTPs @@ -85,64 +87,70 @@ struct ratio { } } -template -[[nodiscard]] constexpr ratio pow(const ratio& r) -{ - if constexpr(N == 0) - return ratio(1); - else if constexpr(N == 1) - return r; - else - return pow(r) * r; -} - namespace detail { -// sqrt_impl avoids overflow and recursion -// from http://www.codecodex.com/wiki/Calculate_an_integer_square_root#C.2B.2B -// if v >= place this will fail (so we can't quite use the last bit) -[[nodiscard]] constexpr std::intmax_t sqrt_impl(std::intmax_t v) +[[nodiscard]] constexpr auto make_exp_align(const ratio& r, std::intmax_t alignment) { - // place = 0x4000 0000 for 32bit - // place = 0x4000 0000 0000 0000 for 64bit - std::intmax_t place = static_cast(1) << (sizeof(std::intmax_t) * 8 - 2); - while (place > v) place /= 4; // optimized by complier as place >>= 2 + Expects(alignment > 0); + const std::intmax_t rem = r.exp % alignment; - std::intmax_t root = 0; - while (place) { - if (v >= root + place) { - v -= root + place; - root += place * 2; - } - root /= 2; - place /= 4; + if (rem == 0) { // already aligned + return std::array{r.num, r.den, r.exp}; } - return root; + + if (r.exp > 0) { // remainder is positive + return std::array{r.num * ipow10(rem), r.den, r.exp - rem}; + } + + // remainder is negative + return std::array{r.num, r.den * ipow10(-rem), r.exp - rem}; } -[[nodiscard]] constexpr auto make_exp_even(const ratio& r) +template +[[nodiscard]] constexpr ratio root(const ratio& r) requires requires { N > 0; } { - if(r.exp % 2 == 0) - return std::array{r.num, r.den, r.exp}; // already even (incl zero) + if constexpr (N == 1) { + return r; + } else { + if (r.num == 0) { + return ratio(0); + } - // safely make exp even, so it can be divided by 2 for square root - if(r.exp > 0) - return std::array{r.num * 10, r.den, r.exp - 1}; - else - return std::array{r.num, r.den * 10, r.exp + 1}; + const auto aligned = make_exp_align(r, N); + return ratio(iroot(aligned[0]), iroot(aligned[1]), aligned[2] / N); + } } } // namespace detail -[[nodiscard]] constexpr ratio sqrt(const ratio& r) +template +[[nodiscard]] constexpr ratio pow(const ratio& r) requires requires { Den != 0; } { - if(r.num == 0) - return ratio(0); - - const auto even = detail::make_exp_even(r); - return ratio(detail::sqrt_impl(even[0]), detail::sqrt_impl(even[1]), even[2] / 2); + if constexpr (Num == 0) { + return ratio(1); + } else if constexpr (Num == Den) { + return r; + } else { + // simplify factors first and compute power for positive exponent + constexpr std::intmax_t gcd = std::gcd(Num, Den); + constexpr std::intmax_t num = detail::abs(Num / gcd); + constexpr std::intmax_t den = detail::abs(Den / gcd); + + // integer root loses precision so do pow first + const ratio result = detail::root(detail::pow_impl(r)); + + if constexpr (Num * Den < 0) { // account for negative exponent + return inverse(result); + } else { + return result; + } + } } +[[nodiscard]] constexpr ratio sqrt(const ratio& r) { return pow<1, 2>(r); } + +[[nodiscard]] constexpr ratio cbrt(const ratio& r) { return pow<1, 3>(r); } + // common_ratio [[nodiscard]] constexpr ratio common_ratio(const ratio& r1, const ratio& r2) { diff --git a/test/unit_test/runtime/math_test.cpp b/test/unit_test/runtime/math_test.cpp index 2dad7f1f..04f8a597 100644 --- a/test/unit_test/runtime/math_test.cpp +++ b/test/unit_test/runtime/math_test.cpp @@ -54,6 +54,16 @@ TEST_CASE("'sqrt()' on quantity changes the value and the dimension accordingly" REQUIRE(sqrt(4_q_m2) == 2_q_m); } +TEST_CASE("'cbrt()' on quantity changes the value and the dimension accordingly", "[math][cbrt]") +{ + REQUIRE(cbrt(8_q_m3) == 2_q_m); +} + +TEST_CASE("'pow()' on quantity changes the value and the dimension accordingly", "[math][pow]") +{ + REQUIRE(pow<1, 4>(16_q_m2) == sqrt(4_q_m)); +} + TEST_CASE("absolute functions on quantity returns the absolute value", "[math][abs][fabs]") { SECTION ("'abs()' on a negative quantity returns the abs") @@ -99,3 +109,74 @@ TEST_CASE("numeric_limits functions", "[limits]") REQUIRE(epsilon().count() != std::numeric_limits::epsilon()); } } + +TEMPLATE_TEST_CASE_SIG("pow() implementation exponentiates values to power N", "[math][pow][exp]", + (std::intmax_t N, N), 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 15, 20, 25) +{ + auto v = GENERATE(range(0.5, 20.0, 0.5)); + REQUIRE(detail::pow_impl(v) == Approx(std::pow(v, N)).epsilon(1e-15).margin(0)); +} + +template +struct Pow { + constexpr static std::intmax_t exponent = N; + template + [[nodiscard]] constexpr static auto pow(T const& v) noexcept + { + return detail::pow_impl(v); + } +}; + +template +struct CompileRoot : Pow { + [[nodiscard]] constexpr static std::intmax_t root(std::intmax_t v) noexcept { return detail::iroot_compile(v); } +}; + +template +struct RuntimeRoot : Pow { + [[nodiscard]] static std::intmax_t root(std::intmax_t v) noexcept { return detail::iroot_runtime(v); } +}; + +// test to make sure precision is not lost when rounding what should be integer roots +template +static void root_test() +{ + SECTION ("Roots are truncated down") { + auto base = GENERATE(range(1.0, 10.0, 1.0)); // doubles to guard against overflow + + if (TestType::pow(base) < static_cast(std::numeric_limits::max())) { + const std::intmax_t x = TestType::pow(static_cast(base)); + const auto expect = static_cast(base); + + REQUIRE(TestType::root(x - 1) == expect - 1); + REQUIRE(TestType::root(x) == expect); + } + } + + SECTION ("Roots are truncated correctly for very large inputs") { + auto exponent = GENERATE(range(10, std::numeric_limits::digits10, 1)); + const auto large_val = static_cast(std::pow(10, exponent)); + const auto expected = static_cast(std::pow(10, exponent / static_cast(TestType::exponent))); + + REQUIRE(TestType::root(large_val) == expected); + } +} + +/* Catch2 uses int for indexing in TEMPLATE_PRODUCT_TEST_CASE_SIG so it does not compile with -Werror=sign-conversion + * https://github.com/catchorg/Catch2/pull/2074 + +TEMPLATE_PRODUCT_TEST_CASE_SIG("detail::iroot()", "[math][pow][iroot]", (std::intmax_t N, N), + (CompileRoot, RuntimeRoot), (1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 15, 20, 25)) +{ + root_test(); +} + */ +#define ROOT_TEST_CASE(Type) \ + TEMPLATE_TEST_CASE_SIG("detail::iroot() - " #Type, "[math][pow][iroot]", (std::intmax_t N, N), 1, 2, 3, 4, 5, 6, \ + 7, 8, 9, 10, 15, 20, 25) \ + { \ + root_test>(); \ + } + +ROOT_TEST_CASE(CompileRoot) +ROOT_TEST_CASE(RuntimeRoot) diff --git a/test/unit_test/static/math_test.cpp b/test/unit_test/static/math_test.cpp index 3fb2b3ce..5de148f7 100644 --- a/test/unit_test/static/math_test.cpp +++ b/test/unit_test/static/math_test.cpp @@ -39,5 +39,16 @@ static_assert(compare(2_q_ft)), decltype(4_q_ft2)>); static_assert(compare); static_assert(compare); static_assert(compare); +static_assert(compare); +static_assert(compare); +static_assert(compare); +static_assert(compare(4_q_m2 * 4_q_m2)), decltype(2_q_m)>); +static_assert(compare(4_q_km2 * 4_q_km2)), decltype(2_q_km)>); +static_assert(compare(4_q_ft2 * 4_q_ft2)), decltype(2_q_ft)>); + +// rational dimensions +static_assert(compare(4_q_m2)), decltype(sqrt(2_q_m))>); +static_assert(compare(4_q_km2)), decltype(sqrt(2_q_km))>); +static_assert(compare(4_q_ft2)), decltype(sqrt(2_q_ft))>); } // namespace diff --git a/test/unit_test/static/ratio_test.cpp b/test/unit_test/static/ratio_test.cpp index 4284439c..3462c483 100644 --- a/test/unit_test/static/ratio_test.cpp +++ b/test/unit_test/static/ratio_test.cpp @@ -63,18 +63,25 @@ static_assert(pow<2>(ratio(1, 2)) == ratio(1, 4)); static_assert(pow<3>(ratio(1, 2)) == ratio(1, 8)); // pow with exponents -static_assert(pow<2>(ratio(1, 2, 3)) == ratio(1, 4, 6)); +static_assert(pow<2>(ratio(1, 2, 3)) == ratio(1, 4, 6)); +static_assert(pow<4, 2>(ratio(1, 2, 3)) == ratio(1, 4, 6)); static_assert(pow<3>(ratio(1, 2, -6)) == ratio(1, 8, -18)); static_assert(sqrt(ratio(9)) == ratio(3)); +static_assert(cbrt(ratio(27)) == ratio(3)); static_assert(sqrt(ratio(4)) == ratio(2)); +static_assert(cbrt(ratio(8)) == ratio(2)); static_assert(sqrt(ratio(1)) == ratio(1)); +static_assert(cbrt(ratio(1)) == ratio(1)); static_assert(sqrt(ratio(0)) == ratio(0)); +static_assert(cbrt(ratio(0)) == ratio(0)); static_assert(sqrt(ratio(1, 4)) == ratio(1, 2)); +static_assert(cbrt(ratio(1, 8)) == ratio(1, 2)); // sqrt with exponents static_assert(sqrt(ratio(9, 1, 2)) == ratio(3, 1, 1)); -static_assert(sqrt(ratio(4)) == ratio(2)); +static_assert(cbrt(ratio(27, 1, 3)) == ratio(3, 1, 1)); +static_assert(cbrt(ratio(27, 1, 2)) == ratio(13, 1, 0)); // common_ratio static_assert(common_ratio(ratio(1), ratio(1000)) == ratio(1));