mirror/guile
mirror/guile
1
Fork 0
mirror of https://git.savannah.gnu.org/git/guile.git synced 2025-05-20 11:40:18 +02:00
Code Activity

Sqrt returns exact results when possible.

Browse source 
* 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.
This commit is contained in:
Mark H Weaver 2013-03-19 18:48:56 -04:00
parent c8248c8ed5
commit 4400266478
2 changed files with 103 additions and 6 deletions
Download patch file Download diff file Expand all files Collapse all files

69
libguile/numbers.c
View file

@ -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);

40
test-suite/tests/numbers.test
View file

@ -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)))