diff --git a/src/include/units/bits/ratio_maths.h b/src/include/units/bits/ratio_maths.h new file mode 100644 index 00000000..b11cdeea --- /dev/null +++ b/src/include/units/bits/ratio_maths.h @@ -0,0 +1,173 @@ +// The MIT License (MIT) +// +// Copyright (c) 2018 Mateusz Pusz, Conor Williams, Oliver Schonrock +// +// Permission is hereby granted, free of charge, to any person obtaining a copy +// of this software and associated documentation files (the "Software"), to deal +// in the Software without restriction, including without limitation the rights +// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +// copies of the Software, and to permit persons to whom the Software is +// furnished to do so, subject to the following conditions: +// +// The above copyright notice and this permission notice shall be included in all +// copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +// SOFTWARE. + +#pragma once + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace units::detail { + +template +[[nodiscard]] constexpr T abs(T v) noexcept +{ + return v < 0 ? -v : v; +} + +// the following functions enable gcd and related computations on ratios +// with exponents. They avoid overflow. Further information here: +// https://github.com/mpusz/units/issues/62#issuecomment-588152833 + +// Computes (a * b) mod m relies on unsigned integer arithmetic, should not +// overflow +constexpr std::uint64_t mulmod(std::uint64_t a, std::uint64_t b, std::uint64_t m) +{ + std::uint64_t res = 0; + + if (b >= m) { + if (m > UINT64_MAX / 2u) { + b -= m; + } else { + b %= m; + } + } + + while (a != 0) { + if (a & 1) { + if (b >= m - res) { + res -= m; + } + res += b; + } + a >>= 1; + + std::uint64_t temp_b = b; + if (b >= m - b) { + temp_b -= m; + } + b += temp_b; + } + + return res; +} + +// Calculates (a ^ e) mod m , should not overflow. +constexpr std::uint64_t modpow(std::uint64_t a, std::uint64_t e, std::uint64_t m) +{ + a %= m; + std::uint64_t result = 1; + + while (e > 0) { + if (e & 1) { + result = mulmod(result, a, m); + } + a = mulmod(a, a, m); + e >>= 1; + } + return result; +} + +// gcd(a * 10 ^ e, b), should not overflow +constexpr std::intmax_t gcdpow(std::intmax_t a, std::intmax_t e, std::intmax_t b) noexcept +{ + assert(a > 0); + assert(e >= 0); + assert(b > 0); + + // gcd(i, j) = gcd(j, i mod j) for j != 0 Euclid; + // + // gcd(a 10^e, b) = gcd(b, a 10^e mod b) + // + // (a 10^e) mod b -> [ (a mod b) (10^e mod b) ] mod b + + return std::gcd( + b, static_cast(mulmod(static_cast(a % b), + modpow(10, static_cast(e), static_cast(b)), + static_cast(b)))); +} + +constexpr void cwap(std::intmax_t& lhs, std::intmax_t& rhs) +{ + std::intmax_t tmp = lhs; + lhs = rhs; + rhs = tmp; +} + +// Computes the rational gcd of n1/d1 x 10^e1 and n2/d2 x 10^e2 +constexpr auto gcd_frac(std::intmax_t n1, std::intmax_t d1, std::intmax_t e1, std::intmax_t n2, std::intmax_t d2, + std::intmax_t e2) noexcept +{ + // Short cut for equal ratios + if (n1 == n2 && d1 == d2 && e1 == e2) { + return std::array{n1, d1, e1}; + } + + if (e2 > e1) { + detail::cwap(n1, n2); + detail::cwap(d1, d2); + detail::cwap(e1, e2); + } + + std::intmax_t exp = e2; // minimum + + // gcd(a/b,c/d) = gcd(a⋅d, c⋅b) / b⋅d + + assert(std::numeric_limits::max() / n1 > d2); + assert(std::numeric_limits::max() / n2 > d1); + + std::intmax_t num = detail::gcdpow(n1 * d2, e1 - e2, n2 * d1); + + assert(std::numeric_limits::max() / d1 > d2); + + std::intmax_t den = d1 * d2; + + std::intmax_t gcd = std::gcd(num, den); + + return std::array{num / gcd, den / gcd, exp}; +} + +constexpr auto normalize(std::intmax_t num, std::intmax_t den, std::intmax_t exp) +{ + std::intmax_t gcd = std::gcd(num, den); + num = num * (den < 0 ? -1 : 1) / gcd; + den = detail::abs(den) / gcd; + + while (num % 10 == 0) { + num /= 10; + ++exp; + } + while (den % 10 == 0) { + den /= 10; + --exp; + } + + return std::array{num, den, exp}; +} + +} // namespace units::detail diff --git a/src/include/units/quantity.h b/src/include/units/quantity.h index 17f9bb69..ba26180f 100644 --- a/src/include/units/quantity.h +++ b/src/include/units/quantity.h @@ -46,7 +46,7 @@ concept safe_convertible = // exposition only template concept safe_divisible = // exposition only treat_as_floating_point || - ratio_divide::den == 1; + ratio_divide::is_integral(); } // namespace detail diff --git a/src/include/units/ratio.h b/src/include/units/ratio.h index 6d2321a4..f5a9bdcb 100644 --- a/src/include/units/ratio.h +++ b/src/include/units/ratio.h @@ -24,6 +24,7 @@ #include #include +#include #include #include #include @@ -32,31 +33,6 @@ namespace units { -namespace detail { - -template -[[nodiscard]] constexpr T abs(T v) noexcept { - return v < 0 ? -v : v; -} - -constexpr std::tuple normalize(std::intmax_t num, std::intmax_t den, std::intmax_t exp) { - std::intmax_t gcd = std::gcd(num, den); - num = num * (den < 0 ? -1 : 1) / gcd; - den = detail::abs(den) / gcd; - - while (num % 10 == 0) { - num /= 10; - ++exp; - } - while (den % 10 == 0) { - den /= 10; - --exp; - } - return std::make_tuple(num, den, exp); -} - -} // namespace detail - template requires(Den != 0) struct ratio { @@ -67,11 +43,19 @@ struct ratio { static constexpr auto norm = detail::normalize(Num, Den, Exp); public: - static constexpr std::intmax_t num = std::get<0>(norm); - static constexpr std::intmax_t den = std::get<1>(norm); - static constexpr std::intmax_t exp = std::get<2>(norm); + static constexpr std::intmax_t num = norm[0]; + static constexpr std::intmax_t den = norm[1]; + static constexpr std::intmax_t exp = norm[2]; using type = ratio; + + static constexpr bool is_integral() { + if constexpr (exp < 0) { + return false; + } else { + return detail::gcdpow(num, exp, den) == den; + } + } }; namespace detail { @@ -79,49 +63,51 @@ namespace detail { template inline constexpr bool is_ratio> = true; -template -constexpr std::tuple ratio_add_detail() { - std::intmax_t num1 = R1::num; - std::intmax_t num2 = R2::num; +// unused, and align exponents process could be subject to overflow in extreme cases - // align exponents - std::intmax_t new_exp = R1::exp; - if constexpr (R1::exp > R2::exp) { - new_exp = R1::exp; - while (new_exp > R2::exp) { - num1 *= 10; - --new_exp; - } - } else if constexpr (R1::exp < R2::exp) { - new_exp = R2::exp; - while (R1::exp < new_exp) { - num2 *= 10; - --new_exp; - } - } +// template +// constexpr auto ratio_add_detail() { +// std::intmax_t num1 = R1::num; +// std::intmax_t num2 = R2::num; - // common denominator - std::intmax_t lcm_den = std::lcm(R1::den, R2::den); - num1 = num1 * (lcm_den / R1::den); - num2 = num2 * (lcm_den / R2::den); +// // align exponents +// std::intmax_t new_exp = R1::exp; +// if constexpr (R1::exp > R2::exp) { +// new_exp = R1::exp; +// while (new_exp > R2::exp) { +// num1 *= 10; +// --new_exp; +// } +// } else if constexpr (R1::exp < R2::exp) { +// new_exp = R2::exp; +// while (R1::exp < new_exp) { +// num2 *= 10; +// --new_exp; +// } +// } - return std::make_tuple(num1 + num2, lcm_den, new_exp); -} +// // common denominator +// std::intmax_t lcm_den = std::lcm(R1::den, R2::den); +// num1 = num1 * (lcm_den / R1::den); +// num2 = num2 * (lcm_den / R2::den); + +// return std::array{num1 + num2, lcm_den, new_exp}; +// } -template -struct ratio_add_impl { - static constexpr auto detail = ratio_add_detail(); - using type = ratio(detail), std::get<1>(detail), std::get<2>(detail)>; -}; +// template +// struct ratio_add_impl { +// static constexpr auto detail = ratio_add_detail(); +// using type = ratio; +// }; }// namespace detail -// ratio_add -template -using ratio_add = detail::ratio_add_impl::type; +// ratio_add : not used +// template +// using ratio_add = detail::ratio_add_impl::type; -// ratio_subtract +// ratio_subtract : not used // TODO implement ratio_subtract // template // using ratio_subtract = detail::ratio_subtract_impl::type; @@ -157,9 +143,9 @@ private: safe_multiply(R1::den / gcd2, R2::den / gcd1), R1::exp + R2::exp); - static constexpr std::intmax_t norm_num = std::get<0>(norm); - static constexpr std::intmax_t norm_den = std::get<1>(norm); - static constexpr std::intmax_t norm_exp = std::get<2>(norm); + static constexpr std::intmax_t norm_num = norm[0]; + static constexpr std::intmax_t norm_den = norm[1]; + static constexpr std::intmax_t norm_exp = norm[2]; public: using type = ratio; @@ -233,22 +219,22 @@ constexpr std::intmax_t sqrt_impl(std::intmax_t v, std::intmax_t l, std::intmax_ static constexpr std::intmax_t sqrt_impl(std::intmax_t v) { return sqrt_impl(v, 1, v); } template -constexpr std::tuple make_exp_even() +constexpr auto make_exp_even() { if constexpr (R::exp % 2 == 0) - return {R::num, R::den, R::exp}; // already even (incl zero) + return std::array{R::num, R::den, R::exp}; // already even (incl zero) // safely make exp even, so it can be divided by 2 for square root if constexpr (R::exp > 0) - return {R::num * 10, R::den, R::exp - 1}; + return std::array{R::num * 10, R::den, R::exp - 1}; else - return {R::num, R::den * 10, R::exp + 1}; + return std::array{R::num, R::den * 10, R::exp + 1}; } template struct ratio_sqrt_impl { constexpr static auto even = detail::make_exp_even(); - using type = ratio(even)), detail::sqrt_impl(std::get<1>(even)), std::get<2>(even) / 2>; + using type = ratio; }; template @@ -268,9 +254,11 @@ namespace detail { // TODO: simplified template struct common_ratio_impl { - static constexpr std::intmax_t gcd_num = std::gcd(R1::num, R2::num); - static constexpr std::intmax_t gcd_den = std::gcd(R1::den, R2::den); - using type = ratio; + static constexpr auto res = gcd_frac(R1::num, R1::den, R1::exp, R2::num, R2::den, R2::exp); + static constexpr std::intmax_t num = res[0]; + static constexpr std::intmax_t den = res[1]; + static constexpr std::intmax_t exp = res[2]; + using type = ratio; }; } // namespace detail diff --git a/test/unit_test/runtime/fmt_units_test.cpp b/test/unit_test/runtime/fmt_units_test.cpp index 22cd75a5..bf43b9f2 100644 --- a/test/unit_test/runtime/fmt_units_test.cpp +++ b/test/unit_test/runtime/fmt_units_test.cpp @@ -170,4 +170,9 @@ TEST_CASE("fmt::format on synthesized unit symbols", "[text][fmt]") { CHECK(fmt::format("{}", 1q_Npm) == "1 N/m"); } + + SECTION("addition with common ratio") + { + CHECK(fmt::format("{}", 1q_in + 1q_yd) == "37 in"); + } } diff --git a/test/unit_test/static/ratio_test.cpp b/test/unit_test/static/ratio_test.cpp index 987bd363..7e68a2fe 100644 --- a/test/unit_test/static/ratio_test.cpp +++ b/test/unit_test/static/ratio_test.cpp @@ -80,17 +80,16 @@ static_assert(std::is_same_v>, ratio<3, 1, 1>>); static_assert(std::is_same_v>, ratio<2>>); - static_assert(std::is_same_v, units::ratio<1, 3, 0>>, units::ratio<5, 6, 0>>); - static_assert(std::is_same_v, units::ratio<1, 3, 1>>, units::ratio<5, 6, 1>>); - static_assert(std::is_same_v, units::ratio<2, 7, 2>>, units::ratio<37, 56, 2>>); - - static_assert(std::is_same_v, units::ratio<2, 7, 1>>, units::ratio<226, 56, 1>>); - static_assert(std::is_same_v, units::ratio<3, 8, 2>>, units::ratio<226, 56, 1>>); - - static_assert(std::is_same_v, units::ratio<2, 7, -1>>, units::ratio<181, 56, -2>>); - static_assert(std::is_same_v, units::ratio<3, 8, -2>>, units::ratio<181, 56, -2>>); + // unused + // static_assert(std::is_same_v, units::ratio<1, 3, 0>>, units::ratio<5, 6, 0>>); + // static_assert(std::is_same_v, units::ratio<1, 3, 1>>, units::ratio<5, 6, 1>>); + // static_assert(std::is_same_v, units::ratio<2, 7, 2>>, units::ratio<37, 56, 2>>); + // static_assert(std::is_same_v, units::ratio<2, 7, 1>>, units::ratio<226, 56, 1>>); + // static_assert(std::is_same_v, units::ratio<3, 8, 2>>, units::ratio<226, 56, 1>>); + // static_assert(std::is_same_v, units::ratio<2, 7, -1>>, units::ratio<181, 56, -2>>); + // static_assert(std::is_same_v, units::ratio<3, 8, -2>>, units::ratio<181, 56, -2>>); // common_ratio