From a6b6afe438030ffe02bfe326e6c2942020f1a755 Mon Sep 17 00:00:00 2001 From: Chip Hogg Date: Fri, 28 Jan 2022 19:06:26 +0000 Subject: [PATCH] Add value and categorization helpers for Magnitude For each Magnitude `m`, we now support: - `is_rational(m)` tells whether it represents a rational number. - `is_integral(m)` tells whether it represents an integer (and `is_integral(m)` implies `is_rational(m)`). - `m.value` represents the value of `m` in the type `T`. If `T` is integral, we only support `m.value` when `is_integral(m)`. We also perform all intermediate computations in the widest type of the same category (floating point, signed, or unsigned). This means we can, for example, give first-class support to embedded users who may have hardware support only for `float`: we can ensure they get the most accurate `float` value we can compute, without burdening them with a dependency on `long double` in their runtime code. We are not yet ready to replace `ratio` as the definition of unit magnitudes, but we're a step closer. The next step will be to decompose a Magnitude into numerator, denominator, and irrational parts, giving us full bidirectional convertibility with existing `ratio` instances. Then we can migrate over in a controlled fashion. --- src/core/include/units/magnitude.h | 115 +++++++++++++++++++++- test/unit_test/runtime/magnitude_test.cpp | 111 +++++++++++++++++++++ 2 files changed, 225 insertions(+), 1 deletion(-) diff --git a/src/core/include/units/magnitude.h b/src/core/include/units/magnitude.h index 8e9c46ca..c33a61d3 100644 --- a/src/core/include/units/magnitude.h +++ b/src/core/include/units/magnitude.h @@ -25,6 +25,7 @@ #include #include #include +#include namespace units { @@ -117,6 +118,98 @@ static constexpr bool is_base_power> = true; template concept BasePower = detail::is_base_power; +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 requires std::is_arithmetic_v +using widen_t = std::conditional_t< + std::is_floating_point_v, + long double, + std::conditional_t, std::intmax_t, std::uintmax_t>>; + +// Raise an arbitrary arithmetic type to a positive integer power at compile time. +template requires std::is_arithmetic_v +consteval T int_power(T base, std::integral auto exp){ + if (exp < 0) { throw std::invalid_argument{"int_power only supports positive integer powers"}; } + + // 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 base * int_power(base, exp - 1); + } + + const auto square_root = int_power(base, exp / 2); + const auto result = square_root * square_root; + + if constexpr(std::is_unsigned_v) { + // As this function can only be called at compile time, the exception functions as a + // "parameter-compatible static_assert", and does not result in exceptions at runtime. + if (result / square_root != square_root) { throw std::overflow_error{"Unsigned wraparound"}; } + } + + return result; +} + + +template requires std::is_arithmetic_v +consteval widen_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 can 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) { + throw std::invalid_argument{"Cannot represent reciprocal as integer"}; + } else { + return 1 / compute_base_power(inverse(bp)); + } + } + + auto power = bp.power.num * int_power(10, bp.power.exp); + return int_power(static_cast>(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 + requires std::is_arithmetic_v + && std::is_arithmetic_v + && (std::is_integral_v == std::is_integral_v) +consteval To checked_static_cast(From x) { + // This function can 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) { + if (std::cmp_less(x, std::numeric_limits::min()) || + std::cmp_greater(x, std::numeric_limits::max())) { + throw std::invalid_argument{"Cannot represent magnitude in this type"}; + } + } else { + if (x < std::numeric_limits::min() || x > std::numeric_limits::max()) { + throw std::invalid_argument{"Cannot represent magnitude in this type"}; + } + } + + return static_cast(x); +} +} // namespace detail + /** * @brief Equality detection for two base powers. */ @@ -215,6 +308,14 @@ static constexpr bool all_bases_in_order = strictly_increasing(BPs.get_base()... template static constexpr bool is_base_power_pack_valid = all_base_powers_valid && all_bases_in_order; + +constexpr bool is_rational(BasePower auto bp) { + return std::is_integral_v && (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 /** @@ -225,7 +326,19 @@ static constexpr bool is_base_power_pack_valid = all_base_powers_valid & */ template requires (detail::is_base_power_pack_valid) -struct magnitude {}; +struct magnitude { + // Whether this magnitude represents an integer. + friend constexpr bool is_integral(magnitude) { return (detail::is_integral(BPs) && ...); } + + // Whether this magnitude represents a rational number. + friend constexpr bool is_rational(magnitude) { return (detail::is_rational(BPs) && ...); } + + // The value of this magnitude, expressed in a given type. + template + requires (is_integral(magnitude{}) || std::is_floating_point_v) + static constexpr T value = detail::checked_static_cast( + (detail::compute_base_power(BPs) * ...)); +}; // Implementation for Magnitude concept (below). namespace detail { diff --git a/test/unit_test/runtime/magnitude_test.cpp b/test/unit_test/runtime/magnitude_test.cpp index 653a6ef8..f41e1955 100644 --- a/test/unit_test/runtime/magnitude_test.cpp +++ b/test/unit_test/runtime/magnitude_test.cpp @@ -37,6 +37,12 @@ struct invalid_negative_base { static constexpr long double value = -1.234L; }; template constexpr auto pi_to_the() { return magnitude{Power}>{}; } +template +void check_same_type_and_value(T actual, U expected) { + CHECK(std::is_same_v); + CHECK(actual == expected); +} + TEST_CASE("base_power") { SECTION("base rep deducible for integral base") @@ -129,6 +135,69 @@ 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(mag_412.value, 412); + check_same_type_and_value(mag_412.value, std::size_t{412}); + check_same_type_and_value(mag_412.value, 412.0f); + check_same_type_and_value(mag_412.value, 412.0); + } + + SECTION("Negative integer powers of integer bases compute correct values") + { + constexpr auto mag_0p125 = as_magnitude(); + check_same_type_and_value(mag_0p125.value, 0.125f); + check_same_type_and_value(mag_0p125.value, 0.125); + } + + SECTION("pi to the 1 supplies correct values") + { + constexpr auto pi = pi_to_the<1>(); + check_same_type_and_value(pi.value, std::numbers::pi_v); + check_same_type_and_value(pi.value, std::numbers::pi_v); + check_same_type_and_value(pi.value, std::numbers::pi_v); + } + + SECTION("pi to arbitrary power performs computations in most accurate type at compile time") + { + 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); + constexpr auto via_long_double = static_cast(cube(std::numbers::pi_v)); + + constexpr auto pi_cubed_value = pi_cubed.value; + 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. + + // (void)as_magnitude<412>().value; + + // Would work for pow<62>: + // (void)pow<63>(as_magnitude<2>()).value; + + // Would work for pow<63>: + // (void)pow<64>(as_magnitude<2>()).value; + + (void)pow<308>(as_magnitude<10>()).value; // Compiles, correctly. + // (void)pow<309>(as_magnitude<10>()).value; + // (void)pow<3099>(as_magnitude<10>()).value; + // (void)pow<3099999>(as_magnitude<10>()).value; + + auto sqrt_2 = pow(as_magnitude<2>()); + CHECK(!is_integral(sqrt_2)); + // (void)sqrt_2.value; + } +} + TEST_CASE("Equality works for magnitudes") { SECTION("Equivalent ratios are equal") @@ -218,8 +287,50 @@ 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(as_magnitude<1>()); + check_rational_and_integral(as_magnitude<3>()); + check_rational_and_integral(as_magnitude<8>()); + check_rational_and_integral(as_magnitude<412>()); + } + + 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()); + check_rational_but_not_integral(as_magnitude()); + } +} + 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(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()") {