summaryrefslogtreecommitdiff
path: root/libguile/numbers.c
diff options
context:
space:
mode:
Diffstat (limited to 'libguile/numbers.c')
-rw-r--r--libguile/numbers.c203
1 files changed, 50 insertions, 153 deletions
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);