summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorAndy Wingo <wingo@pobox.com>2022-01-07 09:57:50 +0100
committerAndy Wingo <wingo@pobox.com>2022-01-13 09:37:17 +0100
commit124d8892274e7293f12ac4d1c4bf0053d6d3a51d (patch)
treec04a78940a6238118eac398b76ee12ad17688950
parent63a18a6c1a2eae570c5de2b9eb866b20bc5e5e5e (diff)
downloadguile-124d8892274e7293f12ac4d1c4bf0053d6d3a51d.tar.gz
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.
-rw-r--r--libguile/integers.c73
-rw-r--r--libguile/integers.h7
-rw-r--r--libguile/numbers.c203
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);