Improvements to `log' and `log10'

* libguile/numbers.c (log_of_shifted_double, log_of_exact_integer,
  log_of_exact_integer_with_size, log_of_fraction): New internal static
  functions used by scm_log and scm_log10.

  (scm_log, scm_log10): Robustly handle large integers, large and small
  fractions, and fractions close to 1.  Previously, computing logarithms
  of fractions close to 1 yielded grossly inaccurate results, and the
  other cases yielded infinities even though the answer could easily fit
  in a double.  (log -0.0) now returns -inf.0+<PI>i, where previously it
  returned -inf.0.  (log 0) now throws a numerical overflow exception,
  where previously it returned -inf.0.  (log 0.0) still returns -inf.0.
  Analogous changes made to `log10'.

* test-suite/tests/numbers.test (log, log10): Add tests.

Signed-off-by: Ludovic Courtès <ludo@gnu.org>

diff --git a/libguile/numbers.c b/libguile/numbers.c
index 7c4ea1b..d0aacb7 100644 (file)
--- a/libguile/numbers.c
+++ b/libguile/numbers.c
@@ -111,6 +111,7 @@ typedef scm_t_signed_bits scm_t_inum;
 
 static SCM flo0;
 static SCM exactly_one_half;
+static SCM flo_log10e;
 
 #define SCM_SWAP(x, y) do { SCM __t = x; x = y; y = __t; } while (0)
 
@@ -9372,6 +9373,62 @@ scm_is_number (SCM z)
 }
 
 
+/* Returns log(x * 2^shift) */
+static SCM
+log_of_shifted_double (double x, long shift)
+{
+  double ans = log (fabs (x)) + shift * M_LN2;
+
+  if (x > 0.0 || double_is_non_negative_zero (x))
+    return scm_from_double (ans);
+  else
+    return scm_c_make_rectangular (ans, M_PI);
+}
+
+/* Returns log(n), for exact integer n of integer-length size */
+static SCM
+log_of_exact_integer_with_size (SCM n, long size)
+{
+  long shift = size - 2 * scm_dblprec[0];
+
+  if (shift > 0)
+    return log_of_shifted_double
+      (scm_to_double (scm_ash (n, scm_from_long(-shift))),
+       shift);
+  else
+    return log_of_shifted_double (scm_to_double (n), 0);
+}
+
+/* Returns log(n), for exact integer n of integer-length size */
+static SCM
+log_of_exact_integer (SCM n)
+{
+  return log_of_exact_integer_with_size
+    (n, scm_to_long (scm_integer_length (n)));
+}
+
+/* Returns log(n/d), for exact non-zero integers n and d */
+static SCM
+log_of_fraction (SCM n, SCM d)
+{
+  long n_size = scm_to_long (scm_integer_length (n));
+  long d_size = scm_to_long (scm_integer_length (d));
+
+  if (abs (n_size - d_size) > 1)
+    return (scm_difference (log_of_exact_integer_with_size (n, n_size),
+                           log_of_exact_integer_with_size (d, d_size)));
+  else if (scm_is_false (scm_negative_p (n)))
+    return scm_from_double
+      (log1p (scm_to_double (scm_divide2real (scm_difference (n, d), d))));
+  else
+    return scm_c_make_rectangular
+      (log1p (scm_to_double (scm_divide2real
+                            (scm_difference (scm_abs (n), d),
+                             d))),
+       M_PI);
+}
+
+
 /* In the following functions we dispatch to the real-arg funcs like log()
    when we know the arg is real, instead of just handing everything to
    clog() for instance.  This is in case clog() doesn't optimize for a
@@ -9394,17 +9451,21 @@ SCM_PRIMITIVE_GENERIC (scm_log, "log", 1, 0, 0,
                                      atan2 (im, re));
 #endif
     }
-  else if (SCM_NUMBERP (z))
+  else if (SCM_REALP (z))
+    return log_of_shifted_double (SCM_REAL_VALUE (z), 0);
+  else if (SCM_I_INUMP (z))
     {
-      /* ENHANCE-ME: When z is a bignum the logarithm will fit a double
-         although the value itself overflows.  */
-      double re = scm_to_double (z);
-      double l = log (fabs (re));
-      if (re >= 0.0)
-        return scm_from_double (l);
-      else
-        return scm_c_make_rectangular (l, M_PI);
+#ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
+      if (scm_is_eq (z, SCM_INUM0))
+       scm_num_overflow (s_scm_log);
+#endif
+      return log_of_shifted_double (SCM_I_INUM (z), 0);
     }
+  else if (SCM_BIGP (z))
+    return log_of_exact_integer (z);
+  else if (SCM_FRACTIONP (z))
+    return log_of_fraction (SCM_FRACTION_NUMERATOR (z),
+                           SCM_FRACTION_DENOMINATOR (z));
   else
     SCM_WTA_DISPATCH_1 (g_scm_log, z, 1, s_scm_log);
 }
@@ -9431,17 +9492,27 @@ SCM_PRIMITIVE_GENERIC (scm_log10, "log10", 1, 0, 0,
                                      M_LOG10E * atan2 (im, re));
 #endif
     }
-  else if (SCM_NUMBERP (z))
+  else if (SCM_REALP (z) || SCM_I_INUMP (z))
     {
-      /* ENHANCE-ME: When z is a bignum the logarithm will fit a double
-         although the value itself overflows.  */
-      double re = scm_to_double (z);
-      double l = log10 (fabs (re));
-      if (re >= 0.0)
-        return scm_from_double (l);
-      else
-        return scm_c_make_rectangular (l, M_LOG10E * M_PI);
+#ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
+      if (scm_is_eq (z, SCM_INUM0))
+       scm_num_overflow (s_scm_log10);
+#endif
+      {
+       double re = scm_to_double (z);
+       double l = log10 (fabs (re));
+       if (re > 0.0 || double_is_non_negative_zero (re))
+         return scm_from_double (l);
+       else
+         return scm_c_make_rectangular (l, M_LOG10E * M_PI);
+      }
     }
+  else if (SCM_BIGP (z))
+    return scm_product (flo_log10e, log_of_exact_integer (z));
+  else if (SCM_FRACTIONP (z))
+    return scm_product (flo_log10e,
+                       log_of_fraction (SCM_FRACTION_NUMERATOR (z),
+                                        SCM_FRACTION_DENOMINATOR (z)));
   else
     SCM_WTA_DISPATCH_1 (g_scm_log10, z, 1, s_scm_log10);
 }
@@ -9536,6 +9607,7 @@ scm_init_numbers ()
   scm_add_feature ("complex");
   scm_add_feature ("inexact");
   flo0 = scm_from_double (0.0);
+  flo_log10e = scm_from_double (M_LOG10E);
 
   /* determine floating point precision */
   for (i=2; i <= SCM_MAX_DBL_RADIX; ++i)

diff --git a/test-suite/tests/numbers.test b/test-suite/tests/numbers.test
index 9e9728f..cb582ed 100644 (file)
--- a/test-suite/tests/numbers.test
+++ b/test-suite/tests/numbers.test
@@ -4323,14 +4323,36 @@
     (log))
   (pass-if-exception "two args" exception:wrong-num-args
     (log 123 456))
-
-  (pass-if (negative-infinity? (log 0)))
-  (pass-if (negative-infinity? (log 0.0)))
-  (pass-if (eqv? 0.0 (log 1)))
-  (pass-if (eqv? 0.0 (log 1.0)))
-  (pass-if (eqv-loosely? 1.0  (log const-e)))
-  (pass-if (eqv-loosely? 2.0  (log const-e^2)))
-  (pass-if (eqv-loosely? -1.0 (log const-1/e)))
+  (pass-if-exception "(log 0)" exception:numerical-overflow
+    (log 0))
+
+  (pass-if (test-eqv? -inf.0 (log 0.0)))
+  (pass-if (test-eqv? +inf.0 (log +inf.0)))
+  (pass-if (test-eqv? -inf.0+3.14159265358979i (log -0.0)))
+  (pass-if (test-eqv? +inf.0+3.14159265358979i (log -inf.0)))
+  (pass-if (test-eqv?  0.0 (log 1  )))
+  (pass-if (test-eqv?  0.0 (log 1.0)))
+  (pass-if (test-eqv?  1.0 (log const-e)))
+  (pass-if (test-eqv?  2.0 (log const-e^2)))
+  (pass-if (test-eqv? -1.0 (log const-1/e)))
+  (pass-if (test-eqv? -1.0+3.14159265358979i (log (- const-1/e))))
+  (pass-if (test-eqv?  2.30258509299405 (log 10)))
+  (pass-if (test-eqv?  2.30258509299405+3.14159265358979i (log -10)))
+
+  (pass-if (test-eqv?  1.0+0.0i (log (+ const-e +0.0i))))
+  (pass-if (test-eqv?  1.0-0.0i (log (+ const-e -0.0i))))
+
+  (pass-if (eqv-loosely?  230258.509299405 (log (expt 10  100000))))
+  (pass-if (eqv-loosely? -230258.509299405 (log (expt 10 -100000))))
+  (pass-if (eqv-loosely?  230257.410687116 (log (/ (expt 10 100000) 3))))
+  (pass-if (eqv-loosely?  230258.509299405+3.14159265358979i
+                          (log (- (expt 10 100000)))))
+  (pass-if (eqv-loosely? -230258.509299405+3.14159265358979i
+                          (log (- (expt 10 -100000)))))
+  (pass-if (eqv-loosely?  230257.410687116+3.14159265358979i
+                          (log (- (/ (expt 10 100000) 3)))))
+  (pass-if (test-eqv?  3.05493636349961e-151
+                       (log (/ (1+ (expt 2 500)) (expt 2 500)))))
 
   (pass-if (eqv-loosely? 1.0+1.57079i (log 0+2.71828i)))
   (pass-if (eqv-loosely? 1.0-1.57079i (log 0-2.71828i)))
@@ -4350,20 +4372,42 @@
     (log10))
   (pass-if-exception "two args" exception:wrong-num-args
     (log10 123 456))
-
-  (pass-if (negative-infinity? (log10 0)))
-  (pass-if (negative-infinity? (log10 0.0)))
-  (pass-if (eqv? 0.0 (log10 1)))
-  (pass-if (eqv? 0.0 (log10 1.0)))
-  (pass-if (eqv-loosely? 1.0  (log10 10.0)))
-  (pass-if (eqv-loosely? 2.0  (log10 100.0)))
-  (pass-if (eqv-loosely? -1.0 (log10 0.1)))
+  (pass-if-exception "(log10 0)" exception:numerical-overflow
+    (log10 0))
+
+  (pass-if (test-eqv? -inf.0 (log10 0.0)))
+  (pass-if (test-eqv? +inf.0 (log10 +inf.0)))
+  (pass-if (test-eqv? -inf.0+1.36437635384184i (log10 -0.0)))
+  (pass-if (test-eqv? +inf.0+1.36437635384184i (log10 -inf.0)))
+  (pass-if (test-eqv?  0.0 (log10   1  )))
+  (pass-if (test-eqv?  0.0 (log10   1.0)))
+  (pass-if (test-eqv?  1.0 (log10  10  )))
+  (pass-if (test-eqv?  1.0 (log10  10.0)))
+  (pass-if (test-eqv?  2.0 (log10 100.0)))
+  (pass-if (test-eqv? -1.0 (log10   0.1)))
+  (pass-if (test-eqv? -1.0+1.36437635384184i (log10  -0.1)))
+  (pass-if (test-eqv?  1.0+1.36437635384184i (log10 -10  )))
+
+  (pass-if (test-eqv?  1.0+0.0i (log10  10.0+0.0i)))
+  (pass-if (test-eqv?  1.0-0.0i (log10  10.0-0.0i)))
+
+  (pass-if (eqv-loosely?  100000.0 (log10 (expt 10  100000))))
+  (pass-if (eqv-loosely? -100000.0 (log10 (expt 10 -100000))))
+  (pass-if (eqv-loosely?   99999.5228787453 (log10 (/ (expt 10 100000) 3))))
+  (pass-if (eqv-loosely?  100000.0+1.36437635384184i
+                          (log10 (- (expt 10 100000)))))
+  (pass-if (eqv-loosely? -100000.0+1.36437635384184i
+                          (log10 (- (expt 10 -100000)))))
+  (pass-if (eqv-loosely?   99999.5228787453+1.36437635384184i
+                          (log10 (- (/ (expt 10 100000) 3)))))
+  (pass-if (test-eqv?  1.32674200523347e-151
+                       (log10 (/ (1+ (expt 2 500)) (expt 2 500)))))
 
   (pass-if (eqv-loosely? 1.0+0.68218i (log10 0+10.0i)))
   (pass-if (eqv-loosely? 1.0-0.68218i (log10 0-10.0i)))
 
-  (pass-if (eqv-loosely? 0.0+1.36437i (log10 -1)))
-  (pass-if (eqv-loosely? 1.0+1.36437i (log10 -10)))
+  (pass-if (eqv-loosely? 0.0+1.36437i (log10   -1)))
+  (pass-if (eqv-loosely? 1.0+1.36437i (log10  -10)))
   (pass-if (eqv-loosely? 2.0+1.36437i (log10 -100))))
 
 ;;;