Experimental: Scale addition operands up to preserve accuracy

This commit is contained in:
Ed Hennis
2026-06-01 21:48:57 -04:00
parent 261508a0ec
commit 35bee87909
3 changed files with 262 additions and 83 deletions

View File

@@ -51,37 +51,43 @@ namespace detail {
* 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>
constexpr std::size_t kUint64Digits = 20;
constexpr std::size_t kUint128Digits = 39;
template <typename T, std::size_t digits>
consteval std::array<T, digits>
buildPowersOfTen()
{
std::array<std::uint64_t, kInt64Digits> result{};
std::array<T, digits> result{};
std::uint64_t power = 1;
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)
if (power > std::numeric_limits<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");
if (power < std::numeric_limits<T>::max() / 10)
throw std::logic_error("Power of 10 table is not big enough for the given type");
return result;
}
} // namespace detail
constexpr std::array<std::uint64_t, detail::kInt64Digits> kPowerOfTen = detail::buildPowersOfTen();
template <typename T = std::uint64_t, std::size_t digits = detail::kUint64Digits>
constexpr std::array<T, digits> kPowerOfTenImpl = detail::buildPowersOfTen<T, digits>();
constexpr auto kPowerOfTen = kPowerOfTenImpl<std::uint64_t, detail::kUint64Digits>;
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);
isPowerOfTen(kPowerOfTen.back()) && *logTen(kPowerOfTen.back()) == detail::kUint64Digits - 1);
/** MantissaRange defines a range for the mantissa of a normalized Number.
*

View File

@@ -720,26 +720,54 @@ Number::operator+=(Number const& y)
uint128_t ym = y.mantissa_;
auto ye = y.exponent_;
Guard g;
auto const& range = kRange.get();
// Bring the exponents of both values into agreement, so the mantissas are on the same scale
// and can be added directly together
// expandM / expandE: First try to expand the mantissa and bring the exponent down
// shringM / shrinkE: Then shrink the mantissa and bring the exponent up, if necessary
auto const adjust = [&g, &range](
uint128_t& expandM, int& expandE, uint128_t& shrinkM, int& shrinkE) {
constexpr uint128_t kSafeLimit = kPowerOfTenImpl<uint128_t, detail::kUint128Digits>[37];
if (range.cuspRoundingFixEnabled == MantissaRange::CuspRoundingFix::Enabled)
{
while (shrinkE < expandE && shrinkM % 10 == 0)
{
g.doDropDigit(shrinkM, shrinkE);
}
// We've got 128 bits of mantissa to work with here. Don't throw away data unless we
// have to
while (shrinkE < expandE && expandE > kMinExponent && expandM < kSafeLimit)
{
expandM *= 10;
--expandE;
}
}
while (shrinkE < expandE)
{
g.doDropDigit(shrinkM, shrinkE);
}
};
if (xe < ye)
{
if (xn)
g.setNegative();
do
{
g.doDropDigit(xm, xe);
} while (xe < ye);
adjust(ym, ye, xm, xe);
}
else if (xe > ye)
{
if (yn)
g.setNegative();
do
{
g.doDropDigit(ym, ye);
} while (xe > ye);
adjust(xm, xe, ym, ye);
}
auto const& range = kRange.get();
auto const& minMantissa = range.min;
auto const& maxMantissa = range.max;
auto const cuspRoundingFixEnabled = range.cuspRoundingFixEnabled;
@@ -747,9 +775,19 @@ Number::operator+=(Number const& y)
if (xn == yn)
{
xm += ym;
if (xm > maxMantissa || xm > kMaxRep)
if (range.cuspRoundingFixEnabled == MantissaRange::CuspRoundingFix::Enabled)
{
g.doDropDigit(xm, xe);
while (xm > maxMantissa || xm > kMaxRep)
{
g.doDropDigit(xm, xe);
}
}
else
{
if (xm > maxMantissa || xm > kMaxRep)
{
g.doDropDigit(xm, xe);
}
}
g.doRoundUp(
xn,
@@ -779,6 +817,25 @@ Number::operator+=(Number const& y)
--xe;
}
g.doRoundDown(xn, xm, xe, minMantissa);
if (range.cuspRoundingFixEnabled == MantissaRange::CuspRoundingFix::Enabled)
{
// make a new guard
Guard g;
if (xn)
g.setNegative();
while (xm > maxMantissa || xm > kMaxRep)
{
g.doDropDigit(xm, xe);
}
g.doRoundUp(
xn,
xm,
xe,
minMantissa,
maxMantissa,
cuspRoundingFixEnabled,
"Number::addition overflow");
}
}
negative_ = xn;

View File

@@ -43,6 +43,15 @@ class Number_test : public beast::unit_test::Suite
return out;
}
static BigInt
toBigInt(Number const& n)
{
BigInt v = n.mantissa();
for (int i = 0; i < n.exponent(); ++i)
v *= 10;
return v;
}
using dec = boost::multiprecision::cpp_dec_float_50;
template <class T = dec>
@@ -169,28 +178,34 @@ public:
auto const scale = Number::getMantissaScale();
testcase << "test_add " << to_string(scale);
using Case = std::tuple<Number, Number, Number>;
using Case = std::tuple<Number, Number, Number, int>;
auto const cSmall = std::to_array<Case>(
{{Number{1'000'000'000'000'000, -15},
Number{6'555'555'555'555'555, -29},
Number{1'000'000'000'000'066, -15}},
Number{1'000'000'000'000'066, -15},
__LINE__},
{Number{-1'000'000'000'000'000, -15},
Number{-6'555'555'555'555'555, -29},
Number{-1'000'000'000'000'066, -15}},
Number{-1'000'000'000'000'066, -15},
__LINE__},
{Number{-1'000'000'000'000'000, -15},
Number{6'555'555'555'555'555, -29},
Number{-9'999'999'999'999'344, -16}},
Number{-9'999'999'999'999'344, -16},
__LINE__},
{Number{-6'555'555'555'555'555, -29},
Number{1'000'000'000'000'000, -15},
Number{9'999'999'999'999'344, -16}},
{Number{}, Number{5}, Number{5}},
{Number{5}, Number{}, Number{5}},
Number{9'999'999'999'999'344, -16},
__LINE__},
{Number{}, Number{5}, Number{5}, __LINE__},
{Number{5}, Number{}, Number{5}, __LINE__},
{Number{5'555'555'555'555'555, -32768},
Number{-5'555'555'555'555'554, -32768},
Number{0}},
Number{0},
__LINE__},
{Number{-9'999'999'999'999'999, -31},
Number{1'000'000'000'000'000, -15},
Number{9'999'999'999'999'990, -16}}});
Number{9'999'999'999'999'990, -16},
__LINE__}});
auto const cLarge = std::to_array<Case>(
// Note that items with extremely large mantissas need to be
// calculated, because otherwise they overflow uint64. Items from C
@@ -198,45 +213,57 @@ public:
{
{Number{1'000'000'000'000'000, -15},
Number{6'555'555'555'555'555, -29},
Number{1'000'000'000'000'065'556, -18}},
Number{1'000'000'000'000'065'556, -18},
__LINE__},
{Number{-1'000'000'000'000'000, -15},
Number{-6'555'555'555'555'555, -29},
Number{-1'000'000'000'000'065'556, -18}},
Number{-1'000'000'000'000'065'556, -18},
__LINE__},
{Number{-1'000'000'000'000'000, -15},
Number{6'555'555'555'555'555, -29},
Number{true, 9'999'999'999'999'344'444ULL, -19, Number::Normalized{}}},
Number{true, 9'999'999'999'999'344'444ULL, -19, Number::Normalized{}},
__LINE__},
{Number{-6'555'555'555'555'555, -29},
Number{1'000'000'000'000'000, -15},
Number{false, 9'999'999'999'999'344'444ULL, -19, Number::Normalized{}}},
{Number{}, Number{5}, Number{5}},
{Number{5}, Number{}, Number{5}},
Number{false, 9'999'999'999'999'344'444ULL, -19, Number::Normalized{}},
__LINE__},
{Number{}, Number{5}, Number{5}, __LINE__},
{Number{5}, Number{}, Number{5}, __LINE__},
{Number{5'555'555'555'555'555'000, -32768},
Number{-5'555'555'555'555'554'000, -32768},
Number{0}},
Number{0},
__LINE__},
{Number{-9'999'999'999'999'999, -31},
Number{1'000'000'000'000'000, -15},
Number{9'999'999'999'999'990, -16}},
Number{9'999'999'999'999'990, -16},
__LINE__},
// Items from cSmall expanded for the larger mantissa
{Number{1'000'000'000'000'000'000, -18},
Number{6'555'555'555'555'555'555, -35},
Number{1'000'000'000'000'000'066, -18}},
Number{1'000'000'000'000'000'066, -18},
__LINE__},
{Number{-1'000'000'000'000'000'000, -18},
Number{-6'555'555'555'555'555'555, -35},
Number{-1'000'000'000'000'000'066, -18}},
Number{-1'000'000'000'000'000'066, -18},
__LINE__},
{Number{-1'000'000'000'000'000'000, -18},
Number{6'555'555'555'555'555'555, -35},
Number{true, 9'999'999'999'999'999'344ULL, -19, Number::Normalized{}}},
Number{true, 9'999'999'999'999'999'344ULL, -19, Number::Normalized{}},
__LINE__},
{Number{-6'555'555'555'555'555'555, -35},
Number{1'000'000'000'000'000'000, -18},
Number{false, 9'999'999'999'999'999'344ULL, -19, Number::Normalized{}}},
{Number{}, Number{5}, Number{5}},
Number{false, 9'999'999'999'999'999'344ULL, -19, Number::Normalized{}},
__LINE__},
{Number{}, Number{5}, Number{5}, __LINE__},
{Number{5'555'555'555'555'555'555, -32768},
Number{-5'555'555'555'555'555'554, -32768},
Number{0}},
Number{0},
__LINE__},
{Number{true, 9'999'999'999'999'999'999ULL, -37, Number::Normalized{}},
Number{1'000'000'000'000'000'000, -18},
Number{false, 9'999'999'999'999'999'990ULL, -19, Number::Normalized{}}},
{Number{Number::kMaxRep - 1}, Number{1, 0}, Number{Number::kMaxRep}},
Number{false, 9'999'999'999'999'999'990ULL, -19, Number::Normalized{}},
__LINE__},
{Number{Number::kMaxRep - 1}, Number{1, 0}, Number{Number::kMaxRep}, __LINE__},
// Test extremes
{
// Each Number operand rounds up, so the actual mantissa is
@@ -244,6 +271,7 @@ public:
Number{false, 9'999'999'999'999'999'999ULL, 0, Number::Normalized{}},
Number{false, 9'999'999'999'999'999'999ULL, 0, Number::Normalized{}},
Number{2, 19},
__LINE__,
},
{
// Does not round. Mantissas are going to be > kMaxRep, so if
@@ -254,21 +282,25 @@ public:
Number{false, 9'999'999'999'999'999'990ULL, 0, Number::Normalized{}},
Number{false, 9'999'999'999'999'999'990ULL, 0, Number::Normalized{}},
Number{false, 1'999'999'999'999'999'998ULL, 1, Number::Normalized{}},
__LINE__,
},
});
auto const cLargeLegacy = std::to_array<Case>({
{Number{Number::kMaxRep}, Number{6, -1}, Number{Number::kMaxRep / 10, 1}},
{Number{Number::kMaxRep}, Number{6, -1}, Number{Number::kMaxRep / 10, 1}, __LINE__},
});
auto const cLargeCorrected = std::to_array<Case>({
{Number{Number::kMaxRep}, Number{6, -1}, Number{(Number::kMaxRep / 10) + 1, 1}},
{Number{Number::kMaxRep},
Number{6, -1},
Number{(Number::kMaxRep / 10) + 1, 1},
__LINE__},
});
auto test = [this](auto const& c) {
for (auto const& [x, y, z] : c)
for (auto const& [x, y, z, line] : c)
{
auto const result = x + y;
std::stringstream ss;
ss << x << " + " << y << " = " << result << ". Expected: " << z;
BEAST_EXPECTS(result == z, ss.str());
expect(result == z, ss.str(), __FILE__, line);
}
};
if (scale == MantissaRange::MantissaScale::Small)
@@ -308,21 +340,28 @@ public:
auto const scale = Number::getMantissaScale();
testcase << "test_sub " << to_string(scale);
using Case = std::tuple<Number, Number, Number>;
using Case = std::tuple<Number, Number, Number, int>;
auto const cSmall = std::to_array<Case>(
{{Number{1'000'000'000'000'000, -15},
Number{6'555'555'555'555'555, -29},
Number{9'999'999'999'999'344, -16}},
Number{9'999'999'999'999'344, -16},
__LINE__},
{Number{6'555'555'555'555'555, -29},
Number{1'000'000'000'000'000, -15},
Number{-9'999'999'999'999'344, -16}},
{Number{1'000'000'000'000'000, -15}, Number{1'000'000'000'000'000, -15}, Number{0}},
Number{-9'999'999'999'999'344, -16},
__LINE__},
{Number{1'000'000'000'000'000, -15},
Number{1'000'000'000'000'000, -15},
Number{0},
__LINE__},
{Number{1'000'000'000'000'000, -15},
Number{1'000'000'000'000'001, -15},
Number{-1'000'000'000'000'000, -30}},
Number{-1'000'000'000'000'000, -30},
__LINE__},
{Number{1'000'000'000'000'001, -15},
Number{1'000'000'000'000'000, -15},
Number{1'000'000'000'000'000, -30}}});
Number{1'000'000'000'000'000, -30},
__LINE__}});
auto const cLarge = std::to_array<Case>(
// Note that items with extremely large mantissas need to be
// calculated, because otherwise they overflow uint64. Items from C
@@ -330,49 +369,63 @@ public:
{
{Number{1'000'000'000'000'000, -15},
Number{6'555'555'555'555'555, -29},
Number{false, 9'999'999'999'999'344'444ULL, -19, Number::Normalized{}}},
Number{false, 9'999'999'999'999'344'444ULL, -19, Number::Normalized{}},
__LINE__},
{Number{6'555'555'555'555'555, -29},
Number{1'000'000'000'000'000, -15},
Number{true, 9'999'999'999'999'344'444ULL, -19, Number::Normalized{}}},
{Number{1'000'000'000'000'000, -15}, Number{1'000'000'000'000'000, -15}, Number{0}},
Number{true, 9'999'999'999'999'344'444ULL, -19, Number::Normalized{}},
__LINE__},
{Number{1'000'000'000'000'000, -15},
Number{1'000'000'000'000'000, -15},
Number{0},
__LINE__},
{Number{1'000'000'000'000'000, -15},
Number{1'000'000'000'000'001, -15},
Number{-1'000'000'000'000'000, -30}},
Number{-1'000'000'000'000'000, -30},
__LINE__},
{Number{1'000'000'000'000'001, -15},
Number{1'000'000'000'000'000, -15},
Number{1'000'000'000'000'000, -30}},
Number{1'000'000'000'000'000, -30},
__LINE__},
// Items from cSmall expanded for the larger mantissa
{Number{1'000'000'000'000'000'000, -18},
Number{6'555'555'555'555'555'555, -32},
Number{false, 9'999'999'999'999'344'444ULL, -19, Number::Normalized{}}},
Number{false, 9'999'999'999'999'344'444ULL, -19, Number::Normalized{}},
__LINE__},
{Number{6'555'555'555'555'555'555, -32},
Number{1'000'000'000'000'000'000, -18},
Number{true, 9'999'999'999'999'344'444ULL, -19, Number::Normalized{}}},
Number{true, 9'999'999'999'999'344'444ULL, -19, Number::Normalized{}},
__LINE__},
{Number{1'000'000'000'000'000'000, -18},
Number{1'000'000'000'000'000'000, -18},
Number{0}},
Number{0},
__LINE__},
{Number{1'000'000'000'000'000'000, -18},
Number{1'000'000'000'000'000'001, -18},
Number{-1'000'000'000'000'000'000, -36}},
Number{-1'000'000'000'000'000'000, -36},
__LINE__},
{Number{1'000'000'000'000'000'001, -18},
Number{1'000'000'000'000'000'000, -18},
Number{1'000'000'000'000'000'000, -36}},
{Number{Number::kMaxRep}, Number{6, -1}, Number{Number::kMaxRep - 1}},
Number{1'000'000'000'000'000'000, -36},
__LINE__},
{Number{Number::kMaxRep}, Number{6, -1}, Number{Number::kMaxRep - 1}, __LINE__},
{Number{false, Number::kMaxRep + 1, 0, Number::Normalized{}},
Number{1, 0},
Number{(Number::kMaxRep / 10) + 1, 1}},
Number{(Number::kMaxRep / 10) + 1, 1},
__LINE__},
{Number{false, Number::kMaxRep + 1, 0, Number::Normalized{}},
Number{3, 0},
Number{Number::kMaxRep}},
{power(2, 63), Number{3, 0}, Number{Number::kMaxRep}},
Number{Number::kMaxRep},
__LINE__},
{power(2, 63), Number{3, 0}, Number{Number::kMaxRep}, __LINE__},
});
auto test = [this](auto const& c) {
for (auto const& [x, y, z] : c)
for (auto const& [x, y, z, line] : c)
{
auto const result = x - y;
std::stringstream ss;
ss << x << " - " << y << " = " << result << ". Expected: " << z;
BEAST_EXPECTS(result == z, ss.str());
expect(result == z, ss.str(), __FILE__, line);
}
};
if (scale == MantissaRange::MantissaScale::Small)
@@ -1644,9 +1697,7 @@ public:
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;
BigInt storedValue = toBigInt(product);
BigInt const signedDifference = storedValue - exactProduct;
@@ -1900,18 +1951,83 @@ public:
<< " Downward = " << downward << "\n\n";
log.flush();
// Upward round UP
BEAST_EXPECT(upward == above);
switch (scale)
{
case MantissaRange::MantissaScale::Small:
// With the small mantissa, everything rounds up
// ToNearest rounds UP when the DOWN neighbor is strictly closer
BEAST_EXPECT(toNearest != above);
BEAST_EXPECT(toNearest == below);
// Upward round UP
BEAST_EXPECT(upward > above);
// Downward undershoots: it returns a value below `below`
BEAST_EXPECT(downward == below);
// ToNearest rounds UP when the DOWN neighbor is strictly closer
BEAST_EXPECT(toNearest > above);
BEAST_EXPECT(toNearest == below);
// Both should have given the same answer, but they differ
BEAST_EXPECT(toNearest == downward);
// Downward undershoots: it returns a value below `below`
BEAST_EXPECT(downward < below);
// Both should have given the same answer, but they differ
BEAST_EXPECT(toNearest > downward);
break;
case MantissaRange::MantissaScale::LargeLegacy:
// Upward round UP
BEAST_EXPECT(upward == above);
// ToNearest rounds UP when the DOWN neighbor is strictly closer
BEAST_EXPECT(toNearest == above);
BEAST_EXPECT(toNearest > below);
// Downward undershoots: it returns a value below `below`
BEAST_EXPECT(downward < below);
// Both should have given the same answer, but they differ
BEAST_EXPECT(toNearest > downward);
break;
default:
// Upward round UP
BEAST_EXPECT(upward == above);
// ToNearest rounds UP when the DOWN neighbor is strictly closer
BEAST_EXPECT(toNearest != above);
BEAST_EXPECT(toNearest == below);
// Downward undershoots: it returns a value below `below`
BEAST_EXPECT(downward == below);
// Both should have given the same answer, but they differ
BEAST_EXPECT(toNearest == downward);
}
}
{
testcase << "operator+ TowardsZero rounds away from zero " << to_string(scale);
Number const a{1LL, 20};
Number const b{-1'000'000'000'000'000'001LL};
BEAST_EXPECT(toBigInt(a) == BigInt{"100000000000000000000"});
if (scale != MantissaRange::MantissaScale::Small)
BEAST_EXPECT(toBigInt(b) == BigInt{"-1000000000000000001"});
else
BEAST_EXPECT(toBigInt(b) == BigInt{"-1000000000000000000"});
Number sum;
{
NumberRoundModeGuard const roundGuard{Number::RoundingMode::TowardsZero};
sum = a + b;
}
BigInt const exact = toBigInt(a) + toBigInt(b);
BigInt const stored = toBigInt(sum);
log << "\n exact a + b = " << exact.str() << "\n TowardsZero = " << stored.str()
<< "\n";
log.flush();
if (scale != MantissaRange::MantissaScale::LargeLegacy)
BEAST_EXPECT(stored == exact);
}
}