diff options
author | Andy Wingo <wingo@pobox.com> | 2022-01-07 09:57:50 +0100 |
---|---|---|
committer | Andy Wingo <wingo@pobox.com> | 2022-01-13 09:37:17 +0100 |
commit | 124d8892274e7293f12ac4d1c4bf0053d6d3a51d (patch) | |
tree | c04a78940a6238118eac398b76ee12ad17688950 | |
parent | 63a18a6c1a2eae570c5de2b9eb866b20bc5e5e5e (diff) | |
download | guile-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.c | 73 | ||||
-rw-r--r-- | libguile/integers.h | 7 | ||||
-rw-r--r-- | 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); |