From da39936ce165ea9462b9e192eb6b608485c94842 Mon Sep 17 00:00:00 2001
From: Nobuyoshi Nakada <nobu@ruby-lang.org>
Date: Sun, 16 Jul 2023 01:24:44 +0900
Subject: [PATCH] Prefer integer as base of intermediate logarithms

As long as "floating point numbers" cannot accurately represent an
irrational number, the result of the natural logarithm cannot be
accurate.  Logarithms with an integer base may have the possibility to
represent more accurately.
---
 math.c | 62 ++++++++++++++++++++++++++++++++++++++--------------------
 1 file changed, 41 insertions(+), 21 deletions(-)

diff --git a/math.c b/math.c
index 51a9eac8b4a661..1e2af7aa75edd3 100644
--- a/math.c
+++ b/math.c
@@ -474,7 +474,6 @@ math_exp(VALUE unused_obj, VALUE x)
 # define M_LN10 2.30258509299404568401799145468436421
 #endif
 
-static double math_log1(VALUE x);
 FUNC_MINIMIZED(static VALUE math_log(int, const VALUE *, VALUE));
 
 /*
@@ -509,20 +508,6 @@ math_log(int argc, const VALUE *argv, VALUE unused_obj)
     return rb_math_log(argc, argv);
 }
 
-VALUE
-rb_math_log(int argc, const VALUE *argv)
-{
-    VALUE x, base;
-    double d;
-
-    rb_scan_args(argc, argv, "11", &x, &base);
-    d = math_log1(x);
-    if (argc == 2) {
-        d /= math_log1(base);
-    }
-    return DBL2NUM(d);
-}
-
 static double
 get_double_rshift(VALUE x, size_t *pnumbits)
 {
@@ -541,16 +526,51 @@ get_double_rshift(VALUE x, size_t *pnumbits)
 }
 
 static double
-math_log1(VALUE x)
+math_log_split(VALUE x, size_t *numbits)
 {
-    size_t numbits;
-    double d = get_double_rshift(x, &numbits);
+    double d = get_double_rshift(x, numbits);
 
     domain_check_min(d, 0.0, "log");
-    /* check for pole error */
-    if (d == 0.0) return -HUGE_VAL;
+    return d;
+}
+
+#if defined(log2) || defined(HAVE_LOG2)
+# define log_intermediate log2
+#else
+# define log_intermediate log10
+#endif
+
+VALUE
+rb_math_log(int argc, const VALUE *argv)
+{
+    VALUE x, base;
+    double d;
+    size_t numbits;
 
-    return log(d) + numbits * M_LN2; /* log(d * 2 ** numbits) */
+    argc = rb_scan_args(argc, argv, "11", &x, &base);
+    d = math_log_split(x, &numbits);
+    if (argc == 2) {
+        size_t numbits_2;
+        double b = math_log_split(base, &numbits_2);
+        /* check for pole error */
+        if (d == 0.0) {
+            if (b > 0.0) return DBL2NUM(HUGE_VAL);
+            if (b < 0.0) return DBL2NUM(-HUGE_VAL);
+            return DBL2NUM(nan(""));
+        }
+        else if (b == 0.0) {
+            return DBL2NUM(-0.0);
+        }
+        d = log_intermediate(d) / log_intermediate(b);
+        numbits -= numbits_2;
+    }
+    else {
+        /* check for pole error */
+        if (d == 0.0) return DBL2NUM(-HUGE_VAL);
+        d = log(d);
+    }
+    d += numbits * M_LN2;
+    return DBL2NUM(d);
 }
 
 #ifndef log2