From 1ef2bf9f28020f29cb86a46548f4f41daf7518a6 Mon Sep 17 00:00:00 2001 From: dkavolis <12998363+dkavolis@users.noreply.github.com> Date: Mon, 2 Nov 2020 13:28:13 +0000 Subject: [PATCH] After review --- src/include/units/bits/constexpr_math.h | 48 +++++++++++++------------ src/include/units/bits/math_concepts.h | 35 ++++++++++++++++++ src/include/units/bits/root.h | 13 ++++--- src/include/units/math.h | 4 +-- src/include/units/ratio.h | 7 ++-- 5 files changed, 76 insertions(+), 31 deletions(-) create mode 100644 src/include/units/bits/math_concepts.h diff --git a/src/include/units/bits/constexpr_math.h b/src/include/units/bits/constexpr_math.h index 9f8c2849..53e4fdf4 100644 --- a/src/include/units/bits/constexpr_math.h +++ b/src/include/units/bits/constexpr_math.h @@ -23,40 +23,41 @@ #pragma once #include +#include #include #include namespace units::detail { struct decimal_fp { - double significand; - std::intmax_t mantissa; + double significant; + std::intmax_t exponent; }; [[nodiscard]] constexpr decimal_fp to_decimal(double v) noexcept { if (v == 0) { - return {.significand = 0.0, .mantissa = 0}; + return {.significant = 0.0, .exponent = 0}; } - double significand = abs(v); - std::intmax_t mantissa = 0; + double significant = abs(v); + std::intmax_t exponent = 0; - while (significand < 1) { - significand *= 10.0; - --mantissa; + while (significant < 1) { + significant *= 10.0; + --exponent; } - while (significand >= 10) { - significand /= 10.0; - ++mantissa; + while (significant >= 10) { + significant /= 10.0; + ++exponent; } if (v < 0) { - significand = -significand; + significant = -significant; } - return {.significand = significand, .mantissa = mantissa}; + return {.significant = significant, .exponent = exponent}; } /* approximate natural log as https://math.stackexchange.com/a/977836 @@ -66,8 +67,8 @@ struct decimal_fp { { 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) + // lookup table to speed up convergence for all significant values + // significant 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: @@ -177,14 +178,14 @@ struct decimal_fp { }; 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]; + // dividing the significant by nearest lower value in [1.0, 1.1, 1.2, ..., 9.9] will greatly improve convergence + x.significant *= 10; + const auto isignificant = static_cast(x.significant); + x.significant /= static_cast(isignificant); + const double result = static_cast(x.exponent - 1) * log_table[9] + log_table[isignificant - 1]; - // 1.0 <= significand < 1.1 converges rapidly - const double y = (x.significand - 1) / (x.significand + 1); + // 1.0 <= significant < 1.1 converges rapidly + const double y = (x.significant - 1) / (x.significant + 1); const double y_squared = y * y; double sum = 0; // 5 terms are needed for convergence to machine precision in the worst case scenario @@ -201,7 +202,8 @@ struct decimal_fp { 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; } + requires gt_zero +[[nodiscard]] constexpr double constexpr_exp(double v) noexcept { if constexpr (N == 0) { return 1.0; diff --git a/src/include/units/bits/math_concepts.h b/src/include/units/bits/math_concepts.h new file mode 100644 index 00000000..e92acef1 --- /dev/null +++ b/src/include/units/bits/math_concepts.h @@ -0,0 +1,35 @@ +// 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 + +namespace units::detail { + +template +concept gt_zero = (N > 0); + +template +concept non_zero = (N != 0); + +} // namespace units::detail diff --git a/src/include/units/bits/root.h b/src/include/units/bits/root.h index c8c18691..cf064a60 100644 --- a/src/include/units/bits/root.h +++ b/src/include/units/bits/root.h @@ -24,6 +24,7 @@ #include #include +#include #include #include #include @@ -31,7 +32,8 @@ namespace units::detail { template -[[nodiscard]] constexpr std::intmax_t iroot_impl(std::intmax_t v, F const& pow_function) noexcept requires requires { N > 0; } + requires gt_zero +[[nodiscard]] constexpr std::intmax_t iroot_impl(std::intmax_t v, F const& pow_function) noexcept { if constexpr (N == 1) { return v; @@ -61,13 +63,15 @@ template // 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; } + requires gt_zero +[[nodiscard]] constexpr std::intmax_t iroot_compile(std::intmax_t v) noexcept { 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; } + requires gt_zero +[[nodiscard]] std::intmax_t iroot_runtime(std::intmax_t v) noexcept { return iroot_impl(v, [](double x, double exponent) { if constexpr (N == 2) { @@ -81,7 +85,8 @@ template } template -[[nodiscard]] constexpr std::intmax_t iroot(std::intmax_t v) noexcept requires requires { N > 0; } + requires gt_zero +[[nodiscard]] constexpr std::intmax_t iroot(std::intmax_t v) noexcept { // compile time version is much slower, use faster version at runtime if (std::is_constant_evaluated()) { diff --git a/src/include/units/math.h b/src/include/units/math.h index 3d5bc059..2c40e8c2 100644 --- a/src/include/units/math.h +++ b/src/include/units/math.h @@ -40,8 +40,8 @@ namespace units { * @return Quantity The result of computation */ template -[[nodiscard]] inline auto pow(const Q& q) noexcept - requires requires { std::pow(q.count(), 1.0); Den != 0; } + requires detail::non_zero +[[nodiscard]] inline auto pow(const Q& q) noexcept requires requires { std::pow(q.count(), 1.0); } { using rep = TYPENAME Q::rep; if constexpr (Num == 0) { diff --git a/src/include/units/ratio.h b/src/include/units/ratio.h index 71ac79dc..2af76f9b 100644 --- a/src/include/units/ratio.h +++ b/src/include/units/ratio.h @@ -23,6 +23,7 @@ #pragma once #include +#include #include #include #include @@ -107,7 +108,8 @@ namespace detail { } template -[[nodiscard]] constexpr ratio root(const ratio& r) requires requires { N > 0; } + requires gt_zero +[[nodiscard]] constexpr ratio root(const ratio& r) { if constexpr (N == 1) { return r; @@ -124,7 +126,8 @@ template } // namespace detail template -[[nodiscard]] constexpr ratio pow(const ratio& r) requires requires { Den != 0; } + requires detail::non_zero +[[nodiscard]] constexpr ratio pow(const ratio& r) { if constexpr (Num == 0) { return ratio(1);