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<T>` represents the value of `m` in the type `T`.

If `T` is integral, we only support `m.value<T>` 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.
This commit is contained in:
Chip Hogg
2022-01-28 19:06:26 +00:00
parent 30b7250033
commit a6b6afe438
2 changed files with 225 additions and 1 deletions

View File

@@ -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>
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<T>) {
// 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<typename T> requires std::is_arithmetic_v<T>
consteval 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 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<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>)
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<To>) {
if (std::cmp_less(x, std::numeric_limits<To>::min()) ||
std::cmp_greater(x, std::numeric_limits<To>::max())) {
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.
*/
@@ -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
/**
@@ -225,7 +326,19 @@ static constexpr bool is_base_power_pack_valid = all_base_powers_valid<BPs...> &
*/
template<BasePower auto... BPs>
requires (detail::is_base_power_pack_valid<BPs...>)
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<typename T>
requires (is_integral(magnitude{}) || std::is_floating_point_v<T>)
static constexpr T value = detail::checked_static_cast<T>(
(detail::compute_base_power<T>(BPs) * ...));
};
// Implementation for Magnitude concept (below).
namespace detail {

View File

@@ -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,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<int>, 412);
check_same_type_and_value(mag_412.value<std::size_t>, std::size_t{412});
check_same_type_and_value(mag_412.value<float>, 412.0f);
check_same_type_and_value(mag_412.value<double>, 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(mag_0p125.value<float>, 0.125f);
check_same_type_and_value(mag_0p125.value<double>, 0.125);
}
SECTION("pi to the 1 supplies correct values")
{
constexpr auto pi = pi_to_the<1>();
check_same_type_and_value(pi.value<float>, std::numbers::pi_v<float>);
check_same_type_and_value(pi.value<double>, std::numbers::pi_v<double>);
check_same_type_and_value(pi.value<long double>, std::numbers::pi_v<long double>);
}
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<float>);
constexpr auto via_long_double = static_cast<float>(cube(std::numbers::pi_v<long double>));
constexpr auto pi_cubed_value = pi_cubed.value<float>;
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<int8_t>;
// Would work for pow<62>:
// (void)pow<63>(as_magnitude<2>()).value<int64_t>;
// Would work for pow<63>:
// (void)pow<64>(as_magnitude<2>()).value<uint64_t>;
(void)pow<308>(as_magnitude<10>()).value<double>; // Compiles, correctly.
// (void)pow<309>(as_magnitude<10>()).value<double>;
// (void)pow<3099>(as_magnitude<10>()).value<double>;
// (void)pow<3099999>(as_magnitude<10>()).value<double>;
auto sqrt_2 = pow<ratio{1, 2}>(as_magnitude<2>());
CHECK(!is_integral(sqrt_2));
// (void)sqrt_2.value<int>;
}
}
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<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()") {