From ea995c6436d3489dbebf36d202f8cb7600b44dab Mon Sep 17 00:00:00 2001 From: Alex Shinn Date: Mon, 27 Jan 2014 23:24:40 +0900 Subject: [PATCH] Moving exact-sqrt into the core, with exact-integer-sqrt a variant that simply wraps in values. --- bignum.c | 11 ++++---- eval.c | 58 +++++++++++++++++++++++++++++++----------- include/chibi/bignum.h | 2 +- include/chibi/eval.h | 2 ++ lib/scheme/extras.scm | 10 ++------ opcodes.c | 1 + 6 files changed, 54 insertions(+), 30 deletions(-) diff --git a/bignum.c b/bignum.c index 30c318ed..c00e82cc 100644 --- a/bignum.c +++ b/bignum.c @@ -583,7 +583,8 @@ sexp sexp_bignum_expt (sexp ctx, sexp a, sexp b) { #define SEXP_MAX_ACCURATE_FLONUM_SQRT 1073741824.0 /* 2^30 */ -sexp sexp_bignum_sqrt (sexp ctx, sexp a) { /* Babylonian method */ +/* Babylonian method */ +sexp sexp_bignum_sqrt (sexp ctx, sexp a, sexp* rem_out) { sexp_gc_var4(res, rem, tmp, tmpa); if (! sexp_bignump(a)) return sexp_type_exception(ctx, NULL, SEXP_BIGNUM, a); sexp_gc_preserve4(ctx, res, rem, tmp, tmpa); @@ -594,10 +595,10 @@ sexp sexp_bignum_sqrt (sexp ctx, sexp a) { /* Babylonian method */ sexp_negate(a); } res = sexp_make_flonum(ctx, sexp_bignum_to_double(a)); - res = sexp_sqrt(ctx, NULL, 1, res); + res = sexp_inexact_sqrt(ctx, NULL, 1, res); if (sexp_flonump(res) && sexp_flonum_value(res) > SEXP_MAX_ACCURATE_FLONUM_SQRT) { - if (!isfinite(sexp_flonum_value(res))) + if (isinf(sexp_flonum_value(res))) res = sexp_make_flonum(ctx, 1e+154); res = sexp_double_to_bignum(ctx, sexp_flonum_value(res)); loop: /* until 0 <= a - res*res < 2*res + 1 */ @@ -615,9 +616,7 @@ sexp sexp_bignum_sqrt (sexp ctx, sexp a) { /* Babylonian method */ goto loop; } /* convert back to inexact if non-zero remainder */ - rem = sexp_bignum_normalize(rem); - if (rem != SEXP_ZERO) - res = sexp_make_flonum(ctx, sexp_fixnump(res) ? sexp_unbox_fixnum(res) : sexp_bignum_to_double(res)); + *rem_out = sexp_bignum_normalize(rem); } sexp_gc_release4(ctx); return sexp_bignum_normalize(res); diff --git a/eval.c b/eval.c index 4fc7191b..0bc18b3d 100644 --- a/eval.c +++ b/eval.c @@ -1452,44 +1452,33 @@ sexp sexp_log (sexp ctx, sexp self, sexp_sint_t n, sexp z) { return sexp_make_flonum(ctx, log(d)); } -sexp sexp_sqrt (sexp ctx, sexp self, sexp_sint_t n, sexp z) { -#if SEXP_USE_COMPLEX || SEXP_USE_BIGNUMS +sexp sexp_inexact_sqrt (sexp ctx, sexp self, sexp_sint_t n, sexp z) { +#if SEXP_USE_COMPLEX int negativep = 0; #endif double d, r; sexp_gc_var1(res); - sexp_gc_preserve1(ctx, 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_ratio(z) /* XXXX add ratio sqrt */ maybe_convert_complex(z, sexp_complex_sqrt) - else { - sexp_gc_release1(ctx); + else return sexp_type_exception(ctx, self, SEXP_NUMBER, z); - } #if SEXP_USE_COMPLEX if (d < 0) { negativep = 1; d = -d; } #endif + sexp_gc_preserve1(ctx, res); r = sqrt(d); if (sexp_fixnump(z) && (((sexp_uint_t)r*(sexp_uint_t)r)==abs(sexp_unbox_fixnum(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); @@ -1498,6 +1487,45 @@ sexp sexp_sqrt (sexp ctx, sexp self, sexp_sint_t n, sexp z) { return res; } +#if SEXP_USE_BIGNUMS +sexp sexp_exact_sqrt (sexp ctx, sexp self, sexp_sint_t n, sexp z) { + sexp_gc_var2(res, rem); + sexp_gc_preserve2(ctx, res, rem); + if (sexp_bignump(z)) { + res = sexp_bignum_sqrt(ctx, z, &rem); + res = sexp_cons(ctx, res, rem); + } else { + res = sexp_inexact_sqrt(ctx, self, n, z); + if (sexp_flonump(res)) { + res = sexp_bignum_normalize(sexp_double_to_bignum(ctx, trunc(sexp_flonum_value(res)))); + } + if (!sexp_exceptionp(res)) { + rem = sexp_mul(ctx, res, res); + rem = sexp_sub(ctx, z, rem); + res = sexp_cons(ctx, res, rem); + } + } + sexp_gc_release2(ctx); + return res; +} +#endif + +sexp sexp_sqrt (sexp ctx, sexp self, sexp_sint_t n, sexp z) { +#if SEXP_USE_BIGNUMS + sexp_gc_var2(res, rem); + if (sexp_bignump(z)) { + sexp_gc_preserve2(ctx, res, rem); + res = sexp_bignum_sqrt(ctx, z, &rem); + rem = sexp_bignum_normalize(rem); + if (rem != SEXP_ZERO) + res = sexp_make_flonum(ctx, sexp_fixnump(res) ? sexp_unbox_fixnum(res) : sexp_bignum_to_double(res)); + sexp_gc_release2(ctx); + return res; + } +#endif + return sexp_inexact_sqrt(ctx, self, n, z); +} + #endif /* SEXP_USE_MATH */ #if SEXP_USE_RATIOS || !SEXP_USE_FLONUMS diff --git a/include/chibi/bignum.h b/include/chibi/bignum.h index 1ae57703..1c8402c6 100644 --- a/include/chibi/bignum.h +++ b/include/chibi/bignum.h @@ -34,7 +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_bignum_sqrt (sexp ctx, sexp a, sexp* rem); 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/include/chibi/eval.h b/include/chibi/eval.h index 6796db84..6efc7a11 100644 --- a/include/chibi/eval.h +++ b/include/chibi/eval.h @@ -162,6 +162,8 @@ SEXP_API sexp sexp_asin(sexp ctx, sexp self, sexp_sint_t n, sexp z); SEXP_API sexp sexp_acos(sexp ctx, sexp self, sexp_sint_t n, sexp z); SEXP_API sexp sexp_atan(sexp ctx, sexp self, sexp_sint_t n, sexp z); SEXP_API sexp sexp_sqrt(sexp ctx, sexp self, sexp_sint_t n, sexp z); +SEXP_API sexp sexp_exact_sqrt(sexp ctx, sexp self, sexp_sint_t n, sexp z); +SEXP_API sexp sexp_inexact_sqrt(sexp ctx, sexp self, sexp_sint_t n, sexp z); SEXP_API sexp sexp_round(sexp ctx, sexp self, sexp_sint_t n, sexp x); SEXP_API sexp sexp_trunc(sexp ctx, sexp self, sexp_sint_t n, sexp x); SEXP_API sexp sexp_floor(sexp ctx, sexp self, sexp_sint_t n, sexp x); diff --git a/lib/scheme/extras.scm b/lib/scheme/extras.scm index c8d0fb00..b717ade2 100644 --- a/lib/scheme/extras.scm +++ b/lib/scheme/extras.scm @@ -57,14 +57,8 @@ (values (floor-quotient n m) (floor-remainder n m))) (define (exact-integer-sqrt x) - (let ((res (sqrt x))) - (if (exact? res) - (values res 0) - (let lp ((res (inexact->exact (truncate res)))) - (let ((rem (- x (* res res)))) - (if (negative? rem) - (lp (quotient (+ res (quotient x res)) 2)) - (values res rem))))))) + (let ((res (exact-sqrt x))) + (values (car res) (cdr res)))) ;; Adapted from Bawden's algorithm. (define (rationalize x e) diff --git a/opcodes.c b/opcodes.c index a752009b..d6a60f3a 100644 --- a/opcodes.c +++ b/opcodes.c @@ -203,6 +203,7 @@ _FN1(_I(SEXP_NUMBER), _I(SEXP_NUMBER), "asin", 0, sexp_asin), _FN1(_I(SEXP_NUMBER), _I(SEXP_NUMBER), "acos", 0, sexp_acos), _FN1(_I(SEXP_NUMBER), _I(SEXP_NUMBER), "atan1", 0, sexp_atan), _FN1(_I(SEXP_NUMBER), _I(SEXP_NUMBER), "sqrt", 0, sexp_sqrt), +_FN1(_I(SEXP_PAIR), _I(SEXP_NUMBER), "exact-sqrt", 0, sexp_exact_sqrt), _FN1(_I(SEXP_NUMBER), _I(SEXP_NUMBER), "round", 0, sexp_round), _FN1(_I(SEXP_NUMBER), _I(SEXP_NUMBER), "truncate", 0, sexp_trunc), _FN1(_I(SEXP_NUMBER), _I(SEXP_NUMBER), "floor", 0, sexp_floor),