diff --git a/src/libxrpl/basics/Number.cpp b/src/libxrpl/basics/Number.cpp index 11f5934b04..0a22d439ea 100644 --- a/src/libxrpl/basics/Number.cpp +++ b/src/libxrpl/basics/Number.cpp @@ -858,31 +858,51 @@ Number::operator/=(Number const& y) auto const& maxMantissa = range.max; auto const cuspRoundingFixEnabled = range.cuspRoundingFixEnabled; + // Cache power of 10 computations to not waste time on subsequent calls + auto getPower10 = [](int exponent) { + static std::unordered_map powersOf10; + if (!powersOf10.contains(exponent)) + { + uint128_t result = 1; + for (int i = 0; i < exponent; ++i) + result *= 10; + powersOf10[exponent] = result; + } + return powersOf10.at(exponent); + }; + // 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 + // uint128_t or the cast back to int64_t for small mantissas. + // + // For large mantissas: + // * log(2^128,10) ~ 38.5 + // * largeRange.log = 18, fits in 10^19 + // * The expanded numerator must fit in 10^38 + // * 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; + auto const factorExponent = small ? 17 : 19; + + uint128_t const f = getPower10(factorExponent); XRPL_ASSERT_PARTS(f >= minMantissa * 10, "Number::operator/=", "factor expected size"); // 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; - auto const numerator = uint128_t(nm) * f; auto zm = numerator / dmu; - auto ze = ne - de - (small ? 17 : 19); + auto ze = ne - de - factorExponent; bool zn = (ns * ds) < 0; if (!small) { + // zm may now hold up to 20 digits. The correctionFactor must be small enough to not cause + // overflow, but as large as possible otherwise + // * zm fits in 10^21 + // * correctionFactor can be up to 10^(38-21) = 10^17 safely + bool const useBigCorrection = + cuspRoundingFixEnabled == MantissaRange::CuspRoundingFix::Enabled; + auto const correctionExponent = useBigCorrection ? 17 : 3; + uint128_t const correctionFactor = getPower10(correctionExponent); + // Virtually multiply numerator by correctionFactor. Since that would // overflow in the existing uint128_t, we'll do that part separately. // The math for this would work for small mantissas, but we need to @@ -906,11 +926,12 @@ Number::operator/=(Number const& y) 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 + // divide by the correctionFactor by moving the exponent, so we don't lose the // integer value we just computed - ze -= 3; + ze -= correctionExponent; + + auto const correction = (remainder * correctionFactor) / dmu; + zm += correction; } } normalize(zn, zm, ze, minMantissa, maxMantissa, cuspRoundingFixEnabled); diff --git a/src/test/basics/Number_test.cpp b/src/test/basics/Number_test.cpp index 52af19ed36..9a80caa8e4 100644 --- a/src/test/basics/Number_test.cpp +++ b/src/test/basics/Number_test.cpp @@ -6,10 +6,12 @@ #include #include +#include #include #include #include +#include #include #include #include @@ -39,6 +41,30 @@ class Number_test : public beast::unit_test::Suite return out; } + using dec = boost::multiprecision::cpp_dec_float_50; + + template + static T + pow10(int n) + { + T p = 1; + if (n >= 0) + for (int i = 0; i < n; ++i) + p *= 10; + else + for (int i = 0; i < -n; ++i) + p /= 10; + return p; + } + + static std::string + fmt(dec const& v) + { + std::ostringstream os; + os << std::setprecision(40) << v; + return os.str(); + } + public: void testZero() @@ -1588,40 +1614,99 @@ public: void testUpwardRoundsDown() { - testcase << "upward rounding produces a value below exact at kMaxRep cusp"; + { + testcase << "upward rounding produces a value below exact at kMaxRep cusp"; - NumberMantissaScaleGuard const mg{MantissaRange::MantissaScale::Large}; - NumberRoundModeGuard const rg{Number::RoundingMode::Upward}; + NumberMantissaScaleGuard const mg{MantissaRange::MantissaScale::Large}; + 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; - // Public conversion operator: STAmount::operator Number() const. - Number const a = kAValue; - Number const b = kBValue; - Number const product = a * b; + // Public conversion operator: STAmount::operator Number() const. + 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"; - BEAST_EXPECT(signedDifference >= 0); - BEAST_EXPECT(product.mantissa() == (std::numeric_limits::max() / 10) + 1); - BEAST_EXPECT(product.exponent() == 19); + 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); + log.flush(); + } + + { + /* Companion to NumberUpwardWrongDirection_test (which targets + * `operator*=` Upward at the kMaxRep cusp on LargeLegacy), but for + * `operator/=` on the cusp-fix-ENABLED `Large` scale. + * + * Under `Large` (`CuspRoundingFix::Enabled`), `operator/=` with Upward + * rounding can return a value STRICTLY LESS than the exact quotient, + * violating Upward's directional invariant. + * + * Mechanism (fix-enabled path): + * 1. `operator/=` computes `numerator = nm * 10^19` and + * `zm = numerator / dm` (integer division, truncates remainder). + * 2. If `remainder != 0`, the correction block runs: + * zm *= 1000 + * correction = (remainder * 1000) / dm // also truncates + * zm += correction + * ze -= 3 + * The truncation in `correction` discards a sub-1/1000 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 "; + + NumberMantissaScaleGuard const scaleGuard{MantissaRange::MantissaScale::Large}; + 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"; + + // Upward invariant: stored >= exact. Bug: stored < exact. + BEAST_EXPECT(stored >= exact); + BEAST_EXPECT(diff < pow10(quotient.exponent())); + } } void