[Guile-commits] GNU Guile branch, master, updated. release_1-9-15-123-g05f167b

[Ludovic Courtès]

This is an automated email from the git hooks/post-receive script. It was
generated because a ref change was pushed to the repository containing
the project "GNU Guile".

http://git.savannah.gnu.org/cgit/guile.git/commit/?id=05f167beb213e34a3ae849e0008c6fad0871e6ae

The branch, master has been updated
       via  05f167beb213e34a3ae849e0008c6fad0871e6ae (commit)
       via  14a01ec1a87e0b21356ccd4d737263b8ec0eba13 (commit)
      from  1d82efbceb6123a4b9130c8d0e2b558f566a023a (commit)

Those revisions listed above that are new to this repository have
not appeared on any other notification email; so we list those
revisions in full, below.

- Log -----------------------------------------------------------------
commit 05f167beb213e34a3ae849e0008c6fad0871e6ae
Author: Mark H Weaver <address@hidden>
Date:   Tue Feb 15 10:37:03 2011 -0500

    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 <address@hidden>

commit 14a01ec1a87e0b21356ccd4d737263b8ec0eba13
Author: Mark H Weaver <address@hidden>
Date:   Mon Feb 14 18:18:52 2011 -0500

    Fix comment above number-theoretic division tests
    
    * test-suite/tests/numbers.test: Fix comment.
    
    Signed-off-by: Ludovic Courtès <address@hidden>

-----------------------------------------------------------------------

Summary of changes:
 libguile/numbers.c            |  108 ++++++++++++++++++++++++++++++++++-------
 test-suite/tests/numbers.test |   94 +++++++++++++++++++++++++++++-------
 2 files changed, 166 insertions(+), 36 deletions(-)

diff --git a/libguile/numbers.c b/libguile/numbers.c
index 7c4ea1b..d0aacb7 100644
--- 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 1f2ee03..cb582ed 100644
--- 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))))
 
 ;;;
@@ -4512,12 +4556,26 @@
 
 
 ;;;
+;;; Tests for number-theoretic division operators:
+;;;
 ;;; euclidean/
 ;;; euclidean-quotient
 ;;; euclidean-remainder
+;;; floor/
+;;; floor-quotient
+;;; floor-remainder
+;;; ceiling/
+;;; ceiling-quotient
+;;; ceiling-remainder
+;;; truncate/
+;;; truncate-quotient
+;;; truncate-remainder
 ;;; centered/
 ;;; centered-quotient
 ;;; centered-remainder
+;;; round/
+;;; round-quotient
+;;; round-remainder
 ;;;
 
 (with-test-prefix "Number-theoretic division"


hooks/post-receive
-- 
GNU Guile