Merge pull request #329 from chiphogg/mag-value-t

Add value and categorization helpers for Magnitude

src/core/include/units/magnitude.h

@@ -25,6 +25,7 @@
 #include <units/ratio.h>
 #include <cstdint>
 #include <numbers>
+#include <stdexcept>
 
 namespace units {
 
@@ -117,6 +118,98 @@ static constexpr bool is_base_power<base_power<T>> = true;
 template<typename T>
 concept BasePower = detail::is_base_power<T>;
 
+namespace detail
+{
+constexpr auto inverse(BasePower auto bp) {
+  bp.power.num *= -1;
+  return bp;
+}
+
+// `widen_t` gives the widest arithmetic type in the same category, for intermediate computations.
+template<typename T> requires std::is_arithmetic_v<T>
+using widen_t = std::conditional_t<
+  std::is_floating_point_v<T>,
+  long double,
+  std::conditional_t<std::is_signed_v<T>, std::intmax_t, std::uintmax_t>>;
+
+// Raise an arbitrary arithmetic type to a positive integer power at compile time.
+template<typename T> requires std::is_arithmetic_v<T>
+constexpr T int_power(T base, std::integral auto exp){
+  // As this function should only be called at compile time, the exceptions herein function as
+  // "parameter-compatible static_asserts", and should not result in exceptions at runtime.
+  if (exp < 0) { throw std::invalid_argument{"int_power only supports positive integer powers"}; }
+
+  constexpr auto checked_multiply = [] (auto a, auto b) {
+    const auto result = a * b;
+    if (result / a != b) { throw std::overflow_error{"Wraparound detected"}; }
+    return result;
+  };
+
+  constexpr auto checked_square = [checked_multiply] (auto a) { return checked_multiply(a, a); };
+
+  // TODO(chogg): Unify this implementation with the one in pow.h.  That one takes its exponent as a
+  // template parameter, rather than a function parameter.
+
+  if (exp == 0) {
+    return 1;
+  }
+
+  if (exp % 2 == 1) {
+    return checked_multiply(base, int_power(base, exp - 1));
+  }
+
+  return checked_square(int_power(base, exp / 2));
+}
+
+
+template<typename T> requires std::is_arithmetic_v<T>
+constexpr widen_t<T> compute_base_power(BasePower auto bp)
+{
+  // This utility can only handle integer powers.  To compute rational powers at compile time, we'll
+  // need to write a custom function.
+  //
+  // Note that since this function should only be called at compile time, the point of these
+  // exceptions is to act as "static_assert substitutes", not to throw actual exceptions at runtime.
+  if (bp.power.den != 1) { throw std::invalid_argument{"Rational powers not yet supported"}; }
+  if (bp.power.exp < 0) { throw std::invalid_argument{"Unsupported exp value"}; }
+
+  if (bp.power.num < 0) {
+    if constexpr (std::is_integral_v<T>) {
+      throw std::invalid_argument{"Cannot represent reciprocal as integer"};
+    } else {
+      return 1 / compute_base_power<T>(inverse(bp));
+    }
+  }
+
+  auto power = bp.power.num * int_power(10, bp.power.exp);
+  return int_power(static_cast<widen_t<T>>(bp.get_base()), power);
+}
+
+// A converter for the value member variable of magnitude (below).
+//
+// The input is the desired result, but in a (wider) intermediate type.  The point of this function
+// is to cast to the desired type, but avoid overflow in doing so.
+template<typename To, typename From>
+  requires std::is_arithmetic_v<To>
+    && std::is_arithmetic_v<From>
+    && (std::is_integral_v<To> == std::is_integral_v<From>)
+constexpr To checked_static_cast(From x) {
+  // This function should only ever be called at compile time.  The purpose of these exceptions is
+  // to produce compiler errors, because we cannot `static_assert` on function arguments.
+  if constexpr (std::is_integral_v<To>) {
+    if (!std::in_range<To>(x)) {
+      throw std::invalid_argument{"Cannot represent magnitude in this type"};
+    }
+  } else {
+    if (x < std::numeric_limits<To>::min() || x > std::numeric_limits<To>::max()) {
+      throw std::invalid_argument{"Cannot represent magnitude in this type"};
+    }
+  }
+
+  return static_cast<To>(x);
+}
+} // namespace detail
+
 /**
  * @brief  Equality detection for two base powers.
  */
@@ -202,7 +295,7 @@ pairwise_all(T) -> pairwise_all<T>;
 
 // Check whether a sequence of (possibly heterogeneously typed) values are strictly increasing.
 template<typename... Ts>
-  requires ((std::is_signed_v<Ts> && ...))
+  requires (std::is_signed_v<Ts> && ...)
 constexpr bool strictly_increasing(Ts&&... ts) {
   return pairwise_all{std::less{}}(std::forward<Ts>(ts)...);
 }
@@ -215,6 +308,14 @@ static constexpr bool all_bases_in_order = strictly_increasing(BPs.get_base()...
 
 template<BasePower auto... BPs>
 static constexpr bool is_base_power_pack_valid = all_base_powers_valid<BPs...> && all_bases_in_order<BPs...>;
+
+constexpr bool is_rational(BasePower auto bp) {
+  return std::is_integral_v<decltype(bp.get_base())> && (bp.power.den == 1) && (bp.power.exp >= 0);
+}
+
+constexpr bool is_integral(BasePower auto bp) {
+  return is_rational(bp) && bp.power.num > 0;
+}
 } // namespace detail
 
 /**
@@ -224,8 +325,14 @@ static constexpr bool is_base_power_pack_valid = all_base_powers_valid<BPs...> &
  * rational powers, and compare for equality.
  */
 template<BasePower auto... BPs>
-  requires (detail::is_base_power_pack_valid<BPs...>)
-struct magnitude {};
+  requires detail::is_base_power_pack_valid<BPs...>
+struct magnitude {
+  // Whether this magnitude represents an integer.
+  friend constexpr bool is_integral(const magnitude&) { return (detail::is_integral(BPs) && ...); }
+
+  // Whether this magnitude represents a rational number.
+  friend constexpr bool is_rational(const magnitude&) { return (detail::is_rational(BPs) && ...); }
+};
 
 // Implementation for Magnitude concept (below).
 namespace detail {
@@ -241,6 +348,18 @@ static constexpr bool is_magnitude<magnitude<BPs...>> = true;
 template<typename T>
 concept Magnitude = detail::is_magnitude<T>;
 
+/**
+ * @brief  The value of a Magnitude in a desired type T.
+ */
+template<typename T, BasePower auto... BPs>
+  requires std::is_floating_point_v<T> || (std::is_integral_v<T> && is_integral(magnitude<BPs...>{}))
+constexpr T get_value(const magnitude<BPs...> &) {
+  // Force the expression to be evaluated in a constexpr context, to catch, e.g., overflow.
+  constexpr auto result = detail::checked_static_cast<T>((detail::compute_base_power<T>(BPs) * ...));
+
+  return result;
+}
+
 /**
  * @brief  A base to represent pi.
  */

test/unit_test/runtime/magnitude_test.cpp

@@ -37,6 +37,12 @@ struct invalid_negative_base { static constexpr long double value = -1.234L; };
 template<ratio Power>
 constexpr auto pi_to_the() { return magnitude<base_power<pi_base>{Power}>{}; }
 
+template<typename T, typename U>
+void check_same_type_and_value(T actual, U expected) {
+  CHECK(std::is_same_v<T, U>);
+  CHECK(actual == expected);
+}
+
 TEST_CASE("base_power")
 {
   SECTION("base rep deducible for integral base")
@@ -129,6 +135,71 @@ TEST_CASE("make_ratio performs prime factorization correctly")
   }
 }
 
+TEST_CASE("magnitude converts to numerical value")
+{
+  SECTION("Positive integer powers of integer bases give integer values")
+  {
+    constexpr auto mag_412 = as_magnitude<412>();
+    check_same_type_and_value(get_value<int>(mag_412), 412);
+    check_same_type_and_value(get_value<std::size_t>(mag_412), std::size_t{412});
+    check_same_type_and_value(get_value<float>(mag_412), 412.0f);
+    check_same_type_and_value(get_value<double>(mag_412), 412.0);
+  }
+
+  SECTION("Negative integer powers of integer bases compute correct values")
+  {
+    constexpr auto mag_0p125 = as_magnitude<ratio{1, 8}>();
+    check_same_type_and_value(get_value<float>(mag_0p125), 0.125f);
+    check_same_type_and_value(get_value<double>(mag_0p125), 0.125);
+  }
+
+  SECTION("pi to the 1 supplies correct values")
+  {
+    constexpr auto pi = pi_to_the<1>();
+    check_same_type_and_value(get_value<float>(pi), std::numbers::pi_v<float>);
+    check_same_type_and_value(get_value<double>(pi), std::numbers::pi_v<double>);
+    check_same_type_and_value(get_value<long double>(pi), std::numbers::pi_v<long double>);
+  }
+
+  SECTION("pi to arbitrary power performs computations in most accurate type at compile time")
+  {
+    if constexpr (sizeof(float) < sizeof(long double)) {
+      constexpr auto pi_cubed = pi_to_the<3>();
+
+      auto cube = [](auto x) { return x * x * x; };
+      constexpr auto via_float = cube(std::numbers::pi_v<float>);
+      constexpr auto via_long_double = static_cast<float>(cube(std::numbers::pi_v<long double>));
+
+      constexpr auto pi_cubed_value = get_value<float>(pi_cubed);
+      REQUIRE(pi_cubed_value != via_float);
+      CHECK(pi_cubed_value == via_long_double);
+    }
+  }
+
+  SECTION("Impossible requests are prevented at compile time")
+  {
+    // Naturally, we cannot actually write a test to verify a compiler error.  But any of these can
+    // be uncommented if desired to verify that it breaks the build.
+
+    // get_value<int8_t>(as_magnitude<412>());
+
+    // Would work for pow<62>:
+    // get_value<int64_t>(pow<63>(as_magnitude<2>()));
+
+    // Would work for pow<63>:
+    // get_value<uint64_t>(pow<64>(as_magnitude<2>()));
+
+    get_value<double>(pow<308>(as_magnitude<10>())); // Compiles, correctly.
+    // get_value<double>(pow<309>(as_magnitude<10>()));
+    // get_value<double>(pow<3099>(as_magnitude<10>()));
+    // get_value<double>(pow<3099999>(as_magnitude<10>()));
+
+    auto sqrt_2 = pow<ratio{1, 2}>(as_magnitude<2>());
+    CHECK(!is_integral(sqrt_2));
+    // get_value<int>(sqrt_2);
+  }
+}
+
 TEST_CASE("Equality works for magnitudes")
 {
   SECTION("Equivalent ratios are equal")
@@ -218,8 +289,52 @@ TEST_CASE("Can raise Magnitudes to rational powers")
   }
 }
 
+TEST_CASE("can distinguish integral, rational, and irrational magnitudes")
+{
+  SECTION("Integer magnitudes are integral and rational")
+  {
+    auto check_rational_and_integral = [](Magnitude auto m) {
+      CHECK(is_integral(m));
+      CHECK(is_rational(m));
+    };
+    check_rational_and_integral(magnitude<>{});
+    check_rational_and_integral(as_magnitude<1>());
+    check_rational_and_integral(as_magnitude<3>());
+    check_rational_and_integral(as_magnitude<8>());
+    check_rational_and_integral(as_magnitude<412>());
+    check_rational_and_integral(as_magnitude<ratio{1, 1}>());
+  }
+
+  SECTION("Fractional magnitudes are rational, but not integral")
+  {
+    auto check_rational_but_not_integral = [](Magnitude auto m) {
+      CHECK(!is_integral(m));
+      CHECK(is_rational(m));
+    };
+    check_rational_but_not_integral(as_magnitude<ratio{1, 2}>());
+    check_rational_but_not_integral(as_magnitude<ratio{5, 8}>());
+  }
+}
+
 namespace detail {
 
+TEST_CASE("int_power computes integer powers") {
+  SECTION("handles floating point") {
+    check_same_type_and_value(int_power(0.123L, 0), 1.0L);
+    check_same_type_and_value(int_power(0.246f, 1), 0.246f);
+    check_same_type_and_value(int_power(0.5f, 3), 0.125f);
+    check_same_type_and_value(int_power(2.5, 4), 39.0625);
+
+    CHECK(std::is_same_v<long double, decltype(compute_base_power<double>(base_power{10, 20}))>);
+  }
+
+  SECTION("handles integral") {
+    check_same_type_and_value(int_power(8, 0), 1);
+    check_same_type_and_value(int_power(9L, 1), 9L);
+    check_same_type_and_value(int_power(2, 10), 1024);
+  }
+}
+
 TEST_CASE("Prime helper functions")
 {
   SECTION("find_first_factor()") {