Improve overflow bounds for rational scaling factors by using modular decomposition

au/abstract_operations.hh

@@ -78,6 +78,13 @@ struct OpSequenceImpl;
 template <typename... Ops>
 using OpSequence = FlattenAs<OpSequenceImpl, Ops...>;
 
+//
+// `ScaleByRational<T, Num, Den>` represents an operation that scales a value of type `T` by the
+// rational number `Num/Den`, taking care to avoid overflow if possible.
+//
+template <typename T, typename Num, typename Den>
+struct ScaleByRational;
+
 ////////////////////////////////////////////////////////////////////////////////////////////////////
 // IMPLEMENTATION DETAILS (`abstract_operations.hh`):
 ////////////////////////////////////////////////////////////////////////////////////////////////////
@@ -188,5 +195,63 @@ struct OpSequenceImpl<Op, Ops...> {
     }
 };
 
+////////////////////////////////////////////////////////////////////////////////////////////////////
+// `ScaleByRational<T, Num, Den>` implementation.
+//
+// Modular decomposition avoids the intermediate `value * Num` overflow that occurs with
+// (Num*T)/Den. GCD reduction of Num and Den maximizes the safe input range.
+
+namespace internal {
+
+template <typename T>
+AU_DEVICE_FUNC constexpr T au_gcd(T a, T b) {
+    // Rely on truncation towards zero to make this work for negative numbers as well.
+    return b == T{0} ? a : au_gcd(b, a % b);
+}
+
+}  // namespace internal
+
+template <typename T, typename Num, typename Den, bool IsIntegral>
+struct ScaleByRationalImpl;
+
+template <typename T, typename Num, typename Den>
+struct OpInputImpl<ScaleByRational<T, Num, Den>> : stdx::type_identity<T> {};
+template <typename T, typename Num, typename Den>
+struct OpOutputImpl<ScaleByRational<T, Num, Den>> : stdx::type_identity<T> {};
+
+template <typename T, typename Num, typename Den>
+struct ScaleByRationalImpl<T, Num, Den, true> {
+    static AU_DEVICE_FUNC constexpr T apply_to(T value) {
+        using RealT = RealPart<T>;
+        constexpr auto p_raw = get_value<RealT>(Num{});
+        constexpr auto q_raw = get_value<RealT>(Den{});
+        constexpr auto g = internal::au_gcd(p_raw, q_raw);
+
+        constexpr auto p = p_raw / g;
+        constexpr auto q = q_raw / g;
+
+        return static_cast<T>(p * (value / q) + (p * (value % q)) / q);
+    }
+};
+
+template <typename T, typename Num, typename Den>
+struct ScaleByRationalImpl<T, Num, Den, false> {
+    static AU_DEVICE_FUNC constexpr T apply_to(T value) {
+        return static_cast<T>(
+            value * get_value<RealPart<T>>(MagProductT<Num, MagInverseT<Den>>{})
+        );
+    }
+};
+
+template <typename T, typename Num, typename Den>
+struct ScaleByRational {
+    static AU_DEVICE_FUNC constexpr T apply_to(T value) {
+        return ScaleByRationalImpl<
+            T, Num, Den,
+            std::is_integral<RealPart<T>>::value
+        >::apply_to(value);
+    }
+};
+
 }  // namespace detail
 }  // namespace au

au/abstract_operations_test.cc

@@ -14,6 +14,8 @@
 
 #include "au/abstract_operations.hh"
 
+#include <limits>
+
 #include "au/testing.hh"
 #include "gmock/gmock.h"
 #include "gtest/gtest.h"
@@ -146,6 +148,34 @@ TEST(OpSequence, EliminatesRedundantOperations) {
         OpSequence<StaticCast<int, float>, MultiplyTypeBy<float, decltype(mag<2>())>>>();
 }
 
+////////////////////////////////////////////////////////////////////////////////////////////////////
+// `ScaleByRational` section:
+
+TEST(ScaleByRational, HasExpectedInputAndOutputTypes) {
+    StaticAssertTypeEq<OpInput<ScaleByRational<int64_t, decltype(mag<9>()), decltype(mag<10>())>>,
+                       int64_t>();
+    StaticAssertTypeEq<OpOutput<ScaleByRational<int64_t, decltype(mag<9>()), decltype(mag<10>())>>,
+                       int64_t>();
+}
+
+TEST(ScaleByRational, NoOverflowForNegativeInt64NearBoundary) {
+    using Num = decltype(mag<9>());
+    using Den = decltype(mag<10>());
+
+    EXPECT_EQ((ScaleByRational<int64_t, Num, Den>::apply_to(-2000000000000000000LL)),
+              -1800000000000000000LL);
+}
+
+TEST(ScaleByRational, NoOverflowForInt64Min) {
+    using Num = decltype(mag<9>());
+    using Den = decltype(mag<10>());
+
+    const int64_t result =
+        ScaleByRational<int64_t, Num, Den>::apply_to(std::numeric_limits<int64_t>::min());
+
+    EXPECT_EQ(result, -8301034833169298227LL);
+}
+
 }  // namespace
 }  // namespace detail
 }  // namespace au

au/apply_magnitude_test.cc

@@ -270,19 +270,17 @@ TEST(WouldOverflow, AlwaysFalseForIntegerDivide) {
 
 TEST(WouldOverflow, UsesNumeratorWhenApplyingRationalMagnitudeToIntegralType) {
     {
-        using ApplyTwoThirdsToI32 = ApplyMagnitudeT<int32_t, decltype(mag<2>() / mag<3>())>;
-
-        EXPECT_THAT(ApplyTwoThirdsToI32::would_overflow(2'147'483'647), IsTrue());
-        EXPECT_THAT(ApplyTwoThirdsToI32::would_overflow(1'073'741'824), IsTrue());
-
-        EXPECT_THAT(ApplyTwoThirdsToI32::would_overflow(1'073'741'823), IsFalse());
-        EXPECT_THAT(ApplyTwoThirdsToI32::would_overflow(1), IsFalse());
-        EXPECT_THAT(ApplyTwoThirdsToI32::would_overflow(0), IsFalse());
-        EXPECT_THAT(ApplyTwoThirdsToI32::would_overflow(-1), IsFalse());
-        EXPECT_THAT(ApplyTwoThirdsToI32::would_overflow(-1'073'741'824), IsFalse());
-
-        EXPECT_THAT(ApplyTwoThirdsToI32::would_overflow(-1'073'741'825), IsTrue());
-        EXPECT_THAT(ApplyTwoThirdsToI32::would_overflow(-2'147'483'648), IsTrue());
+        using ApplyFractionToI32 = ApplyMagnitudeT<int32_t, decltype(mag<46341>() / mag<46342>())>;
+
+        EXPECT_THAT(ApplyFractionToI32::would_overflow(2'147'483'647), IsTrue());
+        EXPECT_THAT(ApplyFractionToI32::would_overflow(46'341), IsTrue());
+        EXPECT_THAT(ApplyFractionToI32::would_overflow(46'340), IsFalse());
+        EXPECT_THAT(ApplyFractionToI32::would_overflow(1), IsFalse());
+        EXPECT_THAT(ApplyFractionToI32::would_overflow(0), IsFalse());
+        EXPECT_THAT(ApplyFractionToI32::would_overflow(-1), IsFalse());
+        EXPECT_THAT(ApplyFractionToI32::would_overflow(-46'340), IsFalse());
+        EXPECT_THAT(ApplyFractionToI32::would_overflow(-46'341), IsTrue());
+        EXPECT_THAT(ApplyFractionToI32::would_overflow(-2'147'483'648), IsTrue());
     }
 
     {
@@ -294,13 +292,21 @@ TEST(WouldOverflow, UsesNumeratorWhenApplyingRationalMagnitudeToIntegralType) {
 
         EXPECT_THAT(ApplyRoughlyOneThirdToU8::would_overflow(255), IsTrue());
         EXPECT_THAT(ApplyRoughlyOneThirdToU8::would_overflow(22), IsTrue());
-
         EXPECT_THAT(ApplyRoughlyOneThirdToU8::would_overflow(21), IsFalse());
         EXPECT_THAT(ApplyRoughlyOneThirdToU8::would_overflow(1), IsFalse());
         EXPECT_THAT(ApplyRoughlyOneThirdToU8::would_overflow(0), IsFalse());
     }
 }
 
+TEST(WouldOverflow, NegativeRational) {
+    {
+        using ApplyNegativeFraction = ApplyMagnitudeT<int32_t, decltype(-mag<2>() / mag<3>())>;
+
+        EXPECT_THAT(ApplyNegativeFraction::would_overflow(2'147'483'647), IsTrue());
+        EXPECT_THAT(ApplyNegativeFraction::would_overflow(-2'147'483'648), IsTrue());
+    }
+}
+
 TEST(WouldOverflow, UsesFullValueWhenApplyingRationalMagnitudeToFloatingPointType) {
     {
         using ApplyThreeHalvesToF = ApplyMagnitudeT<float, decltype(mag<3>() / mag<2>())>;

au/conversion_strategy.hh

@@ -81,7 +81,7 @@ template <typename T, typename Mag>
 struct ApplicationStrategyForImpl<T, Mag, MagKindHolder<MagKind::NONTRIVIAL_RATIONAL>>
     : std::conditional<
           std::is_integral<RealPart<T>>::value,
-          OpSequence<MultiplyTypeBy<T, Numerator<Mag>>, DivideTypeByInteger<T, Denominator<Mag>>>,
+          ScaleByRational<T, NumeratorT<Mag>, DenominatorT<Mag>>,
           MultiplyTypeBy<T, Mag>> {};
 
 //

au/conversion_strategy_test.cc

@@ -72,13 +72,12 @@ TEST(ConversionForRepsAndFactor, SameRepForPromotingTypeHasCastAtBeginningAndEnd
                                   ImplicitConversion<Promoted, T>>>();
 }
 
-TEST(ConversionForRepsAndFactor, ApplyingNontrivialRationalToIntegralTypeIsMultiplyThenDivide) {
+TEST(ConversionForRepsAndFactor, ApplyingNontrivialRationalToIntegralIsScaleByRational) {
     StaticAssertTypeEq<
         CastStrategyIndependentConversionForRepsAndFactor<uint64_t,
                                                           uint64_t,
                                                           decltype(mag<3>() / mag<4>())>,
-        OpSequence<MultiplyTypeBy<uint64_t, decltype(mag<3>())>,
-                   DivideTypeByInteger<uint64_t, decltype(mag<4>())>>>();
+        ScaleByRational<uint64_t, decltype(mag<3>()), decltype(mag<4>())>>();
 }
 
 TEST(ConversionForRepsAndFactor, ApplyingNontrivialRationalToFloatingPointIsSingleMultiply) {
@@ -90,13 +89,12 @@ TEST(ConversionForRepsAndFactor, ApplyingNontrivialRationalToFloatingPointIsSing
 }
 
 TEST(ConversionForRepsAndFactor,
-     ApplyingNontrivialRationalToComplexIntegralTypeIsMultiplyThenDivide) {
+     ApplyingNontrivialRationalToComplexIntegralTypeIsScaleByRational) {
     StaticAssertTypeEq<
         CastStrategyIndependentConversionForRepsAndFactor<std::complex<int>,
                                                           std::complex<int>,
                                                           decltype(mag<3>() / mag<4>())>,
-        OpSequence<MultiplyTypeBy<std::complex<int>, decltype(mag<3>())>,
-                   DivideTypeByInteger<std::complex<int>, decltype(mag<4>())>>>();
+        ScaleByRational<std::complex<int>, decltype(mag<3>()), decltype(mag<4>())>>();
 }
 
 TEST(ConversionForRepsAndFactor, WhenTargetIsPromotedTypeSkipFinalStaticCast) {

au/overflow_boundary.hh

@@ -735,6 +735,95 @@ template <typename T, typename M, typename Limits>
 struct MaxGoodImpl<DivideTypeByInteger<T, M>, Limits>
     : MaxGoodImplForDivideTypeByIntegerUsingRealPart<RealPart<T>, M, Limits> {};
 
+////////////////////////////////////////////////////////////////////////////////////////////////////
+// `ScaleByRational<T, Num, Den>` implementation.
+
+template <typename Num, typename Den>
+struct IsMagnifyingOrNegativeRational {
+    using M = MagProductT<Num, MagInverseT<Den>>;
+    using P = NumeratorT<M>;
+    using Q = DenominatorT<M>;
+    static constexpr auto P_RESULT = get_value_result<uintmax_t>(P{});
+    static constexpr auto Q_RESULT = get_value_result<uintmax_t>(Q{});
+    static constexpr bool value =
+        (P_RESULT.outcome != MagRepresentationOutcome::OK) ||
+        (Q_RESULT.outcome == MagRepresentationOutcome::OK &&
+         P_RESULT.value > Q_RESULT.value);
+};
+
+// Return true when the intermediate modular remainder Num * (value % Den)
+// overflows type T for some values within T's range.
+template <typename T, typename Num, typename Den>
+struct IsScaleByRationalIntermediateOverflowPossible {
+    using M = MagProductT<Num, MagInverseT<Den>>;
+    using P = NumeratorT<M>;
+    using Q = DenominatorT<M>;
+    static constexpr auto P_RESULT = get_value_result<T>(P{});
+    static constexpr auto Q_RESULT = get_value_result<T>(Q{});
+    static constexpr bool value =
+        (Q_RESULT.outcome != MagRepresentationOutcome::OK) ||
+        (Q_RESULT.value > T{1} && P_RESULT.outcome == MagRepresentationOutcome::OK &&
+         P_RESULT.value > std::numeric_limits<T>::max() / (Q_RESULT.value - T{1}));
+};
+
+template <typename T, typename Limits>
+struct ScaleByRationalAtLowerLimit {
+    static constexpr T value() { return LowerLimit<T, Limits>::value(); }
+};
+
+template <typename T, typename Limits>
+struct ScaleByRationalAtUpperLimit {
+    static constexpr T value() { return UpperLimit<T, Limits>::value(); }
+};
+
+template <typename T, typename Num, typename Den, typename Limits>
+struct ScaleByRationalNarrowMinBound {
+    using M = MagProductT<Num, MagInverseT<Den>>;
+    using P = NumeratorT<M>;
+    static constexpr T value() {
+        return std::numeric_limits<T>::lowest() / get_value<T>(P{});
+    }
+};
+
+template <typename T, typename Num, typename Den, typename Limits>
+struct ScaleByRationalNarrowMaxBound {
+    using M = MagProductT<Num, MagInverseT<Den>>;
+    using P = NumeratorT<M>;
+    static constexpr T value() {
+        return std::numeric_limits<T>::max() / get_value<T>(P{});
+    }
+};
+
+template <typename T, typename Num, typename Den, typename Limits>
+struct MinGoodImpl<ScaleByRational<T, Num, Den>, Limits> {
+    using RT = RealPart<T>;
+    using M = MagProductT<Num, MagInverseT<Den>>;
+    using P = NumeratorT<M>;
+    using Q = DenominatorT<M>;
+    using type = std::conditional_t<
+        IsMagnifyingOrNegativeRational<Num, Den>::value,
+        typename MinGoodImpl<OpSequence<MultiplyTypeBy<RT, P>, DivideTypeByInteger<RT, Q>>,
+                             Limits>::type,
+        std::conditional_t<IsScaleByRationalIntermediateOverflowPossible<RT, Num, Den>::value,
+                           ScaleByRationalNarrowMinBound<RT, Num, Den, Limits>,
+                           ScaleByRationalAtLowerLimit<RT, Limits>>>;
+};
+
+template <typename T, typename Num, typename Den, typename Limits>
+struct MaxGoodImpl<ScaleByRational<T, Num, Den>, Limits> {
+    using RT = RealPart<T>;
+    using M = MagProductT<Num, MagInverseT<Den>>;
+    using P = NumeratorT<M>;
+    using Q = DenominatorT<M>;
+    using type = std::conditional_t<
+        IsMagnifyingOrNegativeRational<Num, Den>::value,
+        typename MaxGoodImpl<OpSequence<MultiplyTypeBy<RT, P>, DivideTypeByInteger<RT, Q>>,
+                             Limits>::type,
+        std::conditional_t<IsScaleByRationalIntermediateOverflowPossible<RT, Num, Den>::value,
+                           ScaleByRationalNarrowMaxBound<RT, Num, Den, Limits>,
+                           ScaleByRationalAtUpperLimit<RT, Limits>>>;
+};
+
 ////////////////////////////////////////////////////////////////////////////////////////////////////
 // `OpSequence<Ops...>` implementation.
 

au/truncation_risk.hh

@@ -135,6 +135,10 @@ template <typename T, typename M>
 struct TruncationRiskForImpl<MultiplyTypeBy<T, M>>
     : TruncationRiskForMultiplyByAssumingScalar<RealPart<T>, M> {};
 
+template <typename T, typename Num, typename Den>
+struct TruncationRiskForImpl<ScaleByRational<T, Num, Den>>
+    : TruncationRiskForImpl<MultiplyTypeBy<T, MagProductT<Num, MagInverseT<Den>>>> {};
+
 ////////////////////////////////////////////////////////////////////////////////////////////////////
 // `DivideTypeByInteger<T, M>` section:
 
@@ -216,6 +220,19 @@ template <template <class> class Risk, typename T, typename M>
 struct UpdateRiskImpl<DivideTypeByInteger<T, M>, Risk<RealPart<T>>>
     : stdx::type_identity<Risk<RealPart<T>>> {};
 
+template <template <class> class Risk, typename T, typename Num, typename Den>
+struct UpdateRiskImpl<ScaleByRational<T, Num, Den>, Risk<RealPart<T>>>
+    : stdx::type_identity<Risk<RealPart<T>>> {};
+
+template <typename T, typename Num, typename Den, typename M>
+struct UpdateRiskImpl<ScaleByRational<T, Num, Den>, ValueTimesRatioIsNotInteger<RealPart<T>, M>> {
+    using ScaleMag = MagProductT<Num, MagInverseT<Den>>;
+    using type = typename std::conditional<IsRational<ScaleMag>::value,
+                                           ReduceValueTimesRatioIsNotInteger<RealPart<T>,
+                                                                             MagProductT<ScaleMag, M>>,
+                                           ValueIsNotZero<RealPart<T>>>::type;
+};
+
 template <typename T, typename M1, typename M2>
 struct UpdateRiskImpl<MultiplyTypeBy<T, M1>, ValueTimesRatioIsNotInteger<RealPart<T>, M2>>
     : std::conditional<IsRational<M1>::value,