employing more mathemtaically correct ratio_gcd calc (for common ratio)

really finds the maximum common ration as opposed to previous algo which simplified on the exp part of the ratio by using std::min most of new code credit to Conor Williams discussion and additional doc here: https://github.com/mpusz/units/issues/62#issuecomment-588152833 test case was 1yd + 1in = 37in => added as a test commenting out unusued ratio_add and its tests if to be reintroduced, should also use the new gcd routines additonal change was required to check in `safe_divisible` concept den=1 is not sufficient anymore. reusing new gcd routines moved ratio nomalize and new gcd routines into new, separate bits/ratio_maths.h this resolves #62
2020-02-20 19:59:12 +01:00
parent 1280b7d4be
commit 39a2c2de0e
5 changed files with 248 additions and 83 deletions
@@ -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 <units/bits/external/hacks.h>
+#include <units/concepts.h>
+#include <algorithm>
+#include <cassert>
+#include <cmath>
+#include <cstdint>
+#include <numeric>
+#include <tuple>
+#include <type_traits>
+
+namespace units::detail {
+
+template<typename T>
+[[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<std::intmax_t>(mulmod(static_cast<std::uint64_t>(a % b),
+                                           modpow(10, static_cast<std::uint64_t>(e), static_cast<std::uint64_t>(b)),
+                                           static_cast<std::uint64_t>(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<std::intmax_t>::max() / n1 > d2);
+  assert(std::numeric_limits<std::intmax_t>::max() / n2 > d1);
+
+  std::intmax_t num = detail::gcdpow(n1 * d2, e1 - e2, n2 * d1);
+
+  assert(std::numeric_limits<std::intmax_t>::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