From 7ae254fc28cbe17c9fde6475dccf2c05819149de Mon Sep 17 00:00:00 2001
From: Alex Shinn <ashinn@users.noreply.github.com>
Date: Sat, 22 Dec 2012 18:06:22 +0900
Subject: [PATCH] More accurate square roots for bignums - compute via
iteration rather than approximation via flonums for very large bignums.
---
bignum.c | 38 ++++++++++++++++++++++++++++++++++++++
eval.c | 10 +++++++++-
include/chibi/bignum.h | 1 +
tests/r7rs-tests.scm | 4 ++++
4 files changed, 52 insertions(+), 1 deletion(-)
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)))