[Guile-commits] GNU Guile branch, stable-2.0, updated. v2.0.7-211-g4400266

[Mark H Weaver]

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=4400266478b4a477c6747f9eed38f7c6021491d8

The branch, stable-2.0 has been updated
       via  4400266478b4a477c6747f9eed38f7c6021491d8 (commit)
       via  c8248c8ed5459991e7d2d6d8f20f652295c19514 (commit)
      from  1d64b4edb9da4011ad06c0fab1c6225ec20b0876 (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 4400266478b4a477c6747f9eed38f7c6021491d8
Author: Mark H Weaver <address@hidden>
Date:   Tue Mar 19 18:48:56 2013 -0400

    Sqrt returns exact results when possible.
    
    * libguile/numbers.c (scm_sqrt): Handle exact integers and rationals in
      such a way that exact results are returned whenever possible.
    
    * test-suite/tests/numbers.test ("sqrt"): Add tests.

commit c8248c8ed5459991e7d2d6d8f20f652295c19514
Author: Mark H Weaver <address@hidden>
Date:   Tue Mar 19 22:38:45 2013 -0400

    Optimize scm_i_divide2double for integers less than 2^DBL_MANT_DIG.
    
    * libguile/numbers.c (scm_i_divide2double): Optimize for common case
      when both operands are less than 2^DBL_MANT_DIG (normally 2^53).

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

Summary of changes:
 libguile/numbers.c            |   81 ++++++++++++++++++++++++++++++++++++++---
 test-suite/tests/numbers.test |   40 ++++++++++++++++++++-
 2 files changed, 114 insertions(+), 7 deletions(-)

diff --git a/libguile/numbers.c b/libguile/numbers.c
index 1f845a3..9725fe4 100644
--- a/libguile/numbers.c
+++ b/libguile/numbers.c
@@ -475,8 +475,17 @@ scm_i_divide2double (SCM n, SCM d)
   mpz_t nn, dd, lo, hi, x;
   ssize_t e;
 
-  if (SCM_I_INUMP (d))
+  if (SCM_LIKELY (SCM_I_INUMP (d)))
     {
+      if (SCM_LIKELY (SCM_I_INUMP (n)
+                      && (SCM_I_FIXNUM_BIT-1 <= DBL_MANT_DIG
+                          || (SCM_I_INUM (n) < (1L << DBL_MANT_DIG)
+                              && SCM_I_INUM (d) < (1L << DBL_MANT_DIG)))))
+        /* If both N and D can be losslessly converted to doubles, then
+           we can rely on IEEE floating point to do proper rounding much
+           faster than we can. */
+        return ((double) SCM_I_INUM (n)) / ((double) SCM_I_INUM (d));
+
       if (SCM_UNLIKELY (scm_is_eq (d, SCM_INUM0)))
         {
           if (scm_is_true (scm_positive_p (n)))
@@ -486,6 +495,7 @@ scm_i_divide2double (SCM n, SCM d)
           else
             return 0.0 / 0.0;
         }
+
       mpz_init_set_si (dd, SCM_I_INUM (d));
     }
   else
@@ -9978,11 +9988,70 @@ SCM_PRIMITIVE_GENERIC (scm_sqrt, "sqrt", 1, 0, 0,
     }
   else if (SCM_NUMBERP (z))
     {
-      double xx = scm_to_double (z);
-      if (xx < 0)
-        return scm_c_make_rectangular (0.0, sqrt (-xx));
-      else
-        return scm_from_double (sqrt (xx));
+      if (SCM_I_INUMP (z))
+        {
+          if (SCM_I_INUM (z) >= 0)
+            {
+              if (SCM_I_FIXNUM_BIT < DBL_MANT_DIG
+                  || SCM_I_INUM (z) < (1L << (DBL_MANT_DIG - 1)))
+                {
+                  double root = sqrt (SCM_I_INUM (z));
+
+                  /* If 0 <= x < 2^(DBL_MANT_DIG-1) and sqrt(x) is an
+                     integer, then the result is exact. */
+                  if (root == floor (root))
+                    return SCM_I_MAKINUM ((scm_t_inum) root);
+                  else
+                    return scm_from_double (root);
+                }
+              else
+                {
+                  mpz_t x;
+                  scm_t_inum root;
+
+                  mpz_init_set_ui (x, SCM_I_INUM (z));
+                  if (mpz_perfect_square_p (x))
+                    {
+                      mpz_sqrt (x, x);
+                      root = mpz_get_ui (x);
+                      mpz_clear (x);
+                      return SCM_I_MAKINUM (root);
+                    }
+                  else
+                    mpz_clear (x);
+                }
+            }
+        }
+      else if (SCM_BIGP (z))
+        {
+          /* IMPROVE-ME: Handle square roots of very large integers
+             better: (1) integers too large to fit in a double, and
+             (2) integers so large that the roundoff of the original
+             number would significantly reduce precision. */
+
+          if (mpz_sgn (SCM_I_BIG_MPZ (z)) >= 0
+              && mpz_perfect_square_p (SCM_I_BIG_MPZ (z)))
+            {
+              SCM root = scm_i_mkbig ();
+
+              mpz_sqrt (SCM_I_BIG_MPZ (root), SCM_I_BIG_MPZ (z));
+              scm_remember_upto_here_1 (z);
+              return scm_i_normbig (root);
+            }
+        }
+      else if (SCM_FRACTIONP (z))
+        /* FIXME: This loses precision due to double rounding. */
+        return scm_divide (scm_sqrt (SCM_FRACTION_NUMERATOR (z)),
+                           scm_sqrt (SCM_FRACTION_DENOMINATOR (z)));
+
+      /* Fallback method, when the cases above do not apply. */
+      {
+        double xx = scm_to_double (z);
+        if (xx < 0)
+          return scm_c_make_rectangular (0.0, sqrt (-xx));
+        else
+          return scm_from_double (sqrt (xx));
+      }
     }
   else
     SCM_WTA_DISPATCH_1 (g_scm_sqrt, z, 1, s_scm_sqrt);
diff --git a/test-suite/tests/numbers.test b/test-suite/tests/numbers.test
index be2e317..a52e79a 100644
--- a/test-suite/tests/numbers.test
+++ b/test-suite/tests/numbers.test
@@ -4840,7 +4840,45 @@
   (pass-if-exception "two args" exception:wrong-num-args
     (sqrt 123 456))
 
-  (pass-if (eqv? 0.0 (sqrt 0)))
+  (pass-if (eqv? 0 (sqrt 0)))
+  (pass-if (eqv? 1 (sqrt 1)))
+  (pass-if (eqv? 2 (sqrt 4)))
+  (pass-if (eqv? 3 (sqrt 9)))
+  (pass-if (eqv? 4 (sqrt 16)))
+  (pass-if (eqv? fixnum-max (sqrt (expt fixnum-max 2))))
+  (pass-if (eqv? (+ 1 fixnum-max) (sqrt (expt (+ 1 fixnum-max) 2))))
+  (pass-if (eqv? (expt 10 400) (sqrt (expt 10 800))))
+  (pass-if (eqv? (/ (expt 10 1000)
+                    (expt 13 1000))
+                 (sqrt (/ (expt 10 2000)
+                          (expt 13 2000)))))
+
+  (with-test-prefix "exact sqrt"
+
+    (define (test root)
+      (pass-if (list root 'exact)
+        (eqv? root (sqrt (expt root 2))))
+      (pass-if (list root '-1)
+        (let ((r (sqrt (- (expt root 2) 1))))
+          (and (inexact? r)
+               (eqv-loosely? root r))))
+      (pass-if (list root '+1)
+        (let ((r (sqrt (+ (expt root 2) 1))))
+          (and (inexact? r)
+               (eqv-loosely? root r))))
+      (pass-if (list root 'negative)
+        (eqv-loosely? (* +i root) (sqrt (- (expt root 2))))))
+
+    (test (exact-integer-sqrt (+ -1 (expt 2 (+  2 dbl-mant-dig)))))
+    (test (exact-integer-sqrt (+ -1 (expt 2 (+  1 dbl-mant-dig)))))
+    (test (exact-integer-sqrt (+ -1 (expt 2 (+  0 dbl-mant-dig)))))
+    (test (exact-integer-sqrt (+ -1 (expt 2 (+ -1 dbl-mant-dig)))))
+    (test (exact-integer-sqrt (+ -1 (expt 2 (+ -2 dbl-mant-dig))))))
+
+  (pass-if (eqv? +4i (sqrt -16)))
+  (pass-if (eqv-loosely? +1.0e150i (sqrt #e-1e300)))
+  (pass-if (eqv-loosely? +0.7071i (sqrt -1/2)))
+
   (pass-if (eqv? 0.0 (sqrt 0.0)))
   (pass-if (eqv? 1.0 (sqrt 1.0)))
   (pass-if (eqv-loosely? 2.0   (sqrt 4.0)))


hooks/post-receive
-- 
GNU Guile