Compute values for rational magnitude powers

Since this will only ever be done at compile time (as guaranteed by
using `consteval`), we can afford to prioritize precision over speed.
To compute an Nth root, we simply do a binary search over representable
floating point numbers, looking for the number whose Nth power most
closely matches the original number.

Fixes #494.  We have included a test case reproducing the original
problem exactly.  All tests use "within 4 ULPs" as the criterion, which
is (I believe) equivalent to the googletest `EXPECT_DOUBLE_EQ`
criterion.

src/core/include/mp-units/framework/magnitude.h

@@ -307,6 +307,119 @@ template<typename T>
   return checked_square(int_power(base, exp / 2));
 }
 
+template <typename T>
+[[nodiscard]] consteval std::optional<T> checked_int_pow(T base, std::uintmax_t exp) {
+    T result = T{1};
+    while (exp > 0u) {
+        if (exp % 2u == 1u) {
+            if (base > std::numeric_limits<T>::max() / result) {
+                return std::nullopt;
+            }
+            result *= base;
+        }
+
+        exp /= 2u;
+
+        if (base > std::numeric_limits<T>::max() / base) {
+            return (exp == 0u)
+                       ? std::make_optional(result)
+                       : std::nullopt;
+        }
+        base *= base;
+    }
+    return result;
+}
+
+template <typename T>
+[[nodiscard]] consteval std::optional<T> root(T x, std::uintmax_t n) {
+    // The "zeroth root" would be mathematically undefined.
+    if (n == 0) {
+        return std::nullopt;
+    }
+
+    // The "first root" is trivial.
+    if (n == 1) {
+        return x;
+    }
+
+    // We only support nontrivial roots of floating point types.
+    if (!std::is_floating_point<T>::value) {
+        return std::nullopt;
+    }
+
+    // Handle negative numbers: only odd roots are allowed.
+    if (x < 0) {
+        if (n % 2 == 0) {
+            return std::nullopt;
+        } else {
+            const auto negative_result = root(-x, n);
+            if (!negative_result.has_value()) {
+                return std::nullopt;
+            }
+            return static_cast<T>(-negative_result.value());
+        }
+    }
+
+    // Handle special cases of zero and one.
+    if (x == 0 || x == 1) {
+        return x;
+    }
+
+    // Handle numbers bewtween 0 and 1.
+    if (x < 1) {
+        const auto inverse_result = root(T{1} / x, n);
+        if (!inverse_result.has_value()) {
+            return std::nullopt;
+        }
+        return 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 = static_cast<long double>(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.has_value()) {
+            return std::nullopt;
+        }
+
+        // Early return if we get lucky with an exact answer.
+        if (result.value() == x) {
+            return 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 static_cast<T>(lo_diff < hi_diff ? lo : hi);
+}
 
 template<typename T>
 [[nodiscard]] consteval widen_t<T> compute_base_power(MagnitudeSpec auto el)
@@ -317,9 +430,6 @@ template<typename T>
   // Note that since this function should only be called at compile time, the point of these
   // terminations is to act as "static_assert substitutes", not to actually terminate at runtime.
   const auto exp = get_exponent(el);
-  if (exp.den != 1) {
-    std::abort();  // Rational powers not yet supported
-  }
 
   if (exp.num < 0) {
     if constexpr (std::is_integral_v<T>) {
@@ -329,8 +439,19 @@ template<typename T>
     }
   }
 
-  auto power = exp.num;
-  return int_power(static_cast<widen_t<T>>(get_base_value(el)), power);
+  const auto pow_result = checked_int_pow(static_cast<widen_t<T>>(get_base_value(el)), static_cast<uintmax_t>(exp.num));
+  if (pow_result.has_value()) {
+    const auto final_result = (exp.den > 1) ? root(pow_result.value(), static_cast<uintmax_t>(exp.den)) : pow_result;
+    if (final_result.has_value()) {
+      return final_result.value();
+    }
+    else {
+      std::abort(); // Root computation failed.
+    }
+  }
+  else {
+    std::abort(); // Power computation failed.
+  }
 }
 
 // A converter for the value member variable of magnitude (below).

test/static/quantity_test.cpp

@@ -24,6 +24,7 @@
 #include <mp-units/bits/hacks.h>
 #include <mp-units/ext/fixed_string.h>
 #include <mp-units/ext/type_traits.h>
+#include <mp-units/math.h>
 #include <mp-units/systems/isq/mechanics.h>
 #include <mp-units/systems/isq/space_and_time.h>
 #include <mp-units/systems/si.h>
@@ -52,6 +53,24 @@ using namespace mp_units::si::unit_symbols;
 // quantity class invariants
 //////////////////////////////
 
+template <typename T>
+constexpr bool within_4_ulps(T a, T b) {
+  static_assert(std::is_floating_point_v<T>);
+  auto walk_ulps = [](T x, int n) {
+    while (n > 0) {
+      x = std::nextafter(x, std::numeric_limits<T>::infinity());
+      --n;
+    }
+    while (n < 0) {
+      x = std::nextafter(x, -std::numeric_limits<T>::infinity());
+      ++n;
+    }
+    return x;
+  };
+
+  return (walk_ulps(a, -4) <= b) && (b <= walk_ulps(a, 4));
+}
+
 static_assert(sizeof(quantity<si::metre>) == sizeof(double));
 static_assert(sizeof(quantity<isq::length[m]>) == sizeof(double));
 static_assert(sizeof(quantity<si::metre, short>) == sizeof(short));
@@ -199,6 +218,16 @@ static_assert(std::convertible_to<quantity<isq::length[km], int>, quantity<isq::
 static_assert(std::constructible_from<quantity<isq::length[km]>, quantity<isq::length[m], int>>);
 static_assert(std::convertible_to<quantity<isq::length[m], int>, quantity<isq::length[km]>>);
 
+// conversion requiring radical magnitudes
+static_assert(within_4_ulps(sqrt((1.0 * m) * (1.0 * km)).numerical_value_in(m), sqrt(1000.0)));
+
+// Reproducing issue #494 exactly:
+constexpr auto val_issue_494 = 8.0 * si::si2019::boltzmann_constant * 1000.0 * K / (std::numbers::pi * 10 * Da);
+static_assert(
+  within_4_ulps(
+    sqrt(val_issue_494).numerical_value_in(m / s),
+    sqrt(val_issue_494.numerical_value_in(m * m / s / s))));
+
 
 ///////////////////////
 // obtaining a number