}
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);
(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)))