From 4400266478b4a477c6747f9eed38f7c6021491d8 Mon Sep 17 00:00:00 2001
From: Mark H Weaver <mhw@netris.org>
Date: Tue, 19 Mar 2013 18:48:56 -0400
Subject: [PATCH] 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.
---
 libguile/numbers.c            | 69 ++++++++++++++++++++++++++++++++---
 test-suite/tests/numbers.test | 40 +++++++++++++++++++-
 2 files changed, 103 insertions(+), 6 deletions(-)

diff --git a/libguile/numbers.c b/libguile/numbers.c
index a490f5d74..9725fe443 100644
--- a/libguile/numbers.c
+++ b/libguile/numbers.c
@@ -9988,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 be2e31732..a52e79a01 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)))
-- 
2.20.1