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.

1 files changed, 41 insertions, 21 deletions

diff --git a/math.c b/math.c
index 51a9eac8b4..1e2af7aa75 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