[Guile-commits] GNU Guile branch, stable-2.0, updated. v2.0.7-213-gddb7174

[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=ddb717423619cb2c36fb798dc12552b70cd9b0ad

The branch, stable-2.0 has been updated
       via  ddb717423619cb2c36fb798dc12552b70cd9b0ad (commit)
       via  687a87bf012f0c0afa79dd9bebf7d173d1243880 (commit)
      from  4400266478b4a477c6747f9eed38f7c6021491d8 (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 ddb717423619cb2c36fb798dc12552b70cd9b0ad
Author: Mark H Weaver <address@hidden>
Date:   Wed Mar 20 06:15:32 2013 -0400

    Improve sqrt handling of large integers and large and small rationals.
    
    * libguile/numbers.c (exact_integer_is_perfect_square,
      exact_integer_floor_square_root): New static functions.
    
      (scm_sqrt): Use SCM_LIKELY.  Add 'scm_t_inum' variable in inum case to
      reduce the number of uses of SCM_I_INUM.  Rename 'mpz_t' variable.
      Remove unneeded sign check.  Handle bignums too large to fit in a
      double.  Handle fractions too large or too small to fit in a
      normalized double.
    
    * test-suite/tests/numbers.test ("sqrt"): Add tests.

commit 687a87bf012f0c0afa79dd9bebf7d173d1243880
Author: Mark H Weaver <address@hidden>
Date:   Wed Mar 20 02:27:10 2013 -0400

    Optimize inum case of exact-integer-sqrt.
    
    * libguile/numbers.c (scm_exact_integer_sqrt): Use GMP for inum
      case.  It is faster than what we had before.

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

Summary of changes:
 libguile/numbers.c            |  157 +++++++++++++++++++++++++++++++----------
 test-suite/tests/numbers.test |   37 +++++++++-
 2 files changed, 156 insertions(+), 38 deletions(-)

diff --git a/libguile/numbers.c b/libguile/numbers.c
index 9725fe4..a7c0928 100644
--- a/libguile/numbers.c
+++ b/libguile/numbers.c
@@ -9919,25 +9919,17 @@ scm_exact_integer_sqrt (SCM k, SCM *sp, SCM *rp)
 {
   if (SCM_LIKELY (SCM_I_INUMP (k)))
     {
-      scm_t_inum kk = SCM_I_INUM (k);
-      scm_t_inum uu = kk;
-      scm_t_inum ss;
+      mpz_t kk, ss, rr;
 
-      if (SCM_LIKELY (kk > 0))
-       {
-         do
-           {
-             ss = uu;
-             uu = (ss + kk/ss) / 2;
-           } while (uu < ss);
-         *sp = SCM_I_MAKINUM (ss);
-         *rp = SCM_I_MAKINUM (kk - ss*ss);
-       }
-      else if (SCM_LIKELY (kk == 0))
-       *sp = *rp = SCM_INUM0;
-      else
+      if (SCM_I_INUM (k) < 0)
        scm_wrong_type_arg_msg ("exact-integer-sqrt", SCM_ARG1, k,
                                "exact non-negative integer");
+      mpz_init_set_ui (kk, SCM_I_INUM (k));
+      mpz_inits (ss, rr, NULL);
+      mpz_sqrtrem (ss, rr, kk);
+      *sp = SCM_I_MAKINUM (mpz_get_ui (ss));
+      *rp = SCM_I_MAKINUM (mpz_get_ui (rr));
+      mpz_clears (kk, ss, rr, NULL);
     }
   else if (SCM_LIKELY (SCM_BIGP (k)))
     {
@@ -9958,6 +9950,56 @@ scm_exact_integer_sqrt (SCM k, SCM *sp, SCM *rp)
                            "exact non-negative integer");
 }
 
+/* Return true iff K is a perfect square.
+   K must be an exact integer. */
+static int
+exact_integer_is_perfect_square (SCM k)
+{
+  int result;
+
+  if (SCM_LIKELY (SCM_I_INUMP (k)))
+    {
+      mpz_t kk;
+
+      mpz_init_set_si (kk, SCM_I_INUM (k));
+      result = mpz_perfect_square_p (kk);
+      mpz_clear (kk);
+    }
+  else
+    {
+      result = mpz_perfect_square_p (SCM_I_BIG_MPZ (k));
+      scm_remember_upto_here_1 (k);
+    }
+  return result;
+}
+
+/* Return the floor of the square root of K.
+   K must be an exact integer. */
+static SCM
+exact_integer_floor_square_root (SCM k)
+{
+  if (SCM_LIKELY (SCM_I_INUMP (k)))
+    {
+      mpz_t kk;
+      scm_t_inum ss;
+
+      mpz_init_set_ui (kk, SCM_I_INUM (k));
+      mpz_sqrt (kk, kk);
+      ss = mpz_get_ui (kk);
+      mpz_clear (kk);
+      return SCM_I_MAKINUM (ss);
+    }
+  else
+    {
+      SCM s;
+
+      s = scm_i_mkbig ();
+      mpz_sqrt (SCM_I_BIG_MPZ (s), SCM_I_BIG_MPZ (k));
+      scm_remember_upto_here_1 (k);
+      return scm_i_normbig (s);
+    }
+}
+
 
 SCM_PRIMITIVE_GENERIC (scm_sqrt, "sqrt", 1, 0, 0,
                       (SCM z),
@@ -9990,12 +10032,14 @@ SCM_PRIMITIVE_GENERIC (scm_sqrt, "sqrt", 1, 0, 0,
     {
       if (SCM_I_INUMP (z))
         {
-          if (SCM_I_INUM (z) >= 0)
+          scm_t_inum x = SCM_I_INUM (z);
+
+          if (SCM_LIKELY (x >= 0))
             {
-              if (SCM_I_FIXNUM_BIT < DBL_MANT_DIG
-                  || SCM_I_INUM (z) < (1L << (DBL_MANT_DIG - 1)))
+              if (SCM_LIKELY (SCM_I_FIXNUM_BIT < DBL_MANT_DIG
+                              || x < (1L << (DBL_MANT_DIG - 1))))
                 {
-                  double root = sqrt (SCM_I_INUM (z));
+                  double root = sqrt (x);
 
                   /* If 0 <= x < 2^(DBL_MANT_DIG-1) and sqrt(x) is an
                      integer, then the result is exact. */
@@ -10006,31 +10050,25 @@ SCM_PRIMITIVE_GENERIC (scm_sqrt, "sqrt", 1, 0, 0,
                 }
               else
                 {
-                  mpz_t x;
+                  mpz_t xx;
                   scm_t_inum root;
 
-                  mpz_init_set_ui (x, SCM_I_INUM (z));
-                  if (mpz_perfect_square_p (x))
+                  mpz_init_set_ui (xx, x);
+                  if (mpz_perfect_square_p (xx))
                     {
-                      mpz_sqrt (x, x);
-                      root = mpz_get_ui (x);
-                      mpz_clear (x);
+                      mpz_sqrt (xx, xx);
+                      root = mpz_get_ui (xx);
+                      mpz_clear (xx);
                       return SCM_I_MAKINUM (root);
                     }
                   else
-                    mpz_clear (x);
+                    mpz_clear (xx);
                 }
             }
         }
       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)))
+          if (mpz_perfect_square_p (SCM_I_BIG_MPZ (z)))
             {
               SCM root = scm_i_mkbig ();
 
@@ -10038,11 +10076,56 @@ SCM_PRIMITIVE_GENERIC (scm_sqrt, "sqrt", 1, 0, 0,
               scm_remember_upto_here_1 (z);
               return scm_i_normbig (root);
             }
+          else
+            {
+              long expon;
+              double signif = scm_i_big2dbl_2exp (z, &expon);
+
+              if (expon & 1)
+                {
+                  signif *= 2;
+                  expon--;
+                }
+              if (signif < 0)
+                return scm_c_make_rectangular
+                  (0.0, ldexp (sqrt (-signif), expon / 2));
+              else
+                return scm_from_double (ldexp (sqrt (signif), expon / 2));
+            }
         }
       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)));
+        {
+          SCM n = SCM_FRACTION_NUMERATOR (z);
+          SCM d = SCM_FRACTION_DENOMINATOR (z);
+
+          if (exact_integer_is_perfect_square (n)
+              && exact_integer_is_perfect_square (d))
+            return scm_i_make_ratio_already_reduced
+              (exact_integer_floor_square_root (n),
+               exact_integer_floor_square_root (d));
+          else
+            {
+              double xx = scm_i_divide2double (n, d);
+              double abs_xx = fabs (xx);
+              long shift = 0;
+
+              if (SCM_UNLIKELY (abs_xx > DBL_MAX || abs_xx < DBL_MIN))
+                {
+                  shift = (scm_to_long (scm_integer_length (n))
+                           - scm_to_long (scm_integer_length (d))) / 2;
+                  if (shift > 0)
+                    d = left_shift_exact_integer (d, 2 * shift);
+                  else
+                    n = left_shift_exact_integer (n, -2 * shift);
+                  xx = scm_i_divide2double (n, d);
+                }
+
+              if (xx < 0)
+                return scm_c_make_rectangular (0.0, ldexp (sqrt (-xx), shift));
+              else
+                return scm_from_double (ldexp (sqrt (xx), shift));
+            }
+        }
 
       /* Fallback method, when the cases above do not apply. */
       {
diff --git a/test-suite/tests/numbers.test b/test-suite/tests/numbers.test
index a52e79a..7d30392 100644
--- a/test-suite/tests/numbers.test
+++ b/test-suite/tests/numbers.test
@@ -4858,6 +4858,10 @@
     (define (test root)
       (pass-if (list root 'exact)
         (eqv? root (sqrt (expt root 2))))
+      (pass-if (list root '*2)
+        (let ((r (sqrt (* 2 (expt root 2)))))
+          (and (inexact? r)
+               (eqv-loosely? (* (sqrt 2) root) r))))
       (pass-if (list root '-1)
         (let ((r (sqrt (- (expt root 2) 1))))
           (and (inexact? r)
@@ -4873,7 +4877,38 @@
     (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))))))
+    (test (exact-integer-sqrt (+ -1 (expt 2 (+ -2 dbl-mant-dig)))))
+
+    ;; largest finite inexact
+    (test (* (- (expt 2 dbl-mant-dig) 1)
+             (expt 2 (- dbl-max-exp dbl-mant-dig)))))
+
+  (pass-if-equal "smallest inexact"
+      (expt 2.0 (- dbl-min-exp dbl-mant-dig))
+    (sqrt (/ (+ -1 (expt 2 (* 2 (- dbl-mant-dig dbl-min-exp)))))))
+
+  (with-test-prefix "extreme ratios"
+    (define-syntax-rule (test want x)
+      (pass-if 'x
+        (let ((got (sqrt x)))
+          (and (inexact? got)
+               (test-eqv? 1.0 (/ want got))))))
+    (test 1.511139943175573e176   (/ (expt 3 2001) (expt 2 2001)))
+    (test 2.1370746022826034e176  (/ (expt 3 2001) (expt 2 2000)))
+    (test 8.724570529756128e175   (/ (expt 3 2000) (expt 2 2001)))
+    (test 6.6175207962444435e-177 (/ (expt 2 2001) (expt 3 2001)))
+    (test 1.1461882239239027e-176 (/ (expt 2 2001) (expt 3 2000)))
+    (test 4.679293829667447e-177  (/ (expt 2 2000) (expt 3 2001))))
+
+  (pass-if (eqv? (/ (expt 2 1000)
+                    (expt 3 1000))
+                 (sqrt (/ (expt 2 2000)
+                          (expt 3 2000)))))
+
+  (pass-if (eqv? (/ (expt 3 1000)
+                    (expt 2 1000))
+                 (sqrt (/ (expt 3 2000)
+                          (expt 2 2000)))))
 
   (pass-if (eqv? +4i (sqrt -16)))
   (pass-if (eqv-loosely? +1.0e150i (sqrt #e-1e300)))


hooks/post-receive
-- 
GNU Guile