Fcalva/chibi-scheme
1
0
Fork 0
mirror of https://github.com/ashinn/chibi-scheme.git synced 2025-05-21 22:59:16 +02:00
Code Issues Projects Releases Packages Wiki Activity Actions

More accurate square roots for bignums - compute via iteration rather

Browse source 
than approximation via flonums for very large bignums.
This commit is contained in:
Alex Shinn 2012-12-22 18:06:22 +09:00
parent 48209fa1c1
commit 7ae254fc28
4 changed files with 52 additions and 1 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

38
bignum.c
View file

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

10
eval.c
View file

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

1
include/chibi/bignum.h
View file

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

4
tests/r7rs-tests.scm
View file

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