diff --git a/src/core/include/units/magnitude.h b/src/core/include/units/magnitude.h index 5f27f66c..50dfc232 100644 --- a/src/core/include/units/magnitude.h +++ b/src/core/include/units/magnitude.h @@ -32,13 +32,42 @@ #include namespace units { + namespace detail { + // Higher numbers use fewer trial divisions, at the price of more storage space. using factorizer = wheel_factorizer<4>; + } // namespace detail /** - * @brief Any type which can be used as a basis vector in a BasePower. + * @brief A type to represent a standalone constant value. + */ +template +struct constant { + static constexpr auto value = V; +}; + +// is_derived_from_specialization_of_constant +namespace detail { + +template +void to_base_specialization_of_constant(const volatile constant*); + +template +inline constexpr bool is_derived_from_specialization_of_constant = + requires(T * t) { to_base_specialization_of_constant(t); }; + +} // namespace detail + + +template +concept Constant = detail::is_derived_from_specialization_of_constant; + +struct pi_v : constant> {}; + +/** + * @brief Any type which can be used as a basis vector in a PowerV. * * We have two categories. * @@ -52,9 +81,6 @@ using factorizer = wheel_factorizer<4>; * The reason we can't hold the value directly for floating point bases is so that we can support some compilers (e.g., * GCC 10) which don't yet permit floating point NTTPs. */ -template -concept BaseRep = - std::is_same_v || std::is_same_v, long double>; /** * @brief A basis vector in our magnitude representation, raised to some rational power. @@ -80,42 +106,36 @@ concept BaseRep = * _existing_ bases, including both prime numbers and any other irrational bases. For example, even though `sqrt(2)` is * irrational, we must not ever use it as a base; instead, we would use `base_power{2, ratio{1, 2}}`. */ -template -struct base_power { - // The rational power to which the base is raised. - ratio power{1}; +template +concept PowerVBase = Constant || std::integral; - constexpr long double get_base() const { return T::value; } +// TODO Unify with `power` if UTPs (P1985) are accepted by the Committee +template + requires(Num != 0) +struct power_v { + static constexpr auto base = V; + static constexpr ratio exponent{Num, Den...}; + static_assert(exponent != 1); }; -/** - * @brief Specialization for prime number bases. - */ -template<> -struct base_power { - // The value of the basis "vector". Must be prime to be used with `magnitude` (below). - std::intmax_t base; +namespace detail { - // The rational power to which the base is raised. - ratio power{1}; - - constexpr std::intmax_t get_base() const { return base; } -}; /** * @brief Deduction guides for base_power: only permit deducing integral bases. */ -template U> -base_power(T, U) -> base_power; -template -base_power(T) -> base_power; +// template U> +// base_power(T, U) -> base_power; +// template +// base_power(T) -> base_power; -// Implementation for BasePower concept (below). -namespace detail { +// Implementation for PowerV concept (below). template -inline constexpr bool is_base_power = false; -template -inline constexpr bool is_base_power> = true; +inline constexpr bool is_specialization_of_power_v = false; + +template +inline constexpr bool is_specialization_of_power_v> = true; + } // namespace detail /** @@ -125,130 +145,178 @@ inline constexpr bool is_base_power> = true; * `magnitude<...>`. We will defer that second check to the constraints on the `magnitude` template. */ template -concept BasePower = detail::is_base_power; +concept PowerV = detail::is_specialization_of_power_v; namespace detail { -constexpr auto inverse(BasePower auto bp) +template +[[nodiscard]] consteval auto shorten_T() { - bp.power.num *= -1; - return bp; + if constexpr (std::integral) { + if constexpr (V <= std::numeric_limits::max()) { + return static_cast(V); + } else { + return V; + } + } else { + return V; + } } +template +[[nodiscard]] consteval auto power_v_or_T() +{ + constexpr auto shortT = shorten_T(); + + if constexpr (R.den == 1) { + if constexpr (R.num == 1) + return shortT; + else + return power_v{}; + } else { + return power_v{}; + } +}; + +consteval auto inverse(PowerV auto bp) { return power_v_or_T(); } + // `widen_t` gives the widest arithmetic type in the same category, for intermediate computations. -template -using widen_t = - std::conditional_t, - std::conditional_t, long double, - std::conditional_t, std::intmax_t, std::uintmax_t>>, - T>; +// template +// using widen_t = conditional, +// conditional, long double, +// conditional, std::intmax_t, std::uintmax_t>>, +// T>; // Raise an arbitrary arithmetic type to a positive integer power at compile time. -template -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"}; - } +// template +// 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; - UNITS_DIAGNOSTIC_PUSH - UNITS_DIAGNOSTIC_IGNORE_FLOAT_EQUAL - if (result / a != b) { - throw std::overflow_error{"Wraparound detected"}; - } - UNITS_DIAGNOSTIC_POP - return result; - }; +// constexpr auto checked_multiply = [](auto a, auto b) { +// const auto result = a * b; +// UNITS_DIAGNOSTIC_PUSH +// UNITS_DIAGNOSTIC_IGNORE_FLOAT_EQUAL +// if (result / a != b) { +// throw std::overflow_error{"Wraparound detected"}; +// } +// UNITS_DIAGNOSTIC_POP +// return result; +// }; - constexpr auto checked_square = [checked_multiply](auto a) { return checked_multiply(a, a); }; +// 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. +// // 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 T{1}; - } +// if (exp == 0) { +// return T{1}; +// } - if (exp % 2 == 1) { - return checked_multiply(base, int_power(base, exp - 1)); - } +// if (exp % 2 == 1) { +// return checked_multiply(base, int_power(base, exp - 1)); +// } - return checked_square(int_power(base, exp / 2)); -} +// return checked_square(int_power(base, exp / 2)); +// } -template -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"}; - } +// template +// constexpr widen_t compute_base_power(PowerV 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.num < 0) { - if constexpr (std::is_integral_v) { - throw std::invalid_argument{"Cannot represent reciprocal as integer"}; - } else { - return T{1} / compute_base_power(inverse(bp)); - } - } +// if (bp.power.num < 0) { +// if constexpr (std::is_integral_v) { +// throw std::invalid_argument{"Cannot represent reciprocal as integer"}; +// } else { +// return T{1} / compute_base_power(inverse(bp)); +// } +// } - auto power = bp.power.num; - return int_power(static_cast>(bp.get_base()), power); -} +// auto power = bp.power.num; +// 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 -// TODO(chogg): Migrate this to use `treat_as_floating_point`. - requires(!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"}; - } - } +// template +// // TODO(chogg): Migrate this to use `treat_as_floating_point`. +// requires(!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"}; +// } +// } - return static_cast(x); -} +// return static_cast(x); +// } } // namespace detail /** * @brief Equality detection for two base powers. */ -template -constexpr bool operator==(T t, U u) +template +[[nodiscard]] consteval bool operator==(T, U) { - return std::is_same_v && (t.get_base() == u.get_base()) && (t.power == u.power); + return std::is_same_v; } -/** - * @brief A BasePower, raised to a rational power E. - */ -constexpr auto pow(BasePower auto bp, ratio p) -{ - bp.power = bp.power * p; - return bp; -} +template +concept MagnitudeSpec = std::integral || Constant || PowerV; + // A variety of implementation detail helpers. namespace detail { +template +[[nodiscard]] consteval auto get_base(Element element) +{ + if constexpr (PowerV) + return Element::base; + else + return element; +} + +template +[[nodiscard]] consteval auto get_base_value(Element element) +{ + const auto base = get_base(element); + using base_type = decltype(base); + if constexpr (std::integral) + return base; + else + return base_type::value; +} + +template +[[nodiscard]] consteval ratio get_exponent(Element) +{ + if constexpr (PowerV) + return Element::exponent; + else + return ratio{1}; +} + // The exponent of `factor` in the prime factorization of `n`. -constexpr std::intmax_t multiplicity(std::intmax_t factor, std::intmax_t n) +[[nodiscard]] consteval std::intmax_t multiplicity(std::intmax_t factor, std::intmax_t n) { std::intmax_t m = 0; while (n % factor == 0) { @@ -261,7 +329,7 @@ constexpr std::intmax_t multiplicity(std::intmax_t factor, std::intmax_t n) // Divide a number by a given base raised to some power. // // Undefined unless base > 1, pow >= 0, and (base ^ pow) evenly divides n. -constexpr std::intmax_t remove_power(std::intmax_t base, std::intmax_t pow, std::intmax_t n) +[[nodiscard]] consteval std::intmax_t remove_power(std::intmax_t base, std::intmax_t pow, std::intmax_t n) { while (pow-- > 0) { n /= base; @@ -270,15 +338,18 @@ constexpr std::intmax_t remove_power(std::intmax_t base, std::intmax_t pow, std: } // A way to check whether a number is prime at compile time. -constexpr bool is_prime(std::intmax_t n) { return (n >= 0) && factorizer::is_prime(static_cast(n)); } - -constexpr bool is_valid_base_power(const BasePower auto& bp) +[[nodiscard]] consteval bool is_prime(std::intmax_t n) { - if (bp.power == 0) { + return (n >= 0) && factorizer::is_prime(static_cast(n)); +} + +template +[[nodiscard]] consteval bool is_valid_element(Element element) +{ + if (get_exponent(element) == 0) { return false; } - - if constexpr (std::is_same_v) { + if constexpr (std::integral) { // Some prime numbers are so big, that we can't check their primality without exhausting limits on constexpr steps // and/or iterations. We can still _perform_ the factorization for these by using the `known_first_factor` // workaround. However, we can't _check_ that they are prime, because this workaround depends on the input being @@ -289,13 +360,13 @@ constexpr bool is_valid_base_power(const BasePower auto& bp) // // In our case: we simply give up on excluding every possible ill-formed base power, and settle for catching the // most likely and common mistakes. - if (const bool too_big_to_check = (bp.get_base() > 1'000'000'000)) { + if (const bool too_big_to_check = (get_base_value(element) > 1'000'000'000)) { return true; } - return is_prime(bp.get_base()); + return is_prime(get_base_value(element)); } else { - return bp.get_base() > 0; + return get_base_value(element) > 0; } } @@ -305,7 +376,7 @@ struct pairwise_all { Predicate predicate; template - constexpr bool operator()(Ts&&... ts) const + [[nodiscard]] consteval bool operator()(Ts&&... ts) const { // Carefully handle different sizes, avoiding unsigned integer underflow. constexpr auto num_comparisons = [](auto num_elements) { @@ -322,52 +393,58 @@ struct pairwise_all { }; // Deduction guide: permit constructions such as `pairwise_all{std::less{}}`. -template -pairwise_all(T) -> pairwise_all; +// template +// pairwise_all(T) -> pairwise_all; // Check whether a sequence of (possibly heterogeneously typed) values are strictly increasing. template requires(std::is_signed_v && ...) -constexpr bool strictly_increasing(Ts&&... ts) +[[nodiscard]] consteval bool strictly_increasing(Ts&&... ts) { return pairwise_all{std::less{}}(std::forward(ts)...); } -template -inline constexpr bool all_base_powers_valid = (is_valid_base_power(BPs) && ...); +template +inline constexpr bool all_elements_valid = (is_valid_element(Elements) && ...); -template -inline constexpr bool all_bases_in_order = strictly_increasing(BPs.get_base()...); +template +inline constexpr bool all_elements_in_order = strictly_increasing(get_base_value(Elements)...); -template -inline constexpr bool is_base_power_pack_valid = all_base_powers_valid && all_bases_in_order; +template +inline constexpr bool is_element_pack_valid = all_elements_valid && all_elements_in_order; -constexpr bool is_rational(BasePower auto bp) +[[nodiscard]] consteval bool is_rational(MagnitudeSpec auto element) { - return std::is_integral_v && (bp.power.den == 1); + return std::is_integral_v && (get_exponent(element).den == 1); +} + +[[nodiscard]] consteval bool is_integral(MagnitudeSpec auto element) +{ + return is_rational(element) && get_exponent(element).num > 0; } -constexpr bool is_integral(BasePower auto bp) { return is_rational(bp) && bp.power.num > 0; } } // namespace detail + /** * @brief A representation for positive real numbers which optimizes taking products and rational powers. * * Magnitudes can be treated as values. Each type encodes exactly one value. Users can multiply, divide, raise to * rational powers, and compare for equality. */ -template - requires detail::is_base_power_pack_valid +template + requires detail::is_element_pack_valid struct magnitude { // Whether this magnitude represents an integer. - friend constexpr bool is_integral(const magnitude&) { return (detail::is_integral(BPs) && ...); } + [[nodiscard]] friend consteval bool is_integral(const magnitude&) { return (detail::is_integral(Ms) && ...); } // Whether this magnitude represents a rational number. - friend constexpr bool is_rational(const magnitude&) { return (detail::is_rational(BPs) && ...); } + [[nodiscard]] friend consteval bool is_rational(const magnitude&) { return (detail::is_rational(Ms) && ...); } }; -// Implementation for Magnitude concept (below). namespace detail { + +// Implementation for Magnitude concept (below). template inline constexpr bool is_magnitude = false; template @@ -383,28 +460,21 @@ concept Magnitude = detail::is_magnitude; /** * @brief The value of a Magnitude in a desired type T. */ -template -// TODO(chogg): Migrate this to use `treat_as_floating_point`. - requires(!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) * ... * T{1})); +// template +// // TODO(chogg): Migrate this to use `treat_as_floating_point`. +// requires(!std::integral || 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) * ... * T{1})); - return result; -} - -/** - * @brief A base to represent pi. - */ -struct pi_base { - static constexpr long double value = std::numbers::pi_v; -}; +// return result; +// } /** * @brief A convenient Magnitude constant for pi, which we can manipulate like a regular number. */ -inline constexpr Magnitude auto mag_pi = magnitude{}>{}; +inline constexpr Magnitude auto mag_pi = magnitude{}; //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// // Magnitude equality implementation. @@ -422,24 +492,24 @@ constexpr bool operator==(magnitude, magnitude) //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// // Magnitude rational powers implementation. -template -constexpr auto pow(magnitude) +template +constexpr auto pow(magnitude) { if constexpr (E.num == 0) { return magnitude<>{}; } else { - return magnitude{}; + return magnitude()...>{}; } } -template -constexpr auto sqrt(magnitude m) +template +constexpr auto sqrt(magnitude m) { return pow(m); } -template -constexpr auto cbrt(magnitude m) +template +constexpr auto cbrt(magnitude m) { return pow(m); } @@ -456,8 +526,10 @@ constexpr Magnitude auto operator*(Magnitude auto m, magnitude<>) { return m; } template constexpr Magnitude auto operator*(magnitude, magnitude) { + using namespace detail; + // Case for when H1 has the smaller base. - if constexpr (H1.get_base() < H2.get_base()) { + if constexpr (get_base_value(H1) < get_base_value(H2)) { if constexpr (sizeof...(T1) == 0) { // Shortcut for the "pure prepend" case, which makes it easier to implement some of the other cases. return magnitude{}; @@ -467,21 +539,18 @@ constexpr Magnitude auto operator*(magnitude, magnitude) } // Case for when H2 has the smaller base. - if constexpr (H1.get_base() > H2.get_base()) { + if constexpr (get_base_value(H1) > get_base_value(H2)) { return magnitude

{} * (magnitude{} * magnitude{}); } // "Same leading base" case. - if constexpr (H1.get_base() == H2.get_base()) { + if constexpr (get_base(H1) == get_base(H2)) { constexpr auto partial_product = magnitude{} * magnitude{}; - // Make a new base_power with the common base of H1 and H2, whose power is their powers' sum. - constexpr auto new_head = [&](auto head) { - head.power = H1.power + H2.power; - return head; - }(H1); + // Make a new power_v with the common base of H1 and H2, whose power is their powers' sum. + constexpr auto new_head = power_v_or_T(); - if constexpr (new_head.power == 0) { + if constexpr (get_exponent(new_head) == 0) { return partial_product; } else { return magnitude{} * partial_product; @@ -507,7 +576,7 @@ constexpr auto integer_part(magnitude) constexpr auto power_den = BP.power.den; if constexpr (std::is_integral_v && (power_num >= power_den)) { - constexpr auto largest_integer_power = [=](BasePower auto bp) { + constexpr auto largest_integer_power = [=](PowerV auto bp) { bp.power = (power_num / power_den); // Note: integer division intended. return bp; }(BP); // Note: lambda is immediately invoked. @@ -529,14 +598,14 @@ constexpr auto numerator(magnitude) constexpr auto denominator(Magnitude auto m) { return numerator(pow<-1>(m)); } // Implementation of conversion to ratio goes here, because it needs `numerator()` and `denominator()`. -constexpr ratio as_ratio(Magnitude auto m) - requires(is_rational(decltype(m){})) -{ - return ratio{ - get_value(numerator(m)), - get_value(denominator(m)), - }; -} +// constexpr ratio as_ratio(Magnitude auto m) +// requires(is_rational(decltype(m){})) +// { +// return ratio{ +// get_value(numerator(m)), +// get_value(denominator(m)), +// }; +// } //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// @@ -616,6 +685,7 @@ template inline constexpr std::optional known_first_factor = std::nullopt; namespace detail { + // Helper to perform prime factorization at compile time. template requires(N > 0) @@ -634,7 +704,7 @@ struct prime_factorization { static constexpr std::intmax_t remainder = remove_power(first_base, first_power, N); static constexpr auto value = - magnitude{} * prime_factorization::value; + magnitude()>{} * prime_factorization::value; }; // Specialization for the prime factorization of 1 (base case). @@ -666,7 +736,7 @@ template inline constexpr Magnitude auto mag_power = pow(mag); namespace detail { -template +template constexpr ratio get_power(T base, magnitude) { return ((BPs.get_base() == base ? BPs.power : ratio{0}) + ... + ratio{0});