From 16f1057d2ab951204d7c0c4564620caf78273491 Mon Sep 17 00:00:00 2001 From: Howard Hinnant Date: Thu, 22 Dec 2022 12:20:56 -0500 Subject: [PATCH] Replace Number division algorithm * Replace division with faster algorithm. * Correct some rounding bugs in multiplication. * Add tests for rounding bugs. --- src/ripple/basics/impl/Number.cpp | 39 +++---- src/test/basics/Number_test.cpp | 186 +++++++++++++++++++++++++----- 2 files changed, 175 insertions(+), 50 deletions(-) diff --git a/src/ripple/basics/impl/Number.cpp b/src/ripple/basics/impl/Number.cpp index a7a8159fed3..52d2b556d44 100644 --- a/src/ripple/basics/impl/Number.cpp +++ b/src/ripple/basics/impl/Number.cpp @@ -185,6 +185,8 @@ Number::normalize() --exponent_; } Guard g; + if (negative) + g.set_negative(); while (m > maxMantissa) { if (exponent_ >= maxExponent) @@ -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(zm % 10)); @@ -402,8 +406,11 @@ 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; @@ -411,35 +418,19 @@ Number::operator/=(Number const& y) } 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(uint128_t(nm) * f / uint128_t(dm)); + exponent_ = ne - de - 17; + mantissa_ *= np * dp; + normalize(); return *this; } diff --git a/src/test/basics/Number_test.cpp b/src/test/basics/Number_test.cpp index d3ece630e3e..8c12ff7c5e4 100644 --- a/src/test/basics/Number_test.cpp +++ b/src/test/basics/Number_test.cpp @@ -150,25 +150,94 @@ class Number_test : public beast::unit_test::suite { testcase("test_mul"); using Case = std::tuple; - 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 { @@ -187,13 +256,78 @@ class Number_test : public beast::unit_test::suite { testcase("test_div"); using Case = std::tuple; - 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 {