diff --git a/src/core/include/units/magnitude.h b/src/core/include/units/magnitude.h index 61d7e1a7..77de58b2 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 +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 requires std::is_arithmetic_v +constexpr 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 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) { + 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) +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) { + if (!std::in_range(x)) { + 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. */ @@ -202,7 +295,7 @@ pairwise_all(T) -> pairwise_all; // Check whether a sequence of (possibly heterogeneously typed) values are strictly increasing. template - requires ((std::is_signed_v && ...)) + requires (std::is_signed_v && ...) constexpr bool strictly_increasing(Ts&&... ts) { return pairwise_all{std::less{}}(std::forward(ts)...); } @@ -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 /** @@ -224,8 +325,14 @@ static constexpr bool is_base_power_pack_valid = all_base_powers_valid & * rational powers, and compare for equality. */ template - requires (detail::is_base_power_pack_valid) -struct magnitude {}; + requires detail::is_base_power_pack_valid +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> = true; template concept Magnitude = detail::is_magnitude; +/** + * @brief The value of a Magnitude in a desired type T. + */ +template + requires std::is_floating_point_v || (std::is_integral_v && is_integral(magnitude{})) +constexpr T get_value(const magnitude &) { + // Force the expression to be evaluated in a constexpr context, to catch, e.g., overflow. + constexpr auto result = detail::checked_static_cast((detail::compute_base_power(BPs) * ...)); + + return result; +} + /** * @brief A base to represent pi. */ diff --git a/test/unit_test/runtime/magnitude_test.cpp b/test/unit_test/runtime/magnitude_test.cpp index 653a6ef8..76dcb16c 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,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(mag_412), 412); + check_same_type_and_value(get_value(mag_412), std::size_t{412}); + check_same_type_and_value(get_value(mag_412), 412.0f); + check_same_type_and_value(get_value(mag_412), 412.0); + } + + SECTION("Negative integer powers of integer bases compute correct values") + { + constexpr auto mag_0p125 = as_magnitude(); + check_same_type_and_value(get_value(mag_0p125), 0.125f); + check_same_type_and_value(get_value(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(pi), std::numbers::pi_v); + check_same_type_and_value(get_value(pi), std::numbers::pi_v); + check_same_type_and_value(get_value(pi), std::numbers::pi_v); + } + + 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); + constexpr auto via_long_double = static_cast(cube(std::numbers::pi_v)); + + constexpr auto pi_cubed_value = get_value(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(as_magnitude<412>()); + + // Would work for pow<62>: + // get_value(pow<63>(as_magnitude<2>())); + + // Would work for pow<63>: + // get_value(pow<64>(as_magnitude<2>())); + + get_value(pow<308>(as_magnitude<10>())); // Compiles, correctly. + // get_value(pow<309>(as_magnitude<10>())); + // get_value(pow<3099>(as_magnitude<10>())); + // get_value(pow<3099999>(as_magnitude<10>())); + + auto sqrt_2 = pow(as_magnitude<2>()); + CHECK(!is_integral(sqrt_2)); + // get_value(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()); + } + + 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()") {