diff --git a/bignum.c b/bignum.c index 24431949..9dcbf936 100644 --- a/bignum.c +++ b/bignum.c @@ -477,6 +477,44 @@ sexp sexp_bignum_expt (sexp ctx, sexp a, sexp b) { return sexp_bignum_normalize(res); } +#define SEXP_MAX_ACCURATE_FLONUM_SQRT 1.15292150460685e18 /* 2^60 */ + +sexp sexp_bignum_sqrt (sexp ctx, sexp a) { /* Babylonian method */ + sexp_gc_var3(res, rem, tmp); + sexp_gc_preserve3(ctx, res, rem, tmp); + /* initial estimate via flonum, ignoring signs */ + res = sexp_make_flonum(ctx, sexp_bignum_to_double(a)); + if (sexp_negativep(res)) { sexp_negate(res); } + res = sexp_sqrt(ctx, NULL, 1, res); + if (sexp_flonump(res) && + sexp_flonum_value(res) > SEXP_MAX_ACCURATE_FLONUM_SQRT) { + res = sexp_double_to_bignum(ctx, sexp_flonum_value(res)); + loop: /* while abs(a - res*res) < res */ + rem = sexp_mul(ctx, res, res); + tmp = rem = sexp_sub(ctx, a, rem); + if (sexp_negativep(tmp)) { + tmp = sexp_copy_bignum(ctx, NULL, tmp, 0); + sexp_negate(tmp); + } + tmp = sexp_compare(ctx, tmp, res); + if (sexp_positivep(tmp)) { + tmp = sexp_quotient(ctx, a, res); + res = sexp_add(ctx, res, tmp); + res = sexp_quotient(ctx, res, SEXP_TWO); + goto loop; + } + /* add back in an approximate remainder if non-zero */ + rem = sexp_bignum_normalize(rem); + if (rem != SEXP_ZERO) { + rem = sexp_make_flonum(ctx, sexp_fixnump(rem) ? sexp_unbox_fixnum(rem) : sexp_bignum_to_double(rem)); + rem = sexp_div(ctx, rem, a); + res = sexp_add(ctx, res, rem); + } + } + sexp_gc_release3(ctx); + return sexp_bignum_normalize(res); +} + /************************ ratios ******************************/ #if SEXP_USE_RATIOS diff --git a/eval.c b/eval.c index 58a5ee1e..5548b75a 100644 --- a/eval.c +++ b/eval.c @@ -1351,11 +1351,16 @@ sexp sexp_sqrt (sexp ctx, sexp self, sexp_sint_t n, sexp z) { #endif double d, r; sexp_gc_var1(res); +#if SEXP_USE_BIGNUMS + if (sexp_bignump(z)) { + negativep = sexp_bignum_sign(z) < 0; + res = sexp_bignum_sqrt(ctx, z); + } else { +#endif if (sexp_flonump(z)) d = sexp_flonum_value(z); else if (sexp_fixnump(z)) d = (double)sexp_unbox_fixnum(z); - maybe_convert_bignum(z) /* XXXX add bignum sqrt */ maybe_convert_ratio(z) /* XXXX add ratio sqrt */ maybe_convert_complex(z, sexp_complex_sqrt) else @@ -1373,6 +1378,9 @@ sexp sexp_sqrt (sexp ctx, sexp self, sexp_sint_t n, sexp z) { res = sexp_make_fixnum(round(r)); else res = sexp_make_flonum(ctx, r); +#if SEXP_USE_BIGNUMS + } /* !sexp_bignump(z) */ +#endif #if SEXP_USE_COMPLEX if (negativep) res = sexp_make_complex(ctx, SEXP_ZERO, res); diff --git a/include/chibi/bignum.h b/include/chibi/bignum.h index d23262ee..34e6467c 100644 --- a/include/chibi/bignum.h +++ b/include/chibi/bignum.h @@ -34,6 +34,7 @@ SEXP_API sexp sexp_bignum_sub (sexp ctx, sexp dst, sexp a, sexp b); SEXP_API sexp sexp_bignum_mul (sexp ctx, sexp dst, sexp a, sexp b); SEXP_API sexp sexp_bignum_div (sexp ctx, sexp dst, sexp a, sexp b); SEXP_API sexp sexp_bignum_expt (sexp ctx, sexp n, sexp e); +SEXP_API sexp sexp_bignum_sqrt (sexp ctx, sexp a); SEXP_API sexp sexp_add (sexp ctx, sexp a, sexp b); SEXP_API sexp sexp_sub (sexp ctx, sexp a, sexp b); SEXP_API sexp sexp_mul (sexp ctx, sexp a, sexp b); diff --git a/tests/r7rs-tests.scm b/tests/r7rs-tests.scm index 234d55e6..2ea30edb 100644 --- a/tests/r7rs-tests.scm +++ b/tests/r7rs-tests.scm @@ -169,6 +169,10 @@ (let*-values (((root rem) (exact-integer-sqrt 32))) (test 35 (* root rem))) +(let*-values (((root rem) (exact-integer-sqrt (expt 2 140)))) + (test 0 rem) + (test (expt 2 140) (square root))) + (test '(x y x y) (let ((a 'a) (b 'b) (x 'x) (y 'y)) (let*-values (((a b) (values x y)) ((x y) (values a b)))