mirror of
https://github.com/XRPLF/rippled.git
synced 2026-06-02 16:26:48 +00:00
fix: Improve upward rounding edge cases for Number::operator/= (#7328)
Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> Co-authored-by: Vito Tumas <5780819+Tapanito@users.noreply.github.com>
This commit is contained in:
Ed Hennis
@@ -2,12 +2,14 @@
|
||||
|
||||
#include <xrpl/beast/utility/instrumentation.h>
|
||||
|
||||
#include <array>
|
||||
#include <cstdint>
|
||||
#include <functional>
|
||||
#include <limits>
|
||||
#include <optional>
|
||||
#include <ostream>
|
||||
#include <set>
|
||||
#include <stdexcept>
|
||||
#include <string>
|
||||
#include <unordered_map>
|
||||
|
||||
@@ -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<std::uint64_t, kInt64Digits>
|
||||
buildPowersOfTen()
|
||||
{
|
||||
std::array<std::uint64_t, kInt64Digits> 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<std::uint64_t>::max() / 10)
|
||||
throw std::logic_error("Power of 10 table is too big");
|
||||
}
|
||||
result[exponent] = power;
|
||||
if (power < std::numeric_limits<std::uint64_t>::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<std::uint64_t, detail::kInt64Digits> 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<decltype(MinMantissa)>,
|
||||
Integral64 TMax = std::decay_t<decltype(MaxMantissa)>>
|
||||
Integral64 T = std::decay_t<decltype(MinMantissa)>>
|
||||
[[nodiscard]]
|
||||
std::pair<T, int>
|
||||
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 <auto MinMantissa, auto MaxMantissa, Integral64 T, Integral64 TMax>
|
||||
template <auto MinMantissa, auto MaxMantissa, Integral64 T>
|
||||
std::pair<T, int>
|
||||
Number::normalizeToRange() const
|
||||
{
|
||||
static_assert(std::is_same_v<T, std::uint64_t> || std::is_same_v<T, std::int64_t>);
|
||||
static_assert(std::is_same_v<T, TMax>);
|
||||
static_assert(std::is_same_v<T, std::decay_t<decltype(MinMantissa)>>);
|
||||
static_assert(std::is_same_v<T, std::decay_t<decltype(MaxMantissa)>>);
|
||||
auto constexpr kMIN = static_cast<T>(MinMantissa);
|
||||
auto constexpr kMAX = static_cast<T>(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);
|
||||
|
||||
|
||||
@@ -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 <class T>
|
||||
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<uint128_t>(
|
||||
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<unsigned long long>(
|
||||
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<unsigned long>(
|
||||
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<uint128_t>(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<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;
|
||||
// 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<internalrep>(zm);
|
||||
exponent_ = ze;
|
||||
XRPL_ASSERT_PARTS(isnormal(), "xrpl::Number::operator/=", "result is normalized");
|
||||
|
||||
@@ -6,11 +6,14 @@
|
||||
#include <xrpl/protocol/SystemParameters.h>
|
||||
#include <xrpl/protocol/XRPAmount.h>
|
||||
|
||||
// NOLINTNEXTLINE(misc-include-cleaner)
|
||||
#include <boost/multiprecision/cpp_dec_float.hpp>
|
||||
#include <boost/multiprecision/number.hpp>
|
||||
|
||||
#include <array>
|
||||
#include <cctype>
|
||||
#include <cstdint>
|
||||
#include <iomanip>
|
||||
#include <limits>
|
||||
#include <map>
|
||||
#include <sstream>
|
||||
@@ -40,6 +43,40 @@ 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)
|
||||
{
|
||||
if (n == 0)
|
||||
return 1;
|
||||
if (n == 1)
|
||||
return 10;
|
||||
|
||||
if (n > 1)
|
||||
{
|
||||
auto r = pow10<T>(n / 2);
|
||||
r *= r;
|
||||
if (n % 2 != 0)
|
||||
r *= 10;
|
||||
return r;
|
||||
}
|
||||
|
||||
// n < 0
|
||||
T p = 1;
|
||||
p /= pow10<T>(-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<std::int64_t>::max() / 10) + 1);
|
||||
BEAST_EXPECT(product.exponent() == 19);
|
||||
switch (scale)
|
||||
{
|
||||
case MantissaRange::MantissaScale::Large:
|
||||
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);
|
||||
break;
|
||||
|
||||
case MantissaRange::MantissaScale::LargeLegacy:
|
||||
BEAST_EXPECT(signedDifference < 0);
|
||||
BEAST_EXPECT(
|
||||
product.mantissa() ==
|
||||
(std::numeric_limits<std::int64_t>::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<std::int64_t>::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();
|
||||
}
|
||||
};
|
||||
|
||||
|
||||
Reference in New Issue
Block a user