diff --git a/include/xrpl/basics/Number.h b/include/xrpl/basics/Number.h index 3aeb4b5bcd..93bef82a8c 100644 --- a/include/xrpl/basics/Number.h +++ b/include/xrpl/basics/Number.h @@ -2,12 +2,14 @@ #include +#include #include #include #include #include #include #include +#include #include #include @@ -40,6 +42,47 @@ isPowerOfTen(T value) return logTen(value).has_value(); } +namespace detail { + +/** Builds a table of the powers of 10 + * + * This function is marked consteval, so it can only be run in + * a constexpr context. This assures that it is and can only be run at + * compile time. Doing it at runtime would be pretty wasteful and + * inefficient. + */ +constexpr std::size_t kInt64Digits = 20; +consteval std::array +buildPowersOfTen() +{ + std::array result{}; + + std::uint64_t power = 1; + std::size_t exponent = 0; + // end the loop early so it doesn't overflow; + for (; exponent < result.size() - 1; ++exponent, power *= 10) + { + result[exponent] = power; + if (power > std::numeric_limits::max() / 10) + throw std::logic_error("Power of 10 table is too big"); + } + result[exponent] = power; + if (power < std::numeric_limits::max() / 10) + throw std::logic_error("Power of 10 table is not big enough for the uint64_t type"); + + return result; +} + +} // namespace detail + +constexpr std::array kPowerOfTen = detail::buildPowersOfTen(); + +static_assert(kPowerOfTen[0] == 1); +static_assert(kPowerOfTen[1] == 10); +static_assert(kPowerOfTen[10] == 10'000'000'000); +static_assert( + isPowerOfTen(kPowerOfTen.back()) && *logTen(kPowerOfTen.back()) == detail::kInt64Digits - 1); + /** MantissaRange defines a range for the mantissa of a normalized Number. * * The mantissa is in the range [min, max], where @@ -76,6 +119,7 @@ isPowerOfTen(T value) struct MantissaRange final { using rep = std::uint64_t; + enum class MantissaScale { Small, // LargeLegacy can be removed when fixCleanup3_2_0 is retired @@ -89,19 +133,15 @@ struct MantissaRange final Enabled = true, }; - explicit constexpr MantissaRange(MantissaScale scale) - : min(getMin(scale)) - , cuspRoundingFixEnabled(isCuspFixEnabled(scale)) - , log(logTen(min).value_or(-1)) - , scale(scale) + explicit constexpr MantissaRange(MantissaScale sc) : scale(sc) { } - rep min; - rep max{(min * 10) - 1}; - CuspRoundingFix cuspRoundingFixEnabled; - int log; - MantissaScale scale; + MantissaScale const scale; + int const log{getExponent(scale)}; + rep const min{getMin(scale, log)}; + rep const max{(min * 10) - 1}; + CuspRoundingFix const cuspRoundingFixEnabled{isCuspFixEnabled(scale)}; static MantissaRange const& getMantissaRange(MantissaScale scale); @@ -110,23 +150,35 @@ struct MantissaRange final getAllScales(); private: - static constexpr rep - getMin(MantissaScale scale) + static constexpr int + getExponent(MantissaScale scale) { switch (scale) { case MantissaScale::Small: - return 1'000'000'000'000'000ULL; + return 15; case MantissaScale::LargeLegacy: case MantissaScale::Large: - return 1'000'000'000'000'000'000ULL; + return 18; + // LCOV_EXCL_START default: // If called in a constexpr context, this throw assures that the build fails if an // invalid scale is used. - throw std::runtime_error("Unknown mantissa scale"); // LCOV_EXCL_LINE + throw std::runtime_error("Unknown mantissa scale"); + // LCOV_EXCL_STOP } } + // Keep this function for future use with different ways to compute + // the ranges. + static constexpr rep + getMin(MantissaScale scale, int exponent) + { + if (exponent < 0 || exponent >= kPowerOfTen.size()) + throw std::runtime_error("Invalid exponent"); // LCOV_EXCL_LINE + return kPowerOfTen[exponent]; + } + static constexpr CuspRoundingFix isCuspFixEnabled(MantissaScale scale) { @@ -477,8 +529,7 @@ public: template < auto MinMantissa, auto MaxMantissa, - Integral64 T = std::decay_t, - Integral64 TMax = std::decay_t> + Integral64 T = std::decay_t> [[nodiscard]] std::pair normalizeToRange() const; @@ -519,7 +570,8 @@ private: int& exponent, MantissaRange::rep const& minMantissa, MantissaRange::rep const& maxMantissa, - MantissaRange::CuspRoundingFix cuspRoundingFixEnabled); + MantissaRange::CuspRoundingFix cuspRoundingFixEnabled, + bool dropped); [[nodiscard]] bool isnormal() const noexcept; @@ -725,16 +777,18 @@ Number::isnormal() const noexcept kMinExponent <= exponent_ && exponent_ <= kMaxExponent); } -template +template std::pair Number::normalizeToRange() const { static_assert(std::is_same_v || std::is_same_v); - static_assert(std::is_same_v); + static_assert(std::is_same_v>); + static_assert(std::is_same_v>); auto constexpr kMIN = static_cast(MinMantissa); auto constexpr kMAX = static_cast(MaxMantissa); static_assert(kMIN > 0); static_assert(kMIN % 10 == 0); + static_assert(isPowerOfTen(kMIN)); static_assert(kMAX % 10 == 9); static_assert((kMAX + 1) / 10 == kMIN); diff --git a/src/libxrpl/basics/Number.cpp b/src/libxrpl/basics/Number.cpp index 11f5934b04..275d82d8c9 100644 --- a/src/libxrpl/basics/Number.cpp +++ b/src/libxrpl/basics/Number.cpp @@ -178,6 +178,10 @@ public: setPositive() noexcept; void setNegative() noexcept; + // Should only be called by doNormalize, and then only for division + // operations with remainders. + void + setDropped() noexcept; [[nodiscard]] bool isNegative() const noexcept; @@ -250,6 +254,12 @@ Number::Guard::setNegative() noexcept sbit_ = 1; } +inline void +Number::Guard::setDropped() noexcept +{ + xbit_ = 1; +} + inline bool Number::Guard::isNegative() const noexcept { @@ -398,7 +408,7 @@ Number::Guard::doRoundUp( // _don't_ increment the mantissa. Instead, divide and round recursively. It should // be impossible to recurse more than once, because once the mantissa is divided by // 10, it will be _well_ under maxMantissa and kMaxRep, so adding 1 will have no - // change of bringing it back over. + // chance of bringing it back over. doDropDigit(mantissa, exponent); XRPL_ASSERT_PARTS( safeToIncrement(mantissa), @@ -512,8 +522,6 @@ Number::one() return Number{false, range.min, -range.log, Number::Unchecked{}}; } -// Use the member names in this static function for now so the diff is cleaner -// TODO: Rename the function parameters to get rid of the "_" suffix template void doNormalize( @@ -522,7 +530,8 @@ doNormalize( int& exponent, MantissaRange::rep const& minMantissa, MantissaRange::rep const& maxMantissa, - MantissaRange::CuspRoundingFix cuspRoundingFixEnabled) + MantissaRange::CuspRoundingFix cuspRoundingFixEnabled, + bool dropped) { static constexpr auto kMinExponent = Number::kMinExponent; static constexpr auto kMaxExponent = Number::kMaxExponent; @@ -547,6 +556,8 @@ doNormalize( Guard g; if (negative) g.setNegative(); + if (dropped) + g.setDropped(); while (m > maxMantissa) { if (exponent >= kMaxExponent) @@ -611,7 +622,12 @@ Number::normalize( internalrep const& maxMantissa, MantissaRange::CuspRoundingFix cuspRoundingFixEnabled) { - doNormalize(negative, mantissa, exponent, minMantissa, maxMantissa, cuspRoundingFixEnabled); + // Not used by every compiler version, and thus not necessarily + // counted by coverage build + // LCOV_EXCL_START + doNormalize( + negative, mantissa, exponent, minMantissa, maxMantissa, cuspRoundingFixEnabled, false); + // LCOV_EXCL_STOP } template <> @@ -624,7 +640,12 @@ Number::normalize( internalrep const& maxMantissa, MantissaRange::CuspRoundingFix cuspRoundingFixEnabled) { - doNormalize(negative, mantissa, exponent, minMantissa, maxMantissa, cuspRoundingFixEnabled); + // Not used by every compiler version, and thus not necessarily + // counted by coverage build + // LCOV_EXCL_START + doNormalize( + negative, mantissa, exponent, minMantissa, maxMantissa, cuspRoundingFixEnabled, false); + // LCOV_EXCL_STOP } template <> @@ -637,7 +658,8 @@ Number::normalize( internalrep const& maxMantissa, MantissaRange::CuspRoundingFix cuspRoundingFixEnabled) { - doNormalize(negative, mantissa, exponent, minMantissa, maxMantissa, cuspRoundingFixEnabled); + doNormalize( + negative, mantissa, exponent, minMantissa, maxMantissa, cuspRoundingFixEnabled, false); } void @@ -838,7 +860,9 @@ Number::operator/=(Number const& y) return *this; // n* = numerator // d* = denominator - // *p = negative (positive?) + // z* = result (quotient) + // *p = negative (p for positive, even though the value means not + // positive?) // *s = sign // *m = mantissa // *e = exponent @@ -850,71 +874,155 @@ Number::operator/=(Number const& y) bool const dp = y.negative_; int const ds = (dp ? -1 : 1); - auto dm = y.mantissa_; - auto de = y.exponent_; + // Create the denominator as 128-bit unsigned, since that's what we + // need to work with. + uint128_t const dm = static_cast(y.mantissa_); + auto const de = y.exponent_; auto const& range = kRange.get(); auto const& minMantissa = range.min; auto const& maxMantissa = range.max; auto const cuspRoundingFixEnabled = range.cuspRoundingFixEnabled; - // Shift by 10^17 gives greatest precision while not overflowing - // uint128_t or the cast back to int64_t - // TODO: Can/should this be made bigger for largeRange? - // log(2^128,10) ~ 38.5 - // largeRange.log = 18, fits in 10^19 - // f can be up to 10^(38-19) = 10^19 safely - bool const small = Number::getMantissaScale() == MantissaRange::MantissaScale::Small; - uint128_t const f = small ? 100'000'000'000'000'000 : 10'000'000'000'000'000'000ULL; - XRPL_ASSERT_PARTS(f >= minMantissa * 10, "Number::operator/=", "factor expected size"); + // Division operates on two large integers (16-digit for small + // mantissas, 19-digit for large) using integer math. If the values + // were just divided directly, the result would be only ever be one + // digit or zero - not very useful. + // e.g. 9'876'543'210'987'654 / 1'234'567'890'123'456 = 8 + // 1'234'567'890'123'456 / 9'876'543'210'987'654 = 0 + // Introduce a power-of-ten multiplication factor for the numerator + // which will ensure the result has a meaningful number of digits. + // + // Consider numbers with a 2-digit mantissa: + // * Assume both numbers have an exponent of 0, using "ToNearest" rounding + // * 23 / 67 = 0 + // * Use a factor of 10^4 + // * 230'000 / 67 = 3432 with an exponent of -4 + // * The normalized result will be 34, exponent -2, or 0.34 + // + // The most extreme results are 10/99 and 99/10 + // * 100'000 / 99 = 1'010e-4 = 10e-2 or 0.10 + // * 990'000 / 10 = 99'000e-4 = 99e-1 or 9.9 + // + // Note that the computations give 2 or 3 digits after the + // decimal point to determine which way to round for most scenarios. + // + // For small mantissas (where the MantissaRange.log == 15), shifting by 10^17 gives sufficient + // precision while not overflowing uint128_t or the cast back to int64_t. (This is legacy + // behavior, which must not be changed.) + // + // For large mantissas (where the MantissaRange.log == 18), a shift by 10^20 would be optimal + // for most scenarios. However, larger mantissa values would overflow 2^128. + // + // * log(2^128,10) ~ 38.5 + // * largeRange.log = 18, fits in 10^19 + // * The expanded numerator must fit in 10^38 + // * f not be more than 10^(38-19) = 10^19 safely + // + // So, we do the division into stages: + // + // Stage 1: Use the same factor of 10^17, for the initial division. This + // will frequently not result in a whole number quotient. + // + // Stage 2: If there is a remainder from the first step, repeat the + // process with a "correction" factor of 10^5. Shift the + // result of Stage 1 over by 5 places, and add the second result to it. + // This is equivalent to if we had used an initial factor of 10^22, + // a couple digits more than we actually need. + // + // Stage 3: If there is still a remainder, and the CuspRoundingFix + // is enabled, pass a flag indicating such to doNormalize. The Guard + // in doNormalize will treat that flag as if non-zero digits had + // been dropped from the mantissa when shrinking it into range. + // This is only relevant when rounding away from zero (Upward for + // positive numbers, Downward for negative), or if the "regular" + // remainder is exactly 0.5 for "ToNearest". This will give the + // rounding the most accurate result possible, as if infinite + // precision was used in the initial calculation. - // unsigned denominator - auto const dmu = static_cast(dm); - // correctionFactor can be anything between 10 and f, depending on how much - // extra precision we want to only use for rounding with the - // largeRange. Three digits seems like plenty, and is more than - // the smallRange uses. - uint128_t const correctionFactor = 1'000; + // Stage 1: Do the initial division with a factor of 10^17. + auto constexpr factorExponent = 17; + + uint128_t constexpr f = kPowerOfTen[factorExponent]; auto const numerator = uint128_t(nm) * f; - auto zm = numerator / dmu; - auto ze = ne - de - (small ? 17 : 19); - bool zn = (ns * ds) < 0; - if (!small) + auto zm = numerator / dm; + auto ze = ne - de - factorExponent; + bool zp = (ns * ds) < 0; + // dropped is used in the same way as Guard::xbit_. In the case of + // division, it indicates if there's any remainder left over after + // we have been as precise as reasonable. If there is, it would be as + // if we were using infinite precision math, and a non-zero digit + // had been shifted off the end of the result when normalizing. + bool dropped = false; + + if (range.scale != MantissaRange::MantissaScale::Small) { - // Virtually multiply numerator by correctionFactor. Since that would - // overflow in the existing uint128_t, we'll do that part separately. + // Stage 2 + // + // If there is a remainder, treat it as a secondary numerator. + // Multiply by correctionFactor separately from stage 1. // The math for this would work for small mantissas, but we need to - // preserve existing behavior. + // preserve legacy behavior. // // Consider: - // ((numerator * correctionFactor) / dmu) / correctionFactor - // = ((numerator / dmu) * correctionFactor) / correctionFactor) + // ((numerator * correctionFactor) / dm) / correctionFactor + // = ((numerator / dm) * correctionFactor) / correctionFactor) // // But that assumes infinite precision. With integer math, this is // equivalent to // - // = ((numerator / dmu * correctionFactor) - // + ((numerator % dmu) * correctionFactor) / dmu) / correctionFactor + // = ((numerator / dm * correctionFactor) + // + ((numerator % dm) * correctionFactor) / dm) / correctionFactor + // = ((zm * correctionFactor) + // + (remainder * correctionFactor) / dm) / correctionFactor // - // We have already set `mantissa_ = numerator / dmu`. Now we - // compute `remainder = numerator % dmu`, and if it is - // nonzero, we do the rest of the arithmetic. If it's zero, we can skip - // it. - auto const remainder = (numerator % dmu); + // The trick is that multiplication by correctionFactor is done on the mantissa, but + // division by correctionFactor is done by modifying the exponent, so no precision is lost + // until we normalize. + // + // If remainder is zero, we can skip this stage entirely because + // the first stage gave an exact answer. + auto constexpr correctionExponent = 5; + uint128_t constexpr correctionFactor = kPowerOfTen[correctionExponent]; + static_assert(factorExponent + correctionExponent == 22); + + auto const remainder = (numerator % dm); if (remainder != 0) { - zm *= correctionFactor; - auto const correction = remainder * correctionFactor / dmu; - zm += correction; - // divide by 1000 by moving the exponent, so we don't lose the - // integer value we just computed - ze -= 3; + auto const partialNumerator = remainder * correctionFactor; + auto const correction = partialNumerator / dm; + + // If the correction is zero, we do not have to make any + // modifications to z*, because it will not have any + // effect on the final result. (We'd be adding a bunch of + // zeros to the end of zm that would just be removed in + // normalize.) However, if that is the case, then Stage 3 is + // even more important for accuracy. + if (correction != 0) + { + zm *= correctionFactor; + // divide by the correctionFactor by moving the exponent, so we don't lose the + // integer value we just computed + ze -= correctionExponent; + + zm += correction; + } + + // Stage 3: If there's still anything left, and the cusp + // rounding fix is enabled, flag if there is still + // a remainder from stage 2. + bool const useTrailingRemainder = + cuspRoundingFixEnabled == MantissaRange::CuspRoundingFix::Enabled; + if (useTrailingRemainder) + { + dropped = partialNumerator % dm != 0; + } } } - normalize(zn, zm, ze, minMantissa, maxMantissa, cuspRoundingFixEnabled); - negative_ = zn; + doNormalize(zp, zm, ze, minMantissa, maxMantissa, cuspRoundingFixEnabled, dropped); + negative_ = zp; mantissa_ = static_cast(zm); exponent_ = ze; XRPL_ASSERT_PARTS(isnormal(), "xrpl::Number::operator/=", "result is normalized"); diff --git a/src/test/basics/Number_test.cpp b/src/test/basics/Number_test.cpp index 26060c70e9..81019970ad 100644 --- a/src/test/basics/Number_test.cpp +++ b/src/test/basics/Number_test.cpp @@ -6,11 +6,14 @@ #include #include +// NOLINTNEXTLINE(misc-include-cleaner) +#include #include #include #include #include +#include #include #include #include @@ -40,6 +43,40 @@ class Number_test : public beast::unit_test::Suite return out; } + using dec = boost::multiprecision::cpp_dec_float_50; + + template + static T + pow10(int n) + { + if (n == 0) + return 1; + if (n == 1) + return 10; + + if (n > 1) + { + auto r = pow10(n / 2); + r *= r; + if (n % 2 != 0) + r *= 10; + return r; + } + + // n < 0 + T p = 1; + p /= pow10(-n); + return p; + } + + static std::string + fmt(dec const& v) + { + std::ostringstream os; + os << std::setprecision(40) << v; + return os.str(); + } + public: void testZero() @@ -1589,39 +1626,249 @@ public: void testUpwardRoundsDown() { - testcase << "upward rounding produces a value below exact at kMaxRep cusp"; + auto const scale = Number::getMantissaScale(); + { + testcase << "upward rounding produces a value below exact at kMaxRep cusp " + << to_string(scale); - NumberMantissaScaleGuard const mg{MantissaRange::MantissaScale::Large}; - NumberRoundModeGuard const rg{Number::RoundingMode::Upward}; + NumberRoundModeGuard const rg{Number::RoundingMode::Upward}; - constexpr std::int64_t kAValue = 1'000'000'000'000'049'863LL; - constexpr std::int64_t kBValue = 9'223'372'036'854'315'903LL; + constexpr std::int64_t kAValue = 1'000'000'000'000'049'863LL; + constexpr std::int64_t kBValue = 9'223'372'036'854'315'903LL; - Number const a = kAValue; - Number const b = kBValue; - Number const product = a * b; + Number const a = kAValue; + Number const b = kBValue; + Number const product = a * b; - // Exact reference in BigInt. - BigInt const exactProduct = BigInt(kAValue) * BigInt(kBValue); + // Exact reference in BigInt. + BigInt const exactProduct = BigInt(kAValue) * BigInt(kBValue); - // What Number actually stored. - BigInt storedValue = BigInt(product.mantissa()); - for (int i = 0; i < product.exponent(); ++i) - storedValue *= 10; + // What Number actually stored. + BigInt storedValue = BigInt(product.mantissa()); + for (int i = 0; i < product.exponent(); ++i) + storedValue *= 10; - BigInt const signedDifference = storedValue - exactProduct; + BigInt const signedDifference = storedValue - exactProduct; - log << "\n" - << " a = " << fmt(BigInt(kAValue)) << "\n" - << " b = " << fmt(BigInt(kBValue)) << "\n" - << " exact a*b = " << fmt(exactProduct) << "\n" - << " stored = " << fmt(storedValue) << "\n" - << " stored - exact = " << fmt(signedDifference) << "\n" - << " upward = " << (signedDifference >= 0 ? "held" : "VIOLATED") << "\n"; + log << "\n" + << " a = " << fmt(BigInt(kAValue)) << "\n" + << " b = " << fmt(BigInt(kBValue)) << "\n" + << " exact a*b = " << fmt(exactProduct) << "\n" + << " stored = " << fmt(storedValue) << "\n" + << " stored - exact = " << fmt(signedDifference) << "\n" + << " upward = " << (signedDifference >= 0 ? "held" : "VIOLATED") << "\n" + << " stored.mantissa = " << product.mantissa() << "\n" + << " stored.exponent = " << product.exponent() << "\n"; + log.flush(); - BEAST_EXPECT(signedDifference >= 0); - BEAST_EXPECT(product.mantissa() == (std::numeric_limits::max() / 10) + 1); - BEAST_EXPECT(product.exponent() == 19); + switch (scale) + { + case MantissaRange::MantissaScale::Large: + BEAST_EXPECT(signedDifference >= 0); + BEAST_EXPECT(signedDifference < pow10(product.exponent())); + BEAST_EXPECT( + product.mantissa() == (std::numeric_limits::max() / 10) + 1); + BEAST_EXPECT(product.exponent() == 19); + break; + + case MantissaRange::MantissaScale::LargeLegacy: + BEAST_EXPECT(signedDifference < 0); + BEAST_EXPECT( + product.mantissa() == + (std::numeric_limits::max() / 100) * 100); + BEAST_EXPECT(product.exponent() == 18); + break; + + case MantissaRange::MantissaScale::Small: + // The seemingly weird rounding here is because + // a & b are both normalized, and both round up when + // being converted to Number, so you're really + // getting 1_000_000_000_000_050 * 9_223_372_036_854_316 + BEAST_EXPECT(signedDifference >= 0); + BEAST_EXPECT( + product.mantissa() == + (std::numeric_limits::max() / 1000) + 3); + BEAST_EXPECT(product.exponent() == 21); + break; + } + } + + { + /* Companion regression for the kMaxRep cusp behavior, but for + * `operator/=` on the cusp-fix-ENABLED `Large` scale. + * + * Before the dropped-remainder fix, `operator/=` with Upward + * rounding could return a value STRICTLY LESS than the exact quotient, + * violating Upward's directional invariant. + * + * Mechanism (fix-enabled path): + * 1. `operator/=` computes `numerator = nm * 10^17` and + * `zm = numerator / dm` (integer division, truncates remainder). + * 2. If `remainder != 0`, the correction block runs: + * zm *= 100000 + * correction = (remainder * 100000) / dm // also truncates + * zm += correction + * ze -= 5 + * The truncation in `correction` discards a sub-1/100000 residual. + * 3. `normalize`'s shift loop reduces zm to fit, but the discarded + * residual is BELOW the Guard's visibility, so the Guard sees fraction = 0. + * 4. Under Upward + positive, `round()` returns -1 (no round-up), and + * the algorithm returns the truncated zm + */ + testcase << "operator/= Upward on Large returns value < truth " << to_string(scale); + + NumberRoundModeGuard const roundGuard{Number::RoundingMode::Upward}; + + constexpr std::int64_t aValue = 2LL; + constexpr std::int64_t bValue = 1'000'000'000'000'000'007LL; + // bValue = 10^18 + 7 (prime, in [minMantissa, kMaxRep]). + + Number const a{aValue, 0}; + Number const b{bValue, 0}; + Number const quotient = a / b; + + dec const exact = dec(aValue) / dec(bValue); + dec const stored = dec(quotient.mantissa()) * pow10(quotient.exponent()); + dec const diff = stored - exact; + + log << "\n" + << " a = " << aValue << "\n" + << " b = " << bValue << "\n" + << " exact a/b = " << fmt(exact) << "\n" + << " stored a/b = " << fmt(stored) << "\n" + << " stored - exact = " << fmt(diff) + << " (negative => Upward gave value BELOW truth)\n" + << " quotient.mantissa = " << quotient.mantissa() << "\n" + << " quotient.exponent = " << quotient.exponent() << "\n"; + log.flush(); + + // Upward invariant: stored >= exact. Bug: stored < exact. + switch (scale) + { + case MantissaRange::MantissaScale::Large: + BEAST_EXPECT(stored >= exact); + BEAST_EXPECT(diff < pow10(quotient.exponent())); + break; + + case MantissaRange::MantissaScale::LargeLegacy: + BEAST_EXPECT(stored < exact); + BEAST_EXPECT(diff >= -pow10(quotient.exponent())); + break; + + case MantissaRange::MantissaScale::Small: + // Small mantissa doesn't have the correction for + // dropped remainders + BEAST_EXPECT(stored < exact); + break; + } + } + { + /* Companion test case for Upward positive operator/=: Downward negative + */ + testcase << "operator/= Downward on Large returns value < truth " << to_string(scale); + + NumberRoundModeGuard const roundGuard{Number::RoundingMode::Downward}; + + constexpr std::int64_t aValue = -2LL; + constexpr std::int64_t bValue = 1'000'000'000'000'000'007LL; + // bValue = 10^18 + 7 (prime, in [minMantissa, kMaxRep]). + + Number const a{aValue, 0}; + Number const b{bValue, 0}; + Number const quotient = a / b; + + dec const exact = dec(aValue) / dec(bValue); + dec const stored = dec(quotient.mantissa()) * pow10(quotient.exponent()); + dec const diff = stored - exact; + + log << "\n" + << " a = " << aValue << "\n" + << " b = " << bValue << "\n" + << " exact a/b = " << fmt(exact) << "\n" + << " stored a/b = " << fmt(stored) << "\n" + << " stored - exact = " << fmt(diff) + << " (positive => Downward gave value ABOVE truth)\n" + << " quotient.mantissa = " << quotient.mantissa() << "\n" + << " quotient.exponent = " << quotient.exponent() << "\n"; + log.flush(); + + // invariant: stored <= exact. Bug: stored > exact. + switch (scale) + { + case MantissaRange::MantissaScale::Large: + BEAST_EXPECT(stored <= exact); + BEAST_EXPECT(diff > -pow10(quotient.exponent())); + break; + + case MantissaRange::MantissaScale::LargeLegacy: + BEAST_EXPECT(stored > exact); + BEAST_EXPECT(diff <= pow10(quotient.exponent())); + break; + + case MantissaRange::MantissaScale::Small: + // Small mantissa doesn't have the correction for + // dropped remainders + BEAST_EXPECT(stored < exact); + break; + } + } + { + /* Companion test case for Upward positive operator/=: ToNearest + * + * With ToNearest, if the dropped digits are exactly "5", then the mantissa will be + * rounded to even. The numbers below result in a value where the unrounded mantissa + * ends in an even digit, and "infinite precision" would drop + * "500000000000000000145...", but doNormalize only sees "5". Without the rounding fix, + * doNormalize rounds down to the even value. With the rounding fix, doNormalize knows + * there are more digits beyond "5", and so rounds _up_ to the odd value. + */ + testcase << "operator/= ToNearest on Large returns value < truth " << to_string(scale); + + NumberRoundModeGuard const roundGuard{Number::RoundingMode::ToNearest}; + + constexpr std::int64_t aValue = 1'269'917'268'816'087'809LL; + constexpr std::int64_t bValue = 3'458'525'013'821'685'511LL; + // bValue = 10^18 + 7 (prime, in [minMantissa, kMaxRep]). + + Number const a{aValue, 0}; + Number const b{bValue, 0}; + Number const quotient = a / b; + + dec const exact = dec(aValue) / dec(bValue); + dec const stored = dec(quotient.mantissa()) * pow10(quotient.exponent()); + dec const diff = stored - exact; + + log << "\n" + << " a = " << aValue << "\n" + << " b = " << bValue << "\n" + << " exact a/b = " << fmt(exact) << "\n" + << " stored a/b = " << fmt(stored) << "\n" + << " stored - exact = " << fmt(diff) + << " (negative => ToNearest gave value BELOW truth)\n" + << " quotient.mantissa = " << quotient.mantissa() << "\n" + << " quotient.exponent = " << quotient.exponent() << "\n"; + log.flush(); + + // invariant: stored >= exact. Bug: stored < exact. + switch (scale) + { + case MantissaRange::MantissaScale::Large: + BEAST_EXPECT(stored >= exact); + BEAST_EXPECT(diff < pow10(quotient.exponent())); + break; + + case MantissaRange::MantissaScale::LargeLegacy: + BEAST_EXPECT(stored < exact); + BEAST_EXPECT(diff >= -pow10(quotient.exponent())); + break; + + case MantissaRange::MantissaScale::Small: + // Small mantissa doesn't have the correction for + // dropped remainders + BEAST_EXPECT(stored < exact); + break; + } + } } void @@ -1651,9 +1898,9 @@ public: testTruncate(); testRounding(); testInt64(); + + testUpwardRoundsDown(); } - // This test sets its own number range - testUpwardRoundsDown(); } };