diff --git a/au/BUILD.bazel b/au/BUILD.bazel index 042a8c77..1d7dbc85 100644 --- a/au/BUILD.bazel +++ b/au/BUILD.bazel @@ -105,6 +105,24 @@ cc_test( ################################################################################ # Implementation detail libraries and tests +cc_library( + name = "apply_magnitude", + hdrs = ["apply_magnitude.hh"], + deps = [":magnitude"], +) + +cc_test( + name = "apply_magnitude_test", + size = "small", + srcs = ["apply_magnitude_test.cc"], + copts = ["-Iexternal/gtest/include"], + deps = [ + ":apply_magnitude", + ":testing", + "@com_google_googletest//:gtest_main", + ], +) + cc_library( name = "chrono_policy_validation", testonly = True, @@ -304,6 +322,7 @@ cc_library( name = "quantity", hdrs = ["quantity.hh"], deps = [ + ":apply_magnitude", ":conversion_policy", ":operators", ":unit_of_measure", diff --git a/au/apply_magnitude.hh b/au/apply_magnitude.hh new file mode 100644 index 00000000..798831f5 --- /dev/null +++ b/au/apply_magnitude.hh @@ -0,0 +1,115 @@ +// Copyright 2023 Aurora Operations, Inc. +// +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. + +#pragma once + +#include "au/magnitude.hh" + +namespace au { +namespace detail { + +// The various categories by which a magnitude can be applied to a numeric quantity. +enum class ApplyAs { + INTEGER_MULTIPLY, + INTEGER_DIVIDE, + RATIONAL_MULTIPLY, + IRRATIONAL_MULTIPLY, +}; + +template +constexpr ApplyAs categorize_magnitude(Magnitude) { + if (IsInteger>::value) { + return ApplyAs::INTEGER_MULTIPLY; + } + + if (IsInteger>>::value) { + return ApplyAs::INTEGER_DIVIDE; + } + + return IsRational>::value ? ApplyAs::RATIONAL_MULTIPLY + : ApplyAs::IRRATIONAL_MULTIPLY; +} + +template +struct ApplyMagnitudeImpl; + +// Multiplying by an integer, for any type T. +template +struct ApplyMagnitudeImpl { + static_assert(categorize_magnitude(Mag{}) == ApplyAs::INTEGER_MULTIPLY, + "Mismatched instantiation (should never be done manually)"); + static_assert(is_T_integral == std::is_integral::value, + "Mismatched instantiation (should never be done manually)"); + + constexpr T operator()(const T &x) { return x * get_value(Mag{}); } +}; + +// Dividing by an integer, for any type T. +template +struct ApplyMagnitudeImpl { + static_assert(categorize_magnitude(Mag{}) == ApplyAs::INTEGER_DIVIDE, + "Mismatched instantiation (should never be done manually)"); + static_assert(is_T_integral == std::is_integral::value, + "Mismatched instantiation (should never be done manually)"); + + constexpr T operator()(const T &x) { return x / get_value(MagInverseT{}); } +}; + +// Applying a (non-integer, non-inverse-integer) rational, for any integral type T. +template +struct ApplyMagnitudeImpl { + static_assert(categorize_magnitude(Mag{}) == ApplyAs::RATIONAL_MULTIPLY, + "Mismatched instantiation (should never be done manually)"); + static_assert(std::is_integral::value, + "Mismatched instantiation (should never be done manually)"); + + constexpr T operator()(const T &x) { + return x * get_value(numerator(Mag{})) / get_value(denominator(Mag{})); + } +}; + +// Applying a (non-integer, non-inverse-integer) rational, for any non-integral type T. +template +struct ApplyMagnitudeImpl { + static_assert(categorize_magnitude(Mag{}) == ApplyAs::RATIONAL_MULTIPLY, + "Mismatched instantiation (should never be done manually)"); + static_assert(!std::is_integral::value, + "Mismatched instantiation (should never be done manually)"); + + constexpr T operator()(const T &x) { return x * get_value(Mag{}); } +}; + +// Applying an irrational for any type T (although only non-integral T makes sense). +template +struct ApplyMagnitudeImpl { + static_assert(!std::is_integral::value, "Cannot apply irrational magnitude to integer type"); + + static_assert(categorize_magnitude(Mag{}) == ApplyAs::IRRATIONAL_MULTIPLY, + "Mismatched instantiation (should never be done manually)"); + static_assert(is_T_integral == std::is_integral::value, + "Mismatched instantiation (should never be done manually)"); + + constexpr T operator()(const T &x) { return x * get_value(Mag{}); } +}; + +template +constexpr T apply_magnitude(const T &x, Magnitude m) { + return ApplyMagnitudeImpl, + categorize_magnitude(m), + T, + std::is_integral::value>{}(x); +} + +} // namespace detail +} // namespace au diff --git a/au/apply_magnitude_test.cc b/au/apply_magnitude_test.cc new file mode 100644 index 00000000..d2ca1538 --- /dev/null +++ b/au/apply_magnitude_test.cc @@ -0,0 +1,132 @@ +// Copyright 2023 Aurora Operations, Inc. +// +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. + +#include "au/apply_magnitude.hh" + +#include "au/testing.hh" +#include "gtest/gtest.h" + +using ::testing::ElementsAreArray; +using ::testing::Not; + +namespace au { +namespace detail { +namespace { +template +std::vector first_n_positive_values(std::size_t n) { + std::vector result; + result.reserve(n); + for (auto i = 1u; i <= n; ++i) { + result.push_back(static_cast(i)); + } + return result; +} +} // namespace + +TEST(CategorizeMagnitude, FindsIntegerMultiplyInstances) { + EXPECT_EQ(categorize_magnitude(mag<2>()), ApplyAs::INTEGER_MULTIPLY); + EXPECT_EQ(categorize_magnitude(mag<35>()), ApplyAs::INTEGER_MULTIPLY); + + EXPECT_EQ(categorize_magnitude(mag<35>() / mag<7>()), ApplyAs::INTEGER_MULTIPLY); +} + +TEST(CategorizeMagnitude, FindsIntegerDivideInstances) { + EXPECT_EQ(categorize_magnitude(ONE / mag<2>()), ApplyAs::INTEGER_DIVIDE); + EXPECT_EQ(categorize_magnitude(ONE / mag<35>()), ApplyAs::INTEGER_DIVIDE); + + EXPECT_EQ(categorize_magnitude(mag<7>() / mag<35>()), ApplyAs::INTEGER_DIVIDE); +} + +TEST(CategorizeMagnitude, FindsRationalMultiplyInstances) { + EXPECT_EQ(categorize_magnitude(mag<5>() / mag<2>()), ApplyAs::RATIONAL_MULTIPLY); +} + +TEST(CategorizeMagnitude, FindsIrrationalMultiplyInstances) { + EXPECT_EQ(categorize_magnitude(sqrt(mag<2>())), ApplyAs::IRRATIONAL_MULTIPLY); + EXPECT_EQ(categorize_magnitude(PI), ApplyAs::IRRATIONAL_MULTIPLY); +} + +TEST(ApplyMagnitude, MultipliesForIntegerMultiply) { + constexpr auto m = mag<25>(); + ASSERT_EQ(categorize_magnitude(m), ApplyAs::INTEGER_MULTIPLY); + + EXPECT_THAT(apply_magnitude(4, m), SameTypeAndValue(100)); + EXPECT_THAT(apply_magnitude(4.0f, m), SameTypeAndValue(100.0f)); +} + +TEST(ApplyMagnitude, DividesForIntegerDivide) { + constexpr auto one_thirteenth = ONE / mag<13>(); + ASSERT_EQ(categorize_magnitude(one_thirteenth), ApplyAs::INTEGER_DIVIDE); + + // This test would fail if our implementation multiplied by the float representation of (1/13), + // instead of dividing by 13, under the hood. + for (const auto &i : first_n_positive_values(100u)) { + EXPECT_THAT(apply_magnitude(i * 13, one_thirteenth), SameTypeAndValue(i)); + } +} + +TEST(ApplyMagnitude, MultipliesThenDividesForRationalMagnitudeOnInteger) { + // Consider applying the magnitude (3/2) to the value 5. The exact answer is the real number + // 7.5, which becomes 7 when translated (via truncation) to the integer domain. + // + // If we multiply-then-divide, we get (5 * 3) / 2 = 7, which is correct. + // + // If we divide-then-multiply --- say, because we are trying to avoid overflow --- then we get + // (5 / 2) * 3 = 2 * 3 = 6, which is wrong. + constexpr auto three_halves = mag<3>() / mag<2>(); + ASSERT_EQ(categorize_magnitude(three_halves), ApplyAs::RATIONAL_MULTIPLY); + + EXPECT_THAT(apply_magnitude(5, three_halves), SameTypeAndValue(7)); +} + +TEST(ApplyMagnitude, MultipliesSingleNumberForRationalMagnitudeOnFloatingPoint) { + // Helper similar to `std::transform`, but with more convenient interfaces. + auto apply = [](std::vector vals, auto fun) { + for (auto &v : vals) { + v = fun(v); + } + return vals; + }; + + // Create our rational magnitude, (2 / 13). + constexpr auto two_thirteenths = mag<2>() / mag<13>(); + ASSERT_EQ(categorize_magnitude(two_thirteenths), ApplyAs::RATIONAL_MULTIPLY); + + // Test a bunch of values. We are hoping that the two different strategies will yield different + // results for at least some of these strategies (and we'll check that this is the case). + const auto original_vals = first_n_positive_values(10u); + + // Compute expected answers for each possible strategy. + const auto if_we_multiply_and_divide = + apply(original_vals, [](float v) { return v * 2.0f / 13.0f; }); + const auto if_we_use_one_factor = + apply(original_vals, [](float v) { return v * (2.0f / 13.0f); }); + + // The strategies must be different for at least some results! + ASSERT_THAT(if_we_multiply_and_divide, Not(ElementsAreArray(if_we_use_one_factor))); + + // Make sure we follow the single-number strategy, every time. + const auto results = + apply(original_vals, [=](float v) { return apply_magnitude(v, two_thirteenths); }); + EXPECT_THAT(results, ElementsAreArray(if_we_use_one_factor)); + EXPECT_THAT(results, Not(ElementsAreArray(if_we_multiply_and_divide))); +} + +TEST(ApplyMagnitude, MultipliesSingleNumberForIrrationalMagnitudeOnFloatingPoint) { + ASSERT_EQ(categorize_magnitude(PI), ApplyAs::IRRATIONAL_MULTIPLY); + EXPECT_THAT(apply_magnitude(2.0f, PI), SameTypeAndValue(2.0f * static_cast(M_PI))); +} + +} // namespace detail +} // namespace au diff --git a/au/quantity.hh b/au/quantity.hh index 5e216378..eb4f7d4f 100644 --- a/au/quantity.hh +++ b/au/quantity.hh @@ -16,6 +16,7 @@ #include +#include "au/apply_magnitude.hh" #include "au/conversion_policy.hh" #include "au/operators.hh" #include "au/stdx/functional.hh" @@ -143,17 +144,11 @@ class Quantity { typename NewUnit, typename = std::enable_if_t::value>> constexpr auto as(NewUnit) const { - constexpr auto ratio = unit_ratio(unit, NewUnit{}); - using Common = std::common_type_t; - constexpr auto NUM = integer_part(numerator(ratio)); - constexpr auto DEN = integer_part(denominator(ratio)); - constexpr auto num = get_value(NUM); - constexpr auto den = get_value(DEN); - constexpr auto irr = get_value(ratio * DEN / NUM); + using Factor = UnitRatioT; return make_quantity( - static_cast(static_cast(value_) * num / den * irr)); + static_cast(detail::apply_magnitude(static_cast(value_), Factor{}))); } template ::value>> diff --git a/au/quantity_test.cc b/au/quantity_test.cc index 282eb73d..50d0094f 100644 --- a/au/quantity_test.cc +++ b/au/quantity_test.cc @@ -94,6 +94,14 @@ TEST(MakeQuantity, MakesQuantityInGivenUnit) { EXPECT_EQ(make_quantity(99), feet(99)); } +TEST(Quantity, RationalConversionRecoversExactIntegerValues) { + // This test would fail if our implementation multiplied by the float + // representation of (1/13), instead of dividing by 13, under the hood. + for (int i = 1; i < 100; ++i) { + EXPECT_EQ(feet(static_cast(i * 13)).in(feet * mag<13>()), i); + } +} + TEST(QuantityMaker, CreatesAppropriateQuantityIfCalled) { EXPECT_EQ(yards(3.14).in(yards), 3.14); } TEST(QuantityMaker, CanBeMultipliedBySingularUnitToGetMakerOfProductUnit) { diff --git a/docs/discussion/implementation/applying_magnitudes.md b/docs/discussion/implementation/applying_magnitudes.md new file mode 100644 index 00000000..1cdbfca0 --- /dev/null +++ b/docs/discussion/implementation/applying_magnitudes.md @@ -0,0 +1,137 @@ +# Applying Magnitudes + +Every unit conversion factor is a [_magnitude_](../../reference/magnitude.md) --- a positive real +number. When we apply it to a value, conceptually, we're just multiplying the value by that number. +However, that doesn't mean that multiplying by a number is the best _implementation_! Consider +these examples. + +- **Factor: $\frac{5}{8}$; value: `12` (type: `int`).** + - Computationally, we don't want to leave the integral domain. But of course, $\frac{5}{8}$ + can't be represented as an integer! This suggests we should perform _two_ operations: first + multiply by `5`, and then divide by `8`, yielding `7`. +- **Factor: $\frac{1}{13}$; value: `91.0f` (type: `float`).** + - Conceptually, this has an exactly representable answer: $91\left(\frac{1}{13}\right) = 7$. + However, if we multiply by the single number `(1.0f / 13.0f)`, we obtain the approximation + `7.0000004768371582`! This suggests that for $\frac{1}{13}$, at least, it would be better to + divide by `13.0f`. + +Au is thoughtful about how we apply conversion factors. We first compute a _category_ for the +factor, which dictates the best strategy for applying it. We may also take into account whether +we're dealing with an integral or floating point type. + +## Magnitude categories + +We represent conversion factors with [magnitudes](../../reference/magnitude.md). These +representations support _exact symbolic math_ for products and rational powers. They also support +querying for _numeric properties_ of the number, such as whether it is an integer, whether it's +irrational, and so on. + +For purposes of _applying to a value_, we find four useful categories of magnitude. + +1. Integers. +2. Reciprocal integers. +3. Rational numbers (other than the first two categories). +4. Irrational numbers. + +These categories are mutually exclusive and exhaustive. Below, we'll explain the best strategy for +each one. + +### Integers + +Applying an integer magnitude to a type `T` is simple: we multiply by that integer's representation +in `T`. + +This always compiles to a single instruction, and always produces exact answers whenever they are +representable in the type `T`. + +### Reciprocal integers + +If a magnitude is _not_ an integer, but its _reciprocal is_, then we divide by its reciprocal. For +example, in converting a value from `inches` to `feet`, we will divide by $12$, instead of +multiplying by the representation of $\frac{1}{12}$, which would be inexact. + +As with integers, this always compiles to a single instruction, and always produces exact answers +whenever they are representable in the type `T`. + +### Rational numbers + +Again, to be clear, this category only includes rationals that are _neither_ integers _nor_ +reciprocal integers. So, for example, neither $2$ nor $\frac{1}{5}$ falls in this category, but +$\frac{2}{5}$ does. + +This category is interesting, because it's the first instance where our strategy _depends on the +type `T`_ to which we're applying the factor. The best approach differs between integral and +non-integral types. + +#### Integral types + +Here, we multiply by the numerator, then divide by the denominator. This compiles to _two_ +operations instead of one, but it's the only way to get reasonable accuracy. + +There's another issue: the multiplication operation can overflow. This means we can produce wrong +answers in some instances, even when the correct answer is representable in the type! For example, +let's say our value is `std::numeric_limits::max()`, and we apply the magnitude +$\frac{2}{3}$: by the time we divide by $3$, the multiplication by $2$ has already lost our value to +overflow. + +We might be tempted to prevent this by doing the division first. In the above example, this would +certainly give us a much closer result! However, the cost would be reduced accuracy for _smaller_ +values, which are far more common. Consider applying $\frac{2}{3}$ to a smaller number, such as +`5`. The exact rational answer is $\frac{10}{3}$, which truncates to `3`. If we perform the +multiplication first, this is what we get, but doing the division first would give `2`. + +If you know that your final answer is representable, _and_ you have an integer type with more bits +than your type `T`, then you can work around this issue manually by casting to the wider type, +applying the magnitude, and casting back to `T`. However, if you _don't_ have a wider integer types, +we know of no _general_ "solution" that wouldn't do more harm then good. + +#### Floating point types + +Applying a rational magnitude $\frac{N}{D}$ to a value of floating point type `T` presents a genuine +tradeoff. On the one hand, we could take the same approach as for the integers, and perform two +operations: multiplying by $N$, then dividing by $D$. On the other hand, we could simply multiply +by the single number which best represents $\frac{N}{D}$. Here's a summary of the tradeoffs: + +Criterion | Weighting | Multiply-and-divide: `(val * N) / D` | Single number: `val * (N / D)` +---|---|---|--- +Instructions | medium | 2 | 1 +Overflow | low | More vulnerable | Less vulnerable +Exact answers for multiples of $D$ | low | Guaranteed | Not guaranteed + +Overall, we aren't worried much about missing out on exact answers. Users of floating point know +they need to handle the possibility that a calculation's result can be one or two representable +values away from the best possible result. (This is commonly called the "usual floating point +error".) + +We also aren't very worried about overflow. Even float has a range of $10^{38}$, while going from +femtometers to Astronomical Units (AU) spans a range of "only" about $10^{26}$. + +Going from 1 instruction to 2 is a moderate concern, which means that it outweighs the other two +considerations. It represents a runtime penalty relative to the usual approach people take without +a units library, which is to compute a single conversion factor. We always strive to avoid runtime +penalties in units libraries! The reason we don't consider this even more serious is that unit +conversions should never occur in the "hot loop" for a program; thus, this performance hit isn't +really meaningful. + +**Outcome:** we represent a rational conversion factor $\frac{N}{D}$ with a **single number** when +applying it to a floating point variable. + +### Irrational numbers + +There is no reason to try splitting an irrational number into parts to get an exact answer. Since +we're multiplying our variable by an irrational number, we know the result won't be exactly +representable. Therefore, we always simply multiply by the closest representation of this +conversion factor. + +The one difference is that we forbid this operation for integral types, because it makes no sense. + +## Summary and conclusion + +Applying a conversion factor to a numeric variable of type `T` can be a tricky and subtle business. +Au takes a thoughtful, tailored approach, which can be summarized as follows: + +- If the conversion factor multiplies --- _or divides_ --- by an exact integer, then we do that. +- Otherwise, if it's a rational number $\frac{N}{D}$, and `T` is integral, then we multiply by $N$ + and divide by $D$ (each represented in `T`). +- Otherwise, we simply multiply by the nearest representation of the conversion factor in `T` --- + with the exception that if `T` is integral, we raise a compiler error for irrational factors. diff --git a/docs/discussion/implementation/index.md b/docs/discussion/implementation/index.md index 15c1ca28..5c0980ad 100644 --- a/docs/discussion/implementation/index.md +++ b/docs/discussion/implementation/index.md @@ -3,6 +3,11 @@ This section provides in-depth discussion about core implementation details of Au. If you want to understand _how_ the library works, this is a good place to go. Here's a rough guide. +- **[Applying Magnitudes](./applying_magnitudes.md)**. A conversion factor is a _magnitude_: + a positive real number. The best way to apply it to a value of type `T` depends on what kind of + number it is (integer, rational, and so on), the runtime performance, and whether we can get exact + answers. + - **[Vector Space Representations](./vector_space.md)**. We're not talking about position vectors or velocity vectors! There's a different kind of vector at the heart of every units library. This is the core foundational concept on which we built Au's implementation.