diff --git a/au/abstract_operations.hh b/au/abstract_operations.hh index 1129f65f..b7fc285e 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,63 @@ struct OpSequenceImpl { } }; +//////////////////////////////////////////////////////////////////////////////////////////////////// +// `ScaleByRational` 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 +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 +struct ScaleByRationalImpl; + +template +struct OpInputImpl> : stdx::type_identity {}; +template +struct OpOutputImpl> : stdx::type_identity {}; + +template +struct ScaleByRationalImpl { + static AU_DEVICE_FUNC constexpr T apply_to(T value) { + using RealT = RealPart; + constexpr auto p_raw = get_value(Num{}); + constexpr auto q_raw = get_value(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(p * (value / q) + (p * (value % q)) / q); + } +}; + +template +struct ScaleByRationalImpl { + static AU_DEVICE_FUNC constexpr T apply_to(T value) { + return static_cast( + value * get_value>(MagProductT>{}) + ); + } +}; + +template +struct ScaleByRational { + static AU_DEVICE_FUNC constexpr T apply_to(T value) { + return ScaleByRationalImpl< + T, Num, Den, + std::is_integral>::value + >::apply_to(value); + } +}; + } // namespace detail } // namespace au diff --git a/au/abstract_operations_test.cc b/au/abstract_operations_test.cc index c1c0e937..27f72dc9 100644 --- a/au/abstract_operations_test.cc +++ b/au/abstract_operations_test.cc @@ -14,6 +14,8 @@ #include "au/abstract_operations.hh" +#include + #include "au/testing.hh" #include "gmock/gmock.h" #include "gtest/gtest.h" @@ -146,6 +148,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_magnitude_test.cc b/au/apply_magnitude_test.cc index 1308ecd0..acbe234c 100644 --- a/au/apply_magnitude_test.cc +++ b/au/apply_magnitude_test.cc @@ -270,19 +270,17 @@ TEST(WouldOverflow, AlwaysFalseForIntegerDivide) { TEST(WouldOverflow, UsesNumeratorWhenApplyingRationalMagnitudeToIntegralType) { { - using ApplyTwoThirdsToI32 = ApplyMagnitudeT() / 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() / 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() / 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() / mag<2>())>; 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..8a862481 100644 --- a/au/overflow_boundary.hh +++ b/au/overflow_boundary.hh @@ -735,6 +735,95 @@ 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); +}; + +// Return true when the intermediate modular remainder Num * (value % Den) +// overflows type T for some values within T's range. +template +struct IsScaleByRationalIntermediateOverflowPossible { + 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 = + (Q_RESULT.outcome != MagRepresentationOutcome::OK) || + (Q_RESULT.value > T{1} && P_RESULT.outcome == MagRepresentationOutcome::OK && + P_RESULT.value > std::numeric_limits::max() / (Q_RESULT.value - T{1})); +}; + +template +struct ScaleByRationalAtLowerLimit { + static constexpr T value() { return LowerLimit::value(); } +}; + +template +struct ScaleByRationalAtUpperLimit { + static constexpr T value() { return UpperLimit::value(); } +}; + +template +struct ScaleByRationalNarrowMinBound { + using M = MagProductT>; + using P = NumeratorT; + static constexpr T value() { + return std::numeric_limits::lowest() / get_value(P{}); + } +}; + +template +struct ScaleByRationalNarrowMaxBound { + using M = MagProductT>; + using P = NumeratorT; + static constexpr T value() { + return std::numeric_limits::max() / get_value(P{}); + } +}; + +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, + std::conditional_t::value, + ScaleByRationalNarrowMinBound, + 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, + std::conditional_t::value, + ScaleByRationalNarrowMaxBound, + ScaleByRationalAtUpperLimit>>; +}; + //////////////////////////////////////////////////////////////////////////////////////////////////// // `OpSequence` implementation. diff --git a/au/truncation_risk.hh b/au/truncation_risk.hh index 3e153b67..bd04d02f 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,19 @@ template