diff --git a/include/xrpl/basics/Number.h b/include/xrpl/basics/Number.h index 846607066b..a012251cc2 100644 --- a/include/xrpl/basics/Number.h +++ b/include/xrpl/basics/Number.h @@ -753,7 +753,7 @@ power(Number const& f, unsigned n); // logarithm with base 10 Number -lg(Number const& value); +log10(Number const& value, int iterations = 50); // Returns f^(1/d) // Uses Newton–Raphson iterations until the result stops changing diff --git a/src/libxrpl/basics/Number.cpp b/src/libxrpl/basics/Number.cpp index 2fcd88f3b0..6dd7713235 100644 --- a/src/libxrpl/basics/Number.cpp +++ b/src/libxrpl/basics/Number.cpp @@ -984,46 +984,71 @@ power(Number const& f, unsigned n) return r; } -// Continued fraction approximation of ln(x) +// Series expansion method approximation of ln(x) static Number -ln(Number const& x, unsigned iterations = 50) +ln(Number const& x, int iterations = 50) { + static Number const N0(0); + static Number const N2(2, 0); + static Number const N05(5, -1); + static Number const LN2(693'147'180'559'945'309ll, -18); + if (x <= 0) - throw std::runtime_error("Not positive value"); + throw std::runtime_error("Not a positive value"); + else if (x == 1) + return N0; - Number const z = (x - 1) / (x + 1); - Number const zz = z * z; - Number denom = Number(1, -10); + int exponent = 0; + Number mantissa = x; - // Construct the fraction from the bottom up - for (int i = iterations; i > 0; --i) + while (mantissa >= N2) { - Number k(2 * i - 1); - denom = k - (i * i * zz / denom); + mantissa /= 2; + exponent += 1; + } + while (mantissa < N05) + { + mantissa *= 2; + exponent -= 1; } - auto const r = 2 * z / denom; - return r; + Number z = (mantissa - 1) / (mantissa + 1); + Number const zz = z * z; + Number sum; + + for (int i = 1; i <= iterations; ++i) + { + sum = sum + z / (2 * i - 1); + z = z * zz; + } + + return 2 * sum + exponent * LN2; } Number -lg(Number const& x) +log10(Number const& x, int iterations) { - static Number const ln10 = ln(Number(10)); + static Number const N0(0); + static Number const LN10(2'302'585'092'994'046ll, -15); + + if (x <= 0) + throw std::runtime_error("Not a positive value"); + else if (x == 1) + return N0; if (x <= Number(10)) { - auto const r = ln(x) / ln10; + auto const r = ln(x, iterations) / LN10; return r; } - // ln(x) = ln(normX * 10^norm) = ln(normX) + norm * ln(10) + // (1 <= normalX < 10) + // ln(x) = ln(normalX * 10^norm) = ln(normalX) + norm * ln(10) int diffExp = 15 + x.exponent(); - Number const normalX = x / Number(1, diffExp); // (1 <= normalX < 10) - auto const lnX = ln(normalX) + diffExp * ln10; - - auto const r = lnX / ln10; - return r; + Number const normalX = x / Number(1, diffExp); + auto const lnX = ln(normalX, iterations) + diffExp * LN10; + auto const lgX = lnX / LN10; + return lgX; } // Returns f^(1/d) diff --git a/src/test/app/HostFuncImpl_test.cpp b/src/test/app/HostFuncImpl_test.cpp index fae904c0ec..80dcd2970d 100644 --- a/src/test/app/HostFuncImpl_test.cpp +++ b/src/test/app/HostFuncImpl_test.cpp @@ -3067,7 +3067,7 @@ struct HostFuncImpl_test : public beast::unit_test::suite { auto const x = hfs.floatSet(1, -2, 0); // 0.01 auto const y = - hfs.floatSet(-1999999993734431, -15, 0); // almost -2 + hfs.floatSet(-2'000'000'000'000'000ll, -15, 0); // -2 if (BEAST_EXPECT(x && y)) { auto const result = hfs.floatLog(makeSlice(*x), 0); diff --git a/src/test/basics/Number_test.cpp b/src/test/basics/Number_test.cpp index 1fa5ae6e8f..f647a7abd0 100644 --- a/src/test/basics/Number_test.cpp +++ b/src/test/basics/Number_test.cpp @@ -1712,6 +1712,89 @@ public: } } + void + test_log10() + { + auto const scale = Number::getMantissaScale(); + testcase << "test_lg " << to_string(scale); + + using Case = std::tuple; + auto test = [this](auto const& c) { + for (auto const& [x, z] : c) + { + auto const result = log10(x); + std::stringstream ss; + ss << "lg(" << x << ") = " << result << ". Expected: " << z; + // std::cout << ss.str() << std::endl; + BEAST_EXPECTS(result == z, ss.str()); + } + }; + + auto const cSmall = std::to_array( + {{Number{2}, Number{3'010'299'956'639'811ll, -16}}, + {Number{2'000'000}, Number{6'301'029'995'663'985ll, -15}}, + {Number{2, -30}, Number{-2'969'897'000'433'602ll, -14}}, + {Number{1}, Number{0}}, + {Number{1'000'000'000'000'000ll}, Number{15}}, + {Number{5625, -4}, Number{-2'498'774'732'165'998, -16}}}); + + auto const cLarge = std::to_array( + {{Number{ + false, Number::maxMantissa() - 9, -1, Number::normalized{}}, + Number{ + false, + 1'746'901'684'478'673'451ll, + -17, + Number::normalized{}}}, + {Number{false, Number::maxMantissa() - 9, 0, Number::normalized{}}, + Number{ + false, + 1'846'901'684'478'673'451ll, + -17, + Number::normalized{}}}, + {Number{Number::maxRep}, + Number{ + false, + 1'861'728'612'932'620'011ll, + -17, + Number::normalized{}}}}); + + if (Number::getMantissaScale() == MantissaRange::small) + { + test(cSmall); + } + else + { + NumberRoundModeGuard mg(Number::towards_zero); + test(cLarge); + } + + { + bool caught = false; + try + { + log10(Number{-2}); + } + catch (std::runtime_error const&) + { + caught = true; + } + BEAST_EXPECT(caught); + caught = false; + + try + { + log10(Number()); + } + catch (std::runtime_error const&) + { + caught = true; + } + BEAST_EXPECT(caught); + caught = false; + } + } + void run() override { @@ -1739,6 +1822,7 @@ public: test_truncate(); testRounding(); testInt64(); + test_log10(); } } }; diff --git a/src/xrpld/app/wasm/detail/HostFuncImplFloat.cpp b/src/xrpld/app/wasm/detail/HostFuncImplFloat.cpp index e424e3cedf..d0f8f1e8e9 100644 --- a/src/xrpld/app/wasm/detail/HostFuncImplFloat.cpp +++ b/src/xrpld/app/wasm/detail/HostFuncImplFloat.cpp @@ -454,7 +454,7 @@ floatLogImpl(Slice const& x, int32_t mode) if (!xx) return Unexpected(HostFunctionError::FLOAT_INPUT_MALFORMED); - detail::Number2 res(lg(xx)); + detail::Number2 res(log10(xx)); return res.toBytes(); }