From 26adcc8b08e793d21140c889265fb469b2b7a998 Mon Sep 17 00:00:00 2001 From: Jordan Montgomery Date: Fri, 15 May 2026 13:24:29 -0700 Subject: [PATCH 1/4] Add a ScaleByRational struct to mitigate problems with overflow in large numbers Instead of applying (P*T)/Q, instead use modular decomposition to compute the rational without overflow. We can extend this for some values of P,Q beyond the usual limit by using GCD(P,Q) instead. --- au/abstract_operations.hh | 43 ++++++++++++++++ au/abstract_operations_test.cc | 29 +++++++++++ ...ply_rational_magnitude_to_integral_test.cc | 16 +++--- au/conversion_strategy.hh | 2 +- au/conversion_strategy_test.cc | 10 ++-- au/overflow_boundary.hh | 51 +++++++++++++++++++ au/truncation_risk.hh | 15 ++++++ 7 files changed, 151 insertions(+), 15 deletions(-) diff --git a/au/abstract_operations.hh b/au/abstract_operations.hh index 1129f65f..10e30f02 100644 --- a/au/abstract_operations.hh +++ b/au/abstract_operations.hh @@ -78,6 +78,13 @@ struct OpSequenceImpl; template using OpSequence = FlattenAs; +// +// `ScaleByRational` represents an operation that scales a value of type `T` by the +// rational number `Num/Den`, taking care to avoid overflow if possible. +// +template +struct ScaleByRational; + //////////////////////////////////////////////////////////////////////////////////////////////////// // IMPLEMENTATION DETAILS (`abstract_operations.hh`): //////////////////////////////////////////////////////////////////////////////////////////////////// @@ -188,5 +195,41 @@ struct OpSequenceImpl { } }; +//////////////////////////////////////////////////////////////////////////////////////////////////// +// `ScaleByRational` implementation. +// +// Modular decomposition avoids the intermediate `value * numerator` overflow that afflicts the +// naive multiply-then-divide approach. GCD reduction of p and q maximises the safe input range. + +template +AU_DEVICE_FUNC constexpr T au_gcd(T a, T b) { + return b == T{0} ? a : au_gcd(b, a % b); +} + +template +struct OpInputImpl> : stdx::type_identity {}; +template +struct OpOutputImpl> : stdx::type_identity {}; + +template +struct ScaleByRational { + + static AU_DEVICE_FUNC constexpr T apply_to(T value) { + if constexpr (std::is_integral>::value) { + constexpr auto p_raw = get_value>(Num{}); + constexpr auto q_raw = get_value>(Den{}); + constexpr auto g = au_gcd(p_raw, q_raw); + constexpr auto p = p_raw / g; + constexpr auto q = q_raw / g; + static_assert(q <= 1 || p <= std::numeric_limits>::max() / (q - 1), + "p*(q-1) overflows T; unit ratio cannot be applied to this integral type"); + return static_cast(p * (value / q) + p * (value % q) / q); + } else { + return static_cast( + value * get_value>(MagProductT>{})); + } + } +}; + } // namespace detail } // namespace au diff --git a/au/abstract_operations_test.cc b/au/abstract_operations_test.cc index c1c0e937..e03635cf 100644 --- a/au/abstract_operations_test.cc +++ b/au/abstract_operations_test.cc @@ -17,6 +17,7 @@ #include "au/testing.hh" #include "gmock/gmock.h" #include "gtest/gtest.h" +#include using ::testing::StaticAssertTypeEq; @@ -146,6 +147,34 @@ TEST(OpSequence, EliminatesRedundantOperations) { OpSequence, MultiplyTypeBy())>>>(); } +//////////////////////////////////////////////////////////////////////////////////////////////////// +// `ScaleByRational` section: + +TEST(ScaleByRational, HasExpectedInputAndOutputTypes) { + StaticAssertTypeEq()), decltype(mag<10>())>>, + int64_t>(); + StaticAssertTypeEq()), decltype(mag<10>())>>, + int64_t>(); +} + +TEST(ScaleByRational, NoOverflowForNegativeInt64NearBoundary) { + using Num = decltype(mag<9>()); + using Den = decltype(mag<10>()); + + EXPECT_EQ((ScaleByRational::apply_to(-2000000000000000000LL)), + -1800000000000000000LL); +} + +TEST(ScaleByRational, NoOverflowForInt64Min) { + using Num = decltype(mag<9>()); + using Den = decltype(mag<10>()); + + const int64_t result = + ScaleByRational::apply_to(std::numeric_limits::min()); + + EXPECT_EQ(result, -8301034833169298227LL); +} + } // namespace } // namespace detail } // namespace au diff --git a/au/apply_rational_magnitude_to_integral_test.cc b/au/apply_rational_magnitude_to_integral_test.cc index cab193c5..cb51475f 100644 --- a/au/apply_rational_magnitude_to_integral_test.cc +++ b/au/apply_rational_magnitude_to_integral_test.cc @@ -205,7 +205,7 @@ TEST(MaxNonOverflowingValue, IsMaxTDividedByNWhenTIsNotPromotableAndDenomOverflo {IsPromotable::NO, NumFitsInPromotedType::YES, DenFitsInPromotedType::NO}, huge_denom, max_int); - EXPECT_THAT(max_int, Eq(std::numeric_limits::max() / 3)); + EXPECT_THAT(max_int, Eq(std::numeric_limits::max())); } { @@ -214,7 +214,7 @@ TEST(MaxNonOverflowingValue, IsMaxTDividedByNWhenTIsNotPromotableAndDenomOverflo {IsPromotable::NO, NumFitsInPromotedType::YES, DenFitsInPromotedType::NO}, huge_denom, max_u64); - EXPECT_THAT(max_u64, Eq(std::numeric_limits::max() / 3)); + EXPECT_THAT(max_u64, Eq(std::numeric_limits::max())); } } @@ -239,7 +239,7 @@ TEST(MaxNonOverflowingValue, ASSERT_THAT((std::is_same, int32_t>::value), IsTrue()) << "This test will fail on architectures where uint16_t is not promoted to `int32_t`"; - EXPECT_THAT(max_u16, Eq(2'147)); + EXPECT_THAT(max_u16, Eq(std::numeric_limits::max())); } } @@ -262,7 +262,7 @@ TEST(MaxNonOverflowingValue, {IsPromotable::YES, NumFitsInPromotedType::YES, DenFitsInPromotedType::YES}, mag<1'000'000>() / pow<6>(mag<11>()), max_u16); - EXPECT_THAT(max_u16, Eq(2'147)); + EXPECT_THAT(max_u16, Eq(std::numeric_limits::max())); } } @@ -339,7 +339,7 @@ TEST(MinNonOverflowingValue, IsMinTDividedByNWhenTIsNotPromotableAndDenomOverflo {IsPromotable::NO, NumFitsInPromotedType::YES, DenFitsInPromotedType::NO}, huge_denom, min_int); - EXPECT_THAT(min_int, Eq(std::numeric_limits::lowest() / 3)); + EXPECT_THAT(min_int, Eq(std::numeric_limits::lowest())); } { @@ -348,7 +348,7 @@ TEST(MinNonOverflowingValue, IsMinTDividedByNWhenTIsNotPromotableAndDenomOverflo {IsPromotable::NO, NumFitsInPromotedType::YES, DenFitsInPromotedType::NO}, huge_denom, min_i64); - EXPECT_THAT(min_i64, Eq(std::numeric_limits::lowest() / 3)); + EXPECT_THAT(min_i64, Eq(std::numeric_limits::lowest())); } } @@ -372,7 +372,7 @@ TEST(MinNonOverflowingValue, ASSERT_THAT((std::is_same, int32_t>::value), IsTrue()) << "This test will fail on architectures where `int16_t` is not promoted to `int32_t`"; - EXPECT_THAT(min_i16, Eq(-2'147)); + EXPECT_THAT(min_i16, Eq(std::numeric_limits::lowest())); } } @@ -394,7 +394,7 @@ TEST(MinNonOverflowingValue, {IsPromotable::YES, NumFitsInPromotedType::YES, DenFitsInPromotedType::YES}, mag<1'000'000>() / pow<6>(mag<11>()), min_i16); - EXPECT_THAT(min_i16, Eq(-2'147)); + EXPECT_THAT(min_i16, Eq(std::numeric_limits::lowest())); } } diff --git a/au/conversion_strategy.hh b/au/conversion_strategy.hh index 00f27684..cac9404a 100644 --- a/au/conversion_strategy.hh +++ b/au/conversion_strategy.hh @@ -81,7 +81,7 @@ template struct ApplicationStrategyForImpl> : std::conditional< std::is_integral>::value, - OpSequence>, DivideTypeByInteger>>, + ScaleByRational, DenominatorT>, MultiplyTypeBy> {}; // diff --git a/au/conversion_strategy_test.cc b/au/conversion_strategy_test.cc index 33ebb2e4..9a4fbc6f 100644 --- a/au/conversion_strategy_test.cc +++ b/au/conversion_strategy_test.cc @@ -72,13 +72,12 @@ TEST(ConversionForRepsAndFactor, SameRepForPromotingTypeHasCastAtBeginningAndEnd ImplicitConversion>>(); } -TEST(ConversionForRepsAndFactor, ApplyingNontrivialRationalToIntegralTypeIsMultiplyThenDivide) { +TEST(ConversionForRepsAndFactor, ApplyingNontrivialRationalToIntegralIsScaleByRational) { StaticAssertTypeEq< CastStrategyIndependentConversionForRepsAndFactor() / mag<4>())>, - OpSequence())>, - DivideTypeByInteger())>>>(); + ScaleByRational()), decltype(mag<4>())>>(); } TEST(ConversionForRepsAndFactor, ApplyingNontrivialRationalToFloatingPointIsSingleMultiply) { @@ -90,13 +89,12 @@ TEST(ConversionForRepsAndFactor, ApplyingNontrivialRationalToFloatingPointIsSing } TEST(ConversionForRepsAndFactor, - ApplyingNontrivialRationalToComplexIntegralTypeIsMultiplyThenDivide) { + ApplyingNontrivialRationalToComplexIntegralTypeIsScaleByRational) { StaticAssertTypeEq< CastStrategyIndependentConversionForRepsAndFactor, std::complex, decltype(mag<3>() / mag<4>())>, - OpSequence, decltype(mag<3>())>, - DivideTypeByInteger, decltype(mag<4>())>>>(); + ScaleByRational, decltype(mag<3>()), decltype(mag<4>())>>(); } TEST(ConversionForRepsAndFactor, WhenTargetIsPromotedTypeSkipFinalStaticCast) { diff --git a/au/overflow_boundary.hh b/au/overflow_boundary.hh index e38aa209..b99ed2c1 100644 --- a/au/overflow_boundary.hh +++ b/au/overflow_boundary.hh @@ -735,6 +735,57 @@ template struct MaxGoodImpl, Limits> : MaxGoodImplForDivideTypeByIntegerUsingRealPart, M, Limits> {}; +//////////////////////////////////////////////////////////////////////////////////////////////////// +// `ScaleByRational` implementation. +template +struct IsMagnifyingOrNegativeRational { + using M = MagProductT>; + using P = NumeratorT; + using Q = DenominatorT; + static constexpr auto P_RESULT = get_value_result(P{}); + static constexpr auto Q_RESULT = get_value_result(Q{}); + static constexpr bool value = + (P_RESULT.outcome != MagRepresentationOutcome::OK) || + (Q_RESULT.outcome == MagRepresentationOutcome::OK && + P_RESULT.value > Q_RESULT.value); +}; + +template +struct ScaleByRationalAtLowerLimit { + static constexpr T value() { return LowerLimit::value(); } +}; + +template +struct ScaleByRationalAtUpperLimit { + static constexpr T value() { return UpperLimit::value(); } +}; + +template +struct MinGoodImpl, Limits> { + using RT = RealPart; + using M = MagProductT>; + using P = NumeratorT; + using Q = DenominatorT; + using type = std::conditional_t< + IsMagnifyingOrNegativeRational::value, + typename MinGoodImpl, DivideTypeByInteger>, + Limits>::type, + ScaleByRationalAtLowerLimit>; +}; + +template +struct MaxGoodImpl, Limits> { + using RT = RealPart; + using M = MagProductT>; + using P = NumeratorT; + using Q = DenominatorT; + using type = std::conditional_t< + IsMagnifyingOrNegativeRational::value, + typename MaxGoodImpl, DivideTypeByInteger>, + Limits>::type, + ScaleByRationalAtUpperLimit>; +}; + //////////////////////////////////////////////////////////////////////////////////////////////////// // `OpSequence` implementation. diff --git a/au/truncation_risk.hh b/au/truncation_risk.hh index 3e153b67..54dde003 100644 --- a/au/truncation_risk.hh +++ b/au/truncation_risk.hh @@ -135,6 +135,10 @@ template struct TruncationRiskForImpl> : TruncationRiskForMultiplyByAssumingScalar, M> {}; +template +struct TruncationRiskForImpl> + : TruncationRiskForImpl>>> {}; + //////////////////////////////////////////////////////////////////////////////////////////////////// // `DivideTypeByInteger` section: @@ -216,6 +220,17 @@ template