Switch to series expansion method for ln() (#6268)

* Switch to series expansion method for ln()
Add float lg() tests to Number tests;
* Rename lg -> log10
* Add check for 0 to log10()

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

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)

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);

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<Number, Number>;
+        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<Case>(
+            {{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<Case>(
+            {{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();
         }
     }
 };

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();
     }