Support computing values for all rational powers (#213)

chiphogg · web-flow · commit 7d8a5447b8b7 · 2023-12-19T11:48:40.000-05:00
I thought we would need to wait for a `constexpr` version of `std::powl`, but I was wrong. Instead, we can get away with a `root` function. Note that there's no guarantee that `std::powl` would be better in all cases, even if we _could_ just use it. After all, its second argument is `1.0l / n`, which is a lossy representation of the Nth root. On the implementation side, we take advantage of the fact that this function is only ever meant to be called at compile time, which lets us prioritize accuracy over "speed" (because the runtime speed is always essentially infinite). We do a binary search between 1 and `x` after checking that `x > 1` (and `n > 1`). When we end up with two neighboring floating point values, we pick whichever gives the closest result when we run it through `checked_int_pow`. This algorithm is guaranteed to produce the best representable result for a "pure root". For other rational powers, it's technically possible that we could have some lossiness in the `checked_int_pow` computation, if we enter a floating point realm that is high enough such that not all integer values can be represented. That said, I tried a bunch of test cases to compare it to `powl`, and it's passed all the ones that I was able to come up with. Fixes #116.
diff --git a/au/magnitude.hh b/au/magnitude.hh
@@ -270,7 +270,7 @@ namespace detail {
 enum class MagRepresentationOutcome {
     OK,
     ERR_NON_INTEGER_IN_INTEGER_TYPE,
-    ERR_RATIONAL_POWERS,
+    ERR_INVALID_ROOT,
     ERR_CANNOT_FIT,
 };
 
@@ -316,10 +316,101 @@ constexpr MagRepresentationOrError<T> checked_int_pow(T base, std::uintmax_t exp
     return result;
 }
 
-template <typename T, std::intmax_t N, typename B>
+template <typename T>
+constexpr MagRepresentationOrError<T> root(T x, std::uintmax_t n) {
+    // The "zeroth root" would be mathematically undefined.
+    if (n == 0) {
+        return {MagRepresentationOutcome::ERR_INVALID_ROOT};
+    }
+
+    // The "first root" is trivial.
+    if (n == 1) {
+        return {MagRepresentationOutcome::OK, x};
+    }
+
+    // We only support nontrivial roots of floating point types.
+    if (!std::is_floating_point<T>::value) {
+        return {MagRepresentationOutcome::ERR_NON_INTEGER_IN_INTEGER_TYPE};
+    }
+
+    // Handle negative numbers: only odd roots are allowed.
+    if (x < 0) {
+        if (n % 2 == 0) {
+            return {MagRepresentationOutcome::ERR_INVALID_ROOT};
+        } else {
+            const auto negative_result = root(-x, n);
+            if (negative_result.outcome != MagRepresentationOutcome::OK) {
+                return {negative_result.outcome};
+            }
+            return {MagRepresentationOutcome::OK, static_cast<T>(-negative_result.value)};
+        }
+    }
+
+    // Handle special cases of zero and one.
+    if (x == 0 || x == 1) {
+        return {MagRepresentationOutcome::OK, x};
+    }
+
+    // Handle numbers bewtween 0 and 1.
+    if (x < 1) {
+        const auto inverse_result = root(T{1} / x, n);
+        if (inverse_result.outcome != MagRepresentationOutcome::OK) {
+            return {inverse_result.outcome};
+        }
+        return {MagRepresentationOutcome::OK, static_cast<T>(T{1} / inverse_result.value)};
+    }
+
+    //
+    // At this point, error conditions are finished, and we can proceed with the "core" algorithm.
+    //
+
+    // Always use `long double` for intermediate computations.  We don't ever expect people to be
+    // calling this at runtime, so we want maximum accuracy.
+    long double lo = 1.0;
+    long double hi = x;
+
+    // Do a binary search to find the closest value such that `checked_int_pow` recovers the input.
+    //
+    // Because we know `n > 1`, and `x > 1`, and x^n is monotonically increasing, we know that
+    // `checked_int_pow(lo, n) < x < checked_int_pow(hi, n)`.  We will preserve this as an
+    // invariant.
+    while (lo < hi) {
+        long double mid = lo + (hi - lo) / 2;
+
+        auto result = checked_int_pow(mid, n);
+
+        if (result.outcome != MagRepresentationOutcome::OK) {
+            return {result.outcome};
+        }
+
+        // Early return if we get lucky with an exact answer.
+        if (result.value == x) {
+            return {MagRepresentationOutcome::OK, static_cast<T>(mid)};
+        }
+
+        // Check for stagnation.
+        if (mid == lo || mid == hi) {
+            break;
+        }
+
+        // Preserve the invariant that `checked_int_pow(lo, n) < x < checked_int_pow(hi, n)`.
+        if (result.value < x) {
+            lo = mid;
+        } else {
+            hi = mid;
+        }
+    }
+
+    // Pick whichever one gets closer to the target.
+    const auto lo_diff = x - checked_int_pow(lo, n).value;
+    const auto hi_diff = checked_int_pow(hi, n).value - x;
+    return {MagRepresentationOutcome::OK, static_cast<T>(lo_diff < hi_diff ? lo : hi)};
+}
+
+template <typename T, std::intmax_t N, std::uintmax_t D, typename B>
 constexpr MagRepresentationOrError<Widen<T>> base_power_value(B base) {
     if (N < 0) {
-        const auto inverse_result = base_power_value<T, -N>(base);
+        const auto inverse_result = base_power_value<T, -N, D>(base);
         if (inverse_result.outcome != MagRepresentationOutcome::OK) {
             return inverse_result;
         }
@@ -329,7 +420,12 @@ constexpr MagRepresentationOrError<Widen<T>> base_power_value(B base) {
         };
     }
 
-    return checked_int_pow(static_cast<Widen<T>>(base), static_cast<std::uintmax_t>(N));
+    const auto power_result =
+        checked_int_pow(static_cast<Widen<T>>(base), static_cast<std::uintmax_t>(N));
+    if (power_result.outcome != MagRepresentationOutcome::OK) {
+        return {power_result.outcome};
+    }
+    return (D > 1) ? root(power_result.value, D) : power_result;
 }
 
 template <typename T, std::size_t N>
@@ -393,15 +489,10 @@ constexpr MagRepresentationOrError<T> get_value_result(Magnitude<BPs...>) {
         return {MagRepresentationOutcome::ERR_NON_INTEGER_IN_INTEGER_TYPE};
     }
 
-    // Computing values for rational base powers is something we would _like_ to support, but we
-    // need a `constexpr` implementation of `powl()` first.
-    if (!all({(ExpT<BPs>::den == 1)...})) {
-        return {MagRepresentationOutcome::ERR_RATIONAL_POWERS};
-    }
-
     // Force the expression to be evaluated in a constexpr context.
     constexpr auto widened_result =
-        product({base_power_value<T, (ExpT<BPs>::num / ExpT<BPs>::den)>(BaseT<BPs>::value())...});
+        product({base_power_value<T, ExpT<BPs>::num, static_cast<std::uintmax_t>(ExpT<BPs>::den)>(
+            BaseT<BPs>::value())...});
 
     if ((widened_result.outcome != MagRepresentationOutcome::OK) ||
         !safe_to_cast_to<T>(widened_result.value)) {
@@ -433,8 +524,8 @@ constexpr T get_value(Magnitude<BPs...> m) {
 
     static_assert(result.outcome != MagRepresentationOutcome::ERR_NON_INTEGER_IN_INTEGER_TYPE,
                   "Cannot represent non-integer in integral destination type");
-    static_assert(result.outcome != MagRepresentationOutcome::ERR_RATIONAL_POWERS,
-                  "Computing values for rational powers not yet supported");
+    static_assert(result.outcome != MagRepresentationOutcome::ERR_INVALID_ROOT,
+                  "Could not compute root for rational power of base");
     static_assert(result.outcome != MagRepresentationOutcome::ERR_CANNOT_FIT,
                   "Value outside range of destination type");
 
diff --git a/au/magnitude_test.cc b/au/magnitude_test.cc
@@ -20,6 +20,8 @@
 #include "gtest/gtest.h"
 
 using ::testing::DoubleEq;
+using ::testing::Eq;
+using ::testing::FloatEq;
 using ::testing::StaticAssertTypeEq;
 
 namespace au {
@@ -208,6 +210,11 @@ TEST(GetValue, ImpossibleRequestsArePreventedAtCompileTime) {
     // get_value<int>(sqrt_2);
 }
 
+TEST(GetValue, HandlesRoots) {
+    constexpr auto sqrt_2 = get_value<double>(root<2>(mag<2>()));
+    EXPECT_DOUBLE_EQ(sqrt_2 * sqrt_2, 2.0);
+}
+
 TEST(GetValue, WorksForEmptyPack) {
     constexpr auto one = Magnitude<>{};
     EXPECT_THAT(get_value<int>(one), SameTypeAndValue(1));
@@ -270,6 +277,15 @@ MATCHER(CannotFit, "") {
     return (arg.outcome == MagRepresentationOutcome::ERR_CANNOT_FIT) && (arg.value == 0);
 }
 
+MATCHER(NonIntegerInIntegerType, "") {
+    return (arg.outcome == MagRepresentationOutcome::ERR_NON_INTEGER_IN_INTEGER_TYPE) &&
+           (arg.value == 0);
+}
+
+MATCHER(InvalidRoot, "") {
+    return (arg.outcome == MagRepresentationOutcome::ERR_INVALID_ROOT) && (arg.value == 0);
+}
+
 template <typename T, typename ValueMatcher>
 auto FitsAndMatchesValue(ValueMatcher &&matcher) {
     return ::testing::AllOf(
@@ -298,6 +314,106 @@ TEST(CheckedIntPow, FindsAppropriateLimits) {
     EXPECT_THAT(checked_int_pow(10.0, 309), CannotFit());
 }
 
+TEST(Root, ReturnsErrorForIntegralType) {
+    EXPECT_THAT(root(4, 2), NonIntegerInIntegerType());
+    EXPECT_THAT(root(uint8_t{125}, 3), NonIntegerInIntegerType());
+}
+
+TEST(Root, ReturnsErrorForZerothRoot) {
+    EXPECT_THAT(root(4.0, 0), InvalidRoot());
+    EXPECT_THAT(root(125.0, 0), InvalidRoot());
+}
+
+TEST(Root, NegativeRootsWorkForOddPowersOnly) {
+    EXPECT_THAT(root(-4.0, 2), InvalidRoot());
+    EXPECT_THAT(root(-125.0, 3), FitsAndProducesValue(-5.0));
+    EXPECT_THAT(root(-10000.0, 4), InvalidRoot());
+}
+
+TEST(Root, AnyRootOfOneIsOne) {
+    for (const std::uintmax_t r : {1u, 2u, 3u, 4u, 5u, 6u, 7u, 8u, 9u}) {
+        EXPECT_THAT(root(1.0, r), FitsAndProducesValue(1.0));
+    }
+}
+
+TEST(Root, AnyRootOfZeroIsZero) {
+    for (const std::uintmax_t r : {1u, 2u, 3u, 4u, 5u, 6u, 7u, 8u, 9u}) {
+        EXPECT_THAT(root(0.0, r), FitsAndProducesValue(0.0));
+    }
+}
+
+TEST(Root, OddRootOfNegativeOneIsItself) {
+    EXPECT_THAT(root(-1.0, 1), FitsAndProducesValue(-1.0));
+    EXPECT_THAT(root(-1.0, 2), InvalidRoot());
+    EXPECT_THAT(root(-1.0, 3), FitsAndProducesValue(-1.0));
+    EXPECT_THAT(root(-1.0, 4), InvalidRoot());
+    EXPECT_THAT(root(-1.0, 5), FitsAndProducesValue(-1.0));
+}
+
+TEST(Root, RecoversExactValueWherePossible) {
+    {
+        const auto sqrt_4f = root(4.0f, 2);
+        EXPECT_THAT(sqrt_4f.outcome, Eq(MagRepresentationOutcome::OK));
+        EXPECT_THAT(sqrt_4f.value, SameTypeAndValue(2.0f));
+    }
+
+    {
+        const auto cbrt_125L = root(125.0L, 3);
+        EXPECT_THAT(cbrt_125L.outcome, Eq(MagRepresentationOutcome::OK));
+        EXPECT_THAT(cbrt_125L.value, SameTypeAndValue(5.0L));
+    }
+}
+
+TEST(Root, HandlesArgumentsBetweenOneAndZero) {
+    EXPECT_THAT(root(0.25, 2), FitsAndProducesValue(0.5));
+    EXPECT_THAT(root(0.0001, 4), FitsAndMatchesValue<double>(DoubleEq(0.1)));
+}
+
+TEST(Root, ResultIsVeryCloseToStdPowForPureRoots) {
+    for (const double x : {55.5, 123.456, 789.012, 3456.789, 12345.6789, 5.67e25}) {
+        for (const auto r : {2u, 3u, 4u, 5u, 6u, 7u, 8u, 9u}) {
+            const auto double_result = root(x, r);
+            EXPECT_THAT(double_result.outcome, Eq(MagRepresentationOutcome::OK));
+            EXPECT_THAT(double_result.value, DoubleEq(static_cast<double>(std::pow(x, 1.0L / r))));
+
+            const auto float_result = root(static_cast<float>(x), r);
+            EXPECT_THAT(float_result.outcome, Eq(MagRepresentationOutcome::OK));
+            EXPECT_THAT(float_result.value, FloatEq(static_cast<float>(std::pow(x, 1.0L / r))));
+        }
+    }
+}
+
+TEST(Root, ResultAtLeastAsGoodAsStdPowForRationalPowers) {
+    struct RationalPower {
+        std::uintmax_t num;
+        std::uintmax_t den;
+    };
+
+    auto result_via_root = [](double x, RationalPower power) {
+        return static_cast<double>(
+            root(checked_int_pow(static_cast<long double>(x), power.num).value, power.den).value);
+    };
+
+    auto result_via_std_pow = [](double x, RationalPower power) {
+        return static_cast<double>(
+            std::pow(static_cast<long double>(x),
+                     static_cast<long double>(power.num) / static_cast<long double>(power.den)));
+    };
+
+    auto round_trip_error = [](double x, RationalPower power, auto func) {
+        const auto round_trip_result = func(func(x, power), {power.den, power.num});
+        return std::abs(round_trip_result - x);
+    };
+
+    for (const auto base : {2.0, 3.1415, 98.6, 1.2e-10, 5.5e15}) {
+        for (const auto power : std::vector<RationalPower>{{5, 2}, {2, 3}, {7, 4}}) {
+            const auto error_from_root = round_trip_error(base, power, result_via_root);
+            const auto error_from_std_pow = round_trip_error(base, power, result_via_std_pow);
+            EXPECT_LE(error_from_root, error_from_std_pow);
+        }
+    }
+}
+
 TEST(GetValueResult, HandlesNumbersTooBigForUintmax) {
     EXPECT_THAT(get_value_result<std::uintmax_t>(pow<64>(mag<2>())), CannotFit());
 }
diff --git a/au/quantity_test.cc b/au/quantity_test.cc
@@ -263,18 +263,20 @@ TEST(Quantity, HandlesBaseDimensionsWithFractionalExponents) {
 }
 
 TEST(Quantity, HandlesMagnitudesWithFractionalExponents) {
-    using RootKiloFeet = decltype(root<2>(Kilo<Feet>{}));
-    constexpr auto x = make_quantity<RootKiloFeet>(3);
+    constexpr auto x = sqrt(kilo(feet))(3.0);
 
     // We can retrieve the value in the same unit (regardless of the scale's fractional powers).
-    EXPECT_EQ(x.in(RootKiloFeet{}), 3);
+    EXPECT_EQ(x.in(sqrt(kilo(feet))), 3.0);
 
     // We can retrieve the value in a *different* unit, which *also* has fractional powers, as long
     // as their *ratio* has no fractional powers.
-    EXPECT_EQ(x.in(root<2>(Milli<Feet>{})), 3'000);
+    EXPECT_EQ(x.in(sqrt(milli(feet))), 3'000.0);
+
+    // We can also retrieve the value in a different unit whose ratio *does* have fractional powers.
+    EXPECT_NEAR(x.in(sqrt(feet)), 94.86833, 1e-5);
 
     // Squaring the fractional base power gives us an exact non-fractional dimension and scale.
-    EXPECT_EQ(x * x, kilo(feet)(9));
+    EXPECT_EQ(x * x, kilo(feet)(9.0));
 }
 
 // A custom "Quantity-equivalent" type, whose interop with Quantity we'll provide below.
