[Guile-commits] 56/85: Refactor scm_sqrt in terms of integers.[ch]

[Andy Wingo]

wingo pushed a commit to branch main
in repository guile.

commit 124d8892274e7293f12ac4d1c4bf0053d6d3a51d
Author: Andy Wingo <wingo@pobox.com>
AuthorDate: Fri Jan 7 09:57:50 2022 +0100

    Refactor scm_sqrt in terms of integers.[ch]
    
    * libguile/integers.h:
    * libguile/integers.c (scm_is_integer_perfect_square_i):
    (scm_is_integer_perfect_square_z):
    (scm_integer_floor_sqrt_i):
    (scm_integer_floor_sqrt_z):
    (scm_integer_inexact_sqrt_i):
    (scm_integer_inexact_sqrt_z): New internal functions.
    * libguile/numbers.c (scm_sqrt): Reimplement.
---
 libguile/integers.c |  73 +++++++++++++++++++
 libguile/integers.h |   7 ++
 libguile/numbers.c  | 203 +++++++++++++---------------------------------------
 3 files changed, 130 insertions(+), 153 deletions(-)

diff --git a/libguile/integers.c b/libguile/integers.c
index d318fd775..2fde52625 100644
--- a/libguile/integers.c
+++ b/libguile/integers.c
@@ -3074,3 +3074,76 @@ scm_integer_exact_sqrt_z (struct scm_bignum *k, SCM *s, 
SCM *r)
   *s = take_mpz (zs);
   *r = take_mpz (zr);
 }
+
+int
+scm_is_integer_perfect_square_i (scm_t_inum k)
+{
+  if (k < 0)
+    return 0;
+  if (k == 0)
+    return 1;
+  mp_limb_t kk = k;
+  return mpn_perfect_square_p (&kk, 1);
+}
+
+int
+scm_is_integer_perfect_square_z (struct scm_bignum *k)
+{
+  mpz_t zk;
+  alias_bignum_to_mpz (k, zk);
+  int result = mpz_perfect_square_p (zk);
+  scm_remember_upto_here_1 (k);
+  return result;
+}
+
+SCM
+scm_integer_floor_sqrt_i (scm_t_inum k)
+{
+  if (k <= 0)
+    return SCM_INUM0;
+
+  mp_limb_t kk = k, ss;
+  mpn_sqrtrem (&ss, NULL, &kk, 1);
+  return SCM_I_MAKINUM (ss);
+}
+
+SCM
+scm_integer_floor_sqrt_z (struct scm_bignum *k)
+{
+  mpz_t zk, zs;
+  alias_bignum_to_mpz (k, zk);
+  mpz_init (zs);
+  mpz_sqrt (zs, zk);
+  scm_remember_upto_here_1 (k);
+  return take_mpz (zs);
+}
+
+double
+scm_integer_inexact_sqrt_i (scm_t_inum k)
+{
+  if (k < 0)
+    return -sqrt ((double) -k);
+  return sqrt ((double) k);
+}
+
+double
+scm_integer_inexact_sqrt_z (struct scm_bignum *k)
+{
+  mpz_t zk, zs;
+  alias_bignum_to_mpz (k, zk);
+  mpz_init (zs);
+
+  long expon;
+  double signif = bignum_frexp (k, &expon);
+  int negative = signif < 0;
+  if (negative)
+    signif = -signif;
+
+  if (expon & 1)
+    {
+      signif *= 2;
+      expon--;
+    }
+  double result = ldexp (sqrt (signif), expon / 2);
+  return negative ? -result : result;
+}
diff --git a/libguile/integers.h b/libguile/integers.h
index 1bba509ea..00706d553 100644
--- a/libguile/integers.h
+++ b/libguile/integers.h
@@ -216,6 +216,13 @@ SCM_INTERNAL void scm_integer_exact_sqrt_i (scm_t_inum k, 
SCM *s, SCM *r);
 SCM_INTERNAL void scm_integer_exact_sqrt_z (struct scm_bignum *k,
                                             SCM *s, SCM *r);
 
+SCM_INTERNAL int scm_is_integer_perfect_square_i (scm_t_inum k);
+SCM_INTERNAL int scm_is_integer_perfect_square_z (struct scm_bignum *k);
+SCM_INTERNAL SCM scm_integer_floor_sqrt_i (scm_t_inum k);
+SCM_INTERNAL SCM scm_integer_floor_sqrt_z (struct scm_bignum *k);
+SCM_INTERNAL double scm_integer_inexact_sqrt_i (scm_t_inum k);
+SCM_INTERNAL double scm_integer_inexact_sqrt_z (struct scm_bignum *k);
+
 
 
 #endif  /* SCM_INTEGERS_H */
diff --git a/libguile/numbers.c b/libguile/numbers.c
index 7567ced89..68534c02c 100644
--- a/libguile/numbers.c
+++ b/libguile/numbers.c
@@ -7390,62 +7390,6 @@ 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)))
-    {
-      if (SCM_I_INUM (k) > 0)
-        {
-          mp_limb_t kk = SCM_I_INUM (k);
-
-          result = mpn_perfect_square_p (&kk, 1);
-        }
-      else
-        result = (SCM_I_INUM (k) == 0);
-    }
-  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 non-negative integer. */
-static SCM
-exact_integer_floor_square_root (SCM k)
-{
-  if (SCM_LIKELY (SCM_I_INUMP (k)))
-    {
-      if (SCM_I_INUM (k) > 0)
-        {
-          mp_limb_t kk, ss, rr;
-
-          kk = SCM_I_INUM (k);
-          mpn_sqrtrem (&ss, &rr, &kk, 1);
-          return SCM_I_MAKINUM (ss);
-        }
-      else
-        return SCM_INUM0;
-    }
-  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),
        "Return the square root of @var{z}.  Of the two possible roots\n"
@@ -7461,7 +7405,35 @@ SCM_PRIMITIVE_GENERIC (scm_sqrt, "sqrt", 1, 0, 0,
        "@end example")
 #define FUNC_NAME s_scm_sqrt
 {
-  if (SCM_COMPLEXP (z))
+  if (SCM_I_INUMP (z))
+    {
+      scm_t_inum i = SCM_I_INUM (z);
+      if (scm_is_integer_perfect_square_i (i))
+        return scm_integer_floor_sqrt_i (i);
+      double root = scm_integer_inexact_sqrt_i (i);
+      return (root < 0)
+        ? scm_c_make_rectangular (0.0, -root)
+        : scm_i_from_double (root);
+    }
+  else if (SCM_BIGP (z))
+    {
+      struct scm_bignum *k = scm_bignum (z);
+      if (scm_is_integer_perfect_square_z (k))
+        return scm_integer_floor_sqrt_z (k);
+      double root = scm_integer_inexact_sqrt_z (k);
+      return (root < 0)
+        ? scm_c_make_rectangular (0.0, -root)
+        : scm_i_from_double (root);
+    }
+  else if (SCM_REALP (z))
+    {
+      double xx = SCM_REAL_VALUE (z);
+      if (xx < 0)
+        return scm_c_make_rectangular (0.0, sqrt (-xx));
+      else
+        return scm_i_from_double (sqrt (xx));
+    }
+  else if (SCM_COMPLEXP (z))
     {
 #if defined HAVE_COMPLEX_DOUBLE && defined HAVE_USABLE_CSQRT   \
       && defined SCM_COMPLEX_VALUE
@@ -7473,109 +7445,34 @@ SCM_PRIMITIVE_GENERIC (scm_sqrt, "sqrt", 1, 0, 0,
                                0.5 * atan2 (im, re));
 #endif
     }
-  else if (SCM_NUMBERP (z))
+  else if (SCM_FRACTIONP (z))
     {
-      if (SCM_I_INUMP (z))
-        {
-          scm_t_inum x = SCM_I_INUM (z);
+      SCM n = SCM_FRACTION_NUMERATOR (z);
+      SCM d = SCM_FRACTION_DENOMINATOR (z);
+      SCM nr = scm_sqrt (n);
+      SCM dr = scm_sqrt (d);
+      if (scm_is_exact_integer (nr) && scm_is_exact_integer (dr))
+        return scm_i_make_ratio_already_reduced (nr, dr);
 
-          if (SCM_LIKELY (x >= 0))
-            {
-              if (SCM_LIKELY (SCM_I_FIXNUM_BIT < DBL_MANT_DIG
-                              || x < (1L << (DBL_MANT_DIG - 1))))
-                {
-                  double root = sqrt (x);
-
-                  /* 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_i_from_double (root);
-                }
-              else
-                {
-                  mp_limb_t xx, root, rem;
-
-                  assert (x != 0);
-                  xx = x;
-                  if (mpn_perfect_square_p (&xx, 1))
-                    {
-                      mpn_sqrtrem (&root, &rem, &xx, 1);
-                      return SCM_I_MAKINUM (root);
-                    }
-                }
-            }
-        }
-      else if (SCM_BIGP (z))
-        {
-          if (mpz_perfect_square_p (SCM_I_BIG_MPZ (z)))
-            {
-              SCM root = scm_i_mkbig ();
+      double xx = scm_i_divide2double (n, d);
+      double abs_xx = fabs (xx);
+      long shift = 0;
 
-              mpz_sqrt (SCM_I_BIG_MPZ (root), SCM_I_BIG_MPZ (z));
-              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_i_from_double (ldexp (sqrt (signif), expon / 2));
-            }
-        }
-      else if (SCM_FRACTIONP (z))
+      if (abs_xx > DBL_MAX || abs_xx < DBL_MIN)
         {
-          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));
+          shift = (scm_to_long (scm_integer_length (n))
+                   - scm_to_long (scm_integer_length (d))) / 2;
+          if (shift > 0)
+            d = lsh (d, scm_from_long (2 * shift), FUNC_NAME);
           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 = lsh (d, scm_from_long (2 * shift), FUNC_NAME);
-                  else
-                    n = lsh (n, scm_from_long (-2 * shift), FUNC_NAME);
-                  xx = scm_i_divide2double (n, d);
-                }
-
-              if (xx < 0)
-                return scm_c_make_rectangular (0.0, ldexp (sqrt (-xx), shift));
-              else
-                return scm_i_from_double (ldexp (sqrt (xx), shift));
-            }
+            n = lsh (n, scm_from_long (-2 * shift), FUNC_NAME);
+          xx = scm_i_divide2double (n, d);
         }
 
-      /* 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_i_from_double (sqrt (xx));
-      }
+      if (xx < 0)
+        return scm_c_make_rectangular (0.0, ldexp (sqrt (-xx), shift));
+      else
+        return scm_i_from_double (ldexp (sqrt (xx), shift));
     }
   else
     return scm_wta_dispatch_1 (g_scm_sqrt, z, 1, s_scm_sqrt);