Make Number::operator/= significantly more accurate

- Prevents extreme dust rounding from getting lost, especially when
  rounding away from zero. (Upward for positive, downward for negative.)
This commit is contained in:
Ed Hennis
2026-05-23 18:45:36 +01:00
parent 4ab886bcbc
commit 48b1716e6f
2 changed files with 149 additions and 43 deletions

View File

@@ -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<int, uint128_t> 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<uint128_t>(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);

View File

@@ -6,10 +6,12 @@
#include <xrpl/protocol/SystemParameters.h>
#include <xrpl/protocol/XRPAmount.h>
#include <boost/multiprecision/cpp_dec_float.hpp>
#include <boost/multiprecision/number.hpp>
#include <array>
#include <cstdint>
#include <iomanip>
#include <limits>
#include <map>
#include <sstream>
@@ -39,6 +41,30 @@ class Number_test : public beast::unit_test::Suite
return out;
}
using dec = boost::multiprecision::cpp_dec_float_50;
template <class T = dec>
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<std::int64_t>::max() / 10) + 1);
BEAST_EXPECT(product.exponent() == 19);
BEAST_EXPECT(signedDifference >= 0);
BEAST_EXPECT(signedDifference < pow10<BigInt>(product.exponent()));
BEAST_EXPECT(product.mantissa() == (std::numeric_limits<std::int64_t>::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