Make Magnitude sum N-ary and more general (#652)

When we start collecting like terms as part of adding `Constant`
instances to each other, we will naturally form N-ary sums.  To prepare
the way, we introduce a `MagSum<Ms...>` utility to perform this
operation.

Now, we _could_ just implement this in terms of the existing `operator+`
for `Magnitude`.  However, that gets a little weird: it makes the result
depend on the order.  It's possible for a subset of inputs to overflow
(thus making the binary operation return an error), when subsequent
inputs might "rescue" the result (say, by subtracting to bring it back
into range).  Thus, I redid the implementation from scratch, using
modular arithmetic for the result, and keeping track of the (signed)
number of overflows.  This gives us exact information about the true sum
(over a very very wide range of values), although we can't always _use_
this information.  Specifically, we can only use the result when the
overflow is 0 (simple), or when it is -1 _and_ the sum is not 0 (because
if sum _were_ 0 in this case it would actually represent -2^64).

This new, more symmetrical implementation immediately takes the place of
all three `Magnitude`-`Magnitude` implementations of `operator+`.

Helps #607.

---------

Co-authored-by: Michael Hordijk <hordijk@aurora.tech>

Author: chiphogg

au/magnitude.hh

@@ -275,6 +275,19 @@ using CommonMagnitude = typename CommonMagnitudeImpl<Ms...>::type;
 template <typename... Ms>
 using CommonMagnitudeT = CommonMagnitude<Ms...>;
 
+// The sum of arbitrarily many `Magnitude` and/or `Zero` types.
+//
+// We only support this when it is "easy" to compute, where "easy" is defined as:
+// 1) all inputs being expressible as integer multiples of some common factor;
+// 2) each such integer's absolute value fitting in a `uint64_t`;
+// 3) *and*, the absolute value of the sum also fitting in a `uint64_t`.
+//
+// For all other cases, we currently produce a compile time error.
+template <typename... Ms>
+struct MagSumImpl;
+template <typename... Ms>
+using MagSum = typename MagSumImpl<Ms...>::type;
+
 ////////////////////////////////////////////////////////////////////////////////////////////////////
 // Value based interface for Magnitude (and Zero).
 
@@ -513,45 +526,8 @@ constexpr Zero mag_round(Zero) { return {}; }
 
 // Addition:
 template <typename... BP1s, typename... BP2s>
-constexpr auto operator+(Magnitude<BP1s...> m1, Magnitude<BP2s...> m2) {
-    constexpr auto sgn1 = sign(m1);
-    constexpr auto sgn2 = sign(m2);
-
-    constexpr auto abs_common = Abs<CommonMagnitude<Magnitude<BP1s...>, Magnitude<BP2s...>>>{};
-    constexpr auto abs_num1 = abs(m1) / abs_common;
-    constexpr auto abs_num2 = abs(m2) / abs_common;
-
-    // These `get_value` calls automatically check that individual _inputs_ fit in `uint64_t`.
-    constexpr auto abs_num1_u64 = get_value<std::uint64_t>(abs_num1);
-    constexpr auto abs_num2_u64 = get_value<std::uint64_t>(abs_num2);
-
-    // Biggest absolute input determines overall sign.
-    //
-    // Note that when the magnitudes are equal, either the choice doesn't matter (when the inputs
-    // have the same sign), or the outcome should just be `Zero`.  In the latter case, we rely on
-    // the explicit `Negative` overloads below being a better match.
-    constexpr auto sgn =
-        std::conditional_t<(abs_num1_u64 > abs_num2_u64), decltype(sgn1), decltype(sgn2)>{};
-
-    // Here, we are taking advantage of modular arithmetic on unsigned integers.  This actually does
-    // handle all the signs correctly, although it may not be obvious at first glance.
-    constexpr auto num1_u64 = is_positive(sgn1) ? abs_num1_u64 : -abs_num1_u64;
-    constexpr auto num2_u64 = is_positive(sgn2) ? abs_num2_u64 : -abs_num2_u64;
-    constexpr auto abs_sum_u64 = is_positive(sgn) ? (num1_u64 + num2_u64) : -(num1_u64 + num2_u64);
-
-    // Here is where we guard against overflow in the _output_.
-    static_assert((sgn1 != sgn2) || abs_sum_u64 >= abs_num1_u64,
-                  "Magnitude addition overflowed uint64_t");
-
-    return sgn * mag<abs_sum_u64>() * abs_common;
-}
-template <typename... BPs>
-constexpr Zero operator+(Magnitude<Negative, BPs...>, Magnitude<BPs...>) {
-    return {};
-}
-template <typename... BPs>
-constexpr Zero operator+(Magnitude<BPs...>, Magnitude<Negative, BPs...>) {
-    return {};
+constexpr auto operator+(Magnitude<BP1s...>, Magnitude<BP2s...>) {
+    return MagSum<Magnitude<BP1s...>, Magnitude<BP2s...>>{};
 }
 template <typename... BPs>
 constexpr auto operator+(Zero, Magnitude<BPs...> m) {
@@ -1276,4 +1252,80 @@ struct CommonMagnitudeImpl<Zero, M> : stdx::type_identity<M> {};
 template <>
 struct CommonMagnitudeImpl<Zero, Zero> : stdx::type_identity<Zero> {};
 
+////////////////////////////////////////////////////////////////////////////////////////////////////
+// `MagSum` implementation.
+
+namespace detail {
+
+// `U64MagSum<Ms...>` is the `Magnitude` (or `Zero`) equal to the sum of the Magnitudes `Ms...`, as
+// long as these preconditions are met:
+//
+// 1. The absolute value of each member of `Ms...` fits in a `std::uint64_t`.
+// 2. The absolute value of the sum of all members of `Ms...` fits in a `std::uint64_t`.
+template <typename... Ms>
+struct U64MagSumImpl {
+    struct U64SumResult {
+        std::uint64_t sum = 0u;
+        int overflow = 0;
+    };
+
+    static constexpr U64SumResult compute() {
+        const std::uint64_t abs_values[] = {get_value<std::uint64_t>(Abs<Ms>{})...};
+        const int overflows[] = {(IsPositive<Ms>::value ? 0 : -1)...};
+
+        U64SumResult result = {0u, 0};
+        for (std::size_t i = 0u; i < sizeof...(Ms); ++i) {
+            std::uint64_t old_sum = result.sum;
+            result.sum += (overflows[i] >= 0) ? abs_values[i] : -abs_values[i];
+            result.overflow += overflows[i] + (result.sum < old_sum);
+        }
+
+        return result;
+    }
+    static constexpr std::uint64_t sum = compute().sum;
+    static constexpr int overflow = compute().overflow;
+
+    static_assert((overflow == 0) || (overflow == -1 && sum > 0u),
+                  "Magnitude sum overflowed uint64_t");
+
+    using Sign = std::conditional_t<(overflow == -1), Magnitude<Negative>, Magnitude<>>;
+
+    using AbsMag = std::conditional_t<(overflow == 0) && (sum == 0u),
+                                      Zero,
+                                      // The surprising `sum == 0u` avoids asking for `mag<0>()`.
+                                      // It's fine, because it can never actually be used.
+                                      decltype(mag<(overflow == -1 ? -sum : sum) + (sum == 0u)>())>;
+
+    using type = decltype(Sign{} * AbsMag{});
+};
+template <typename... Ms>
+using U64MagSum = typename U64MagSumImpl<Ms...>::type;
+
+template <typename... Ms>
+constexpr std::uint64_t U64MagSumImpl<Ms...>::sum;
+
+template <typename... Ms>
+constexpr int U64MagSumImpl<Ms...>::overflow;
+
+template <typename... Ms>
+struct MagSumImplHelper {
+    using Common = CommonMagnitude<Ms...>;
+    using type = decltype(Common{} * U64MagSum<decltype(Ms{} / Common{})...>{});
+};
+
+// The sum of no things is nothing.
+template <>
+struct MagSumImplHelper<> : stdx::type_identity<Zero> {};
+
+// Keep stripping off zeros until we find at least one nonzero element, so that the common magnitude
+// machinery can find something meaningful and avoid dividing by zero.  (Zeros in the middle will be
+// automatically handled correctly as long as there are nonzero elements.)
+template <typename... Ms>
+struct MagSumImplHelper<Zero, Ms...> : MagSumImplHelper<Ms...> {};
+
+}  // namespace detail
+
+template <typename... Ms>
+struct MagSumImpl : detail::MagSumImplHelper<Ms...> {};
+
 }  // namespace  au

au/magnitude_test.cc

@@ -1216,6 +1216,60 @@ TEST(IsMagnitudeU64RationalCompatible, FalseForIrrationals) {
     EXPECT_THAT(is_magnitude_u64_rational_compatible(sqrt(mag<2>())), IsFalse());
 }
 
+TEST(MagSum, SumsArbitraryNumberOfMagnitudes) {
+    StaticAssertTypeEq<MagSum<decltype(mag<1>()), decltype(mag<2>()), decltype(mag<3>())>,
+                       decltype(mag<6>())>();
+
+    StaticAssertTypeEq<
+        MagSum<decltype(mag<10>()), decltype(mag<20>()), decltype(mag<30>()), decltype(mag<40>())>,
+        decltype(mag<100>())>();
+}
+
+TEST(MagSum, HandlesNegativeMagnitudes) {
+    StaticAssertTypeEq<MagSum<decltype(mag<5>()), decltype(-mag<3>())>, decltype(mag<2>())>();
+
+    StaticAssertTypeEq<MagSum<decltype(-mag<5>()), decltype(mag<3>())>, decltype(-mag<2>())>();
+
+    StaticAssertTypeEq<MagSum<decltype(mag<10>()), decltype(-mag<3>()), decltype(-mag<2>())>,
+                       decltype(mag<5>())>();
+}
+
+TEST(MagSum, ProducesZeroWhenInputsCancelOut) {
+    StaticAssertTypeEq<MagSum<decltype(mag<5>()), decltype(-mag<5>())>, Zero>();
+
+    StaticAssertTypeEq<MagSum<decltype(mag<10>()), decltype(-mag<3>()), decltype(-mag<7>())>,
+                       Zero>();
+}
+
+TEST(MagSum, HandlesFractions) {
+    StaticAssertTypeEq<MagSum<decltype(mag<1>() / mag<6>()),
+                              decltype(mag<1>() / mag<3>()),
+                              decltype(mag<1>() / mag<2>())>,
+                       decltype(mag<1>())>();
+}
+
+TEST(MagSum, IgnoresZeroInputs) {
+    StaticAssertTypeEq<MagSum<Zero, decltype(mag<5>())>, decltype(mag<5>())>();
+    StaticAssertTypeEq<MagSum<decltype(mag<5>()), Zero>, decltype(mag<5>())>();
+    StaticAssertTypeEq<MagSum<Zero, decltype(mag<5>()), Zero>, decltype(mag<5>())>();
+    StaticAssertTypeEq<MagSum<Zero, Zero, decltype(mag<5>())>, decltype(mag<5>())>();
+    StaticAssertTypeEq<MagSum<decltype(mag<2>()), Zero, decltype(mag<3>())>, decltype(mag<5>())>();
+    StaticAssertTypeEq<MagSum<Zero>, Zero>();
+    StaticAssertTypeEq<MagSum<Zero, Zero>, Zero>();
+    StaticAssertTypeEq<MagSum<Zero, Zero, Zero>, Zero>();
+}
+
+TEST(MagSum, ResultIsIndependentOfInputOrder) {
+    // A binary fold could only handle this in certain orders; N-ary handles all orders.
+    using One = decltype(mag<1>());
+    using TwoTo63 = decltype(pow<63>(mag<2>()));
+    using NegTwo = decltype(-mag<2>());
+
+    // 1 + 2^63 + 2^63 - 2 = 2^64 - 1
+    using Expected = decltype(mag<-uint64_t{1}>());
+    StaticAssertTypeEq<MagSum<One, TwoTo63, TwoTo63, NegTwo>, Expected>();
+}
+
 TEST(AssertMagnitudeU64RationalCompatible, NoCompilerErrorForKnownValidInput) {
     (void)AssertMagnitudeU64RationalCompatible<decltype(mag<1>())>{};
     (void)AssertMagnitudeU64RationalCompatible<decltype(mag<1000>())>{};