Skip to content

Commit

Permalink
Replace Number division algorithm
Browse files Browse the repository at this point in the history 
* Replace division with faster algorithm.
* Correct some rounding bugs in multiplication.
* Add tests for rounding bugs.
  • Loading branch information
@HowardHinnant
HowardHinnant committed Dec 22, 2022
1 parent 71e0b70 commit 16f1057
Show file tree
Hide file tree
Showing 2 changed files with 175 additions and 50 deletions.
39 changes: 15 additions & 24 deletions src/ripple/basics/impl/Number.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -185,6 +185,8 @@ Number::normalize()
--exponent_;
}
Guard g;
if (negative)
g.set_negative();
while (m > maxMantissa)
{
if (exponent_ >= maxExponent)
Expand Down Expand Up @@ -364,6 +366,8 @@ Number::operator*=(Number const& y)
auto ze = xe + ye;
auto zn = xn * yn;
Guard g;
if (zn == -1)
g.set_negative();
while (zm > maxMantissa)
{
g.push(static_cast<unsigned>(zm % 10));
Expand Down Expand Up @@ -402,44 +406,31 @@ Number::operator/=(Number const& y)
{
if (y == Number{})
throw std::overflow_error("Number: divide by 0");
if (*this == Number{})
return *this;
int np = 1;
auto nm = mantissa();
auto ne = exponent();
if (nm < 0)
{
nm = -nm;
np = -1;
}
int dp = 1;
auto dm = y.mantissa();
auto de = y.exponent();
if (dm < 0)
{
dm = -dm;
dp = -1;
}
// Divide numerator and denominator such that the
// denominator is in the range [1, 10).
const int offset = -15 - y.exponent();
Number n{nm * (np * dp), exponent() + offset};
Number d{dm, y.exponent() + offset};
// Quadratic least squares fit to 1/x in the range [1, 10]
constexpr Number a0{9178756872006464, -16, unchecked{}};
constexpr Number a1{-2149215784206187, -16, unchecked{}};
constexpr Number a2{1405502114116773, -17, unchecked{}};
static_assert(a0.isnormal());
static_assert(a1.isnormal());
static_assert(a2.isnormal());
Number rm2{};
Number rm1{};
Number r = (a2 * d + a1) * d + a0;
// Newton–Raphson iteration of 1/x - d with initial guess r
// halt when r stops changing, checking for bouncing on the last iteration
do
{
rm2 = rm1;
rm1 = r;
r = r + r * (one - d * r);
} while (r != rm1 && r != rm2);
*this = n * r;
// Shift by 10^17 gives greatest precision while not overflowing uint128_t
// or the cast back to int64_t
const uint128_t f = 100'000'000'000'000'000;
mantissa_ = static_cast<std::int64_t>(uint128_t(nm) * f / uint128_t(dm));
exponent_ = ne - de - 17;
mantissa_ *= np * dp;
normalize();
return *this;
}

Expand Down
186 changes: 160 additions & 26 deletions src/test/basics/Number_test.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -150,25 +150,94 @@ class Number_test : public beast::unit_test::suite
{
testcase("test_mul");
using Case = std::tuple<Number, Number, Number>;
Case c[]{
{Number{7}, Number{8}, Number{56}},
{Number{1414213562373095, -15},
Number{1414213562373095, -15},
Number{2000000000000000, -15}},
{Number{-1414213562373095, -15},
Number{1414213562373095, -15},
Number{-2000000000000000, -15}},
{Number{-1414213562373095, -15},
Number{-1414213562373095, -15},
Number{2000000000000000, -15}},
{Number{3214285714285706, -15},
Number{3111111111111119, -15},
Number{1000000000000000, -14}},
{Number{1000000000000000, -32768},
Number{1000000000000000, -32768},
Number{0}}};
for (auto const& [x, y, z] : c)
BEAST_EXPECT(x * y == z);
saveNumberRoundMode save{Number::setround(Number::to_nearest)};
{
Case c[]{
{Number{7}, Number{8}, Number{56}},
{Number{1414213562373095, -15},
Number{1414213562373095, -15},
Number{2000000000000000, -15}},
{Number{-1414213562373095, -15},
Number{1414213562373095, -15},
Number{-2000000000000000, -15}},
{Number{-1414213562373095, -15},
Number{-1414213562373095, -15},
Number{2000000000000000, -15}},
{Number{3214285714285706, -15},
Number{3111111111111119, -15},
Number{1000000000000000, -14}},
{Number{1000000000000000, -32768},
Number{1000000000000000, -32768},
Number{0}}};
for (auto const& [x, y, z] : c)
BEAST_EXPECT(x * y == z);
}
Number::setround(Number::towards_zero);
{
Case c[]{
{Number{7}, Number{8}, Number{56}},
{Number{1414213562373095, -15},
Number{1414213562373095, -15},
Number{1999999999999999, -15}},
{Number{-1414213562373095, -15},
Number{1414213562373095, -15},
Number{-1999999999999999, -15}},
{Number{-1414213562373095, -15},
Number{-1414213562373095, -15},
Number{1999999999999999, -15}},
{Number{3214285714285706, -15},
Number{3111111111111119, -15},
Number{9999999999999999, -15}},
{Number{1000000000000000, -32768},
Number{1000000000000000, -32768},
Number{0}}};
for (auto const& [x, y, z] : c)
BEAST_EXPECT(x * y == z);
}
Number::setround(Number::downward);
{
Case c[]{
{Number{7}, Number{8}, Number{56}},
{Number{1414213562373095, -15},
Number{1414213562373095, -15},
Number{1999999999999999, -15}},
{Number{-1414213562373095, -15},
Number{1414213562373095, -15},
Number{-2000000000000000, -15}},
{Number{-1414213562373095, -15},
Number{-1414213562373095, -15},
Number{1999999999999999, -15}},
{Number{3214285714285706, -15},
Number{3111111111111119, -15},
Number{9999999999999999, -15}},
{Number{1000000000000000, -32768},
Number{1000000000000000, -32768},
Number{0}}};
for (auto const& [x, y, z] : c)
BEAST_EXPECT(x * y == z);
}
Number::setround(Number::upward);
{
Case c[]{
{Number{7}, Number{8}, Number{56}},
{Number{1414213562373095, -15},
Number{1414213562373095, -15},
Number{2000000000000000, -15}},
{Number{-1414213562373095, -15},
Number{1414213562373095, -15},
Number{-1999999999999999, -15}},
{Number{-1414213562373095, -15},
Number{-1414213562373095, -15},
Number{2000000000000000, -15}},
{Number{3214285714285706, -15},
Number{3111111111111119, -15},
Number{1000000000000000, -14}},
{Number{1000000000000000, -32768},
Number{1000000000000000, -32768},
Number{0}}};
for (auto const& [x, y, z] : c)
BEAST_EXPECT(x * y == z);
}
bool caught = false;
try
{
Expand All @@ -187,13 +256,78 @@ class Number_test : public beast::unit_test::suite
{
testcase("test_div");
using Case = std::tuple<Number, Number, Number>;
Case c[]{
{Number{1}, Number{2}, Number{5, -1}},
{Number{1}, Number{10}, Number{1, -1}},
{Number{1}, Number{-10}, Number{-1, -1}},
{Number{0}, Number{100}, Number{0}}};
for (auto const& [x, y, z] : c)
BEAST_EXPECT(x / y == z);
saveNumberRoundMode save{Number::setround(Number::to_nearest)};
{
Case c[]{
{Number{1}, Number{2}, Number{5, -1}},
{Number{1}, Number{10}, Number{1, -1}},
{Number{1}, Number{-10}, Number{-1, -1}},
{Number{0}, Number{100}, Number{0}},
{Number{1414213562373095, -10},
Number{1414213562373095, -10},
Number{1}},
{Number{9'999'999'999'999'999},
Number{1'000'000'000'000'000},
Number{9'999'999'999'999'999, -15}},
{Number{2}, Number{3}, Number{6'666'666'666'666'667, -16}},
{Number{-2}, Number{3}, Number{-6'666'666'666'666'667, -16}}};
for (auto const& [x, y, z] : c)
BEAST_EXPECT(x / y == z);
}
Number::setround(Number::towards_zero);
{
Case c[]{
{Number{1}, Number{2}, Number{5, -1}},
{Number{1}, Number{10}, Number{1, -1}},
{Number{1}, Number{-10}, Number{-1, -1}},
{Number{0}, Number{100}, Number{0}},
{Number{1414213562373095, -10},
Number{1414213562373095, -10},
Number{1}},
{Number{9'999'999'999'999'999},
Number{1'000'000'000'000'000},
Number{9'999'999'999'999'999, -15}},
{Number{2}, Number{3}, Number{6'666'666'666'666'666, -16}},
{Number{-2}, Number{3}, Number{-6'666'666'666'666'666, -16}}};
for (auto const& [x, y, z] : c)
BEAST_EXPECT(x / y == z);
}
Number::setround(Number::downward);
{
Case c[]{
{Number{1}, Number{2}, Number{5, -1}},
{Number{1}, Number{10}, Number{1, -1}},
{Number{1}, Number{-10}, Number{-1, -1}},
{Number{0}, Number{100}, Number{0}},
{Number{1414213562373095, -10},
Number{1414213562373095, -10},
Number{1}},
{Number{9'999'999'999'999'999},
Number{1'000'000'000'000'000},
Number{9'999'999'999'999'999, -15}},
{Number{2}, Number{3}, Number{6'666'666'666'666'666, -16}},
{Number{-2}, Number{3}, Number{-6'666'666'666'666'667, -16}}};
for (auto const& [x, y, z] : c)
BEAST_EXPECT(x / y == z);
}
Number::setround(Number::upward);
{
Case c[]{
{Number{1}, Number{2}, Number{5, -1}},
{Number{1}, Number{10}, Number{1, -1}},
{Number{1}, Number{-10}, Number{-1, -1}},
{Number{0}, Number{100}, Number{0}},
{Number{1414213562373095, -10},
Number{1414213562373095, -10},
Number{1}},
{Number{9'999'999'999'999'999},
Number{1'000'000'000'000'000},
Number{9'999'999'999'999'999, -15}},
{Number{2}, Number{3}, Number{6'666'666'666'666'667, -16}},
{Number{-2}, Number{3}, Number{-6'666'666'666'666'666, -16}}};
for (auto const& [x, y, z] : c)
BEAST_EXPECT(x / y == z);
}
bool caught = false;
try
{
Expand Down

0 comments on commit 16f1057

Please sign in to comment.