Fcalva/chibi-scheme
1
0
Fork 0
mirror of https://github.com/ashinn/chibi-scheme.git synced 2025-07-08 13:37:35 +02:00
Code Issues Projects Releases Packages Wiki Activity Actions 1

Fixing termination condition in certain bignum sqrts, handling negative inputs.

Browse source
This commit is contained in:
Alex Shinn 2013-01-03 23:49:31 +09:00
parent 1012f7e129
commit b42379539d
2 changed files with 17 additions and 17 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

30
bignum.c
View file

@ -480,38 +480,38 @@ sexp sexp_bignum_expt (sexp ctx, sexp a, sexp b) {
#define SEXP_MAX_ACCURATE_FLONUM_SQRT 1.15292150460685e18 /* 2^60 */ #define SEXP_MAX_ACCURATE_FLONUM_SQRT 1.15292150460685e18 /* 2^60 */
sexp sexp_bignum_sqrt (sexp ctx, sexp a) { /* Babylonian method */ sexp sexp_bignum_sqrt (sexp ctx, sexp a) { /* Babylonian method */
sexp_gc_var3(res, rem, tmp); sexp_gc_var4(res, rem, tmp, tmpa);
sexp_gc_preserve3(ctx, res, rem, tmp); sexp_gc_preserve4(ctx, res, rem, tmp, tmpa);
/* initial estimate via flonum, ignoring signs */ /* initial estimate via flonum, ignoring signs */
if (sexp_negativep(a)) {
a = tmpa = sexp_copy_bignum(ctx, NULL, a, 0);
sexp_negate(a);
}
res = sexp_make_flonum(ctx, sexp_bignum_to_double(a)); res = sexp_make_flonum(ctx, sexp_bignum_to_double(a));
if (sexp_negativep(res)) { sexp_negate(res); }
res = sexp_sqrt(ctx, NULL, 1, res); res = sexp_sqrt(ctx, NULL, 1, res);
if (sexp_flonump(res) && if (sexp_flonump(res) &&
sexp_flonum_value(res) > SEXP_MAX_ACCURATE_FLONUM_SQRT) { sexp_flonum_value(res) > SEXP_MAX_ACCURATE_FLONUM_SQRT) {
res = sexp_double_to_bignum(ctx, sexp_flonum_value(res)); res = sexp_double_to_bignum(ctx, sexp_flonum_value(res));
loop: /* while abs(a - res*res) < res */ loop: /* until 0 <= a - res*res < 2*res + 1 */
rem = sexp_mul(ctx, res, res); rem = sexp_mul(ctx, res, res);
tmp = rem = sexp_sub(ctx, a, rem); tmp = rem = sexp_sub(ctx, a, rem);
if (sexp_negativep(tmp)) { if (!sexp_positivep(tmp)) goto adjust;
tmp = sexp_copy_bignum(ctx, NULL, tmp, 0); tmp = sexp_sub(ctx, tmp, SEXP_ONE);
sexp_negate(tmp); tmp = sexp_quotient(ctx, tmp, SEXP_TWO);
}
tmp = sexp_compare(ctx, tmp, res); tmp = sexp_compare(ctx, tmp, res);
if (sexp_positivep(tmp)) { if (sexp_positivep(tmp)) {
adjust:
tmp = sexp_quotient(ctx, a, res); tmp = sexp_quotient(ctx, a, res);
res = sexp_add(ctx, res, tmp); res = sexp_add(ctx, res, tmp);
res = sexp_quotient(ctx, res, SEXP_TWO); res = sexp_quotient(ctx, res, SEXP_TWO);
goto loop; goto loop;
} }
/* add back in an approximate remainder if non-zero */ /* convert back to inexact if non-zero remainder */
rem = sexp_bignum_normalize(rem); rem = sexp_bignum_normalize(rem);
if (rem != SEXP_ZERO) { if (rem != SEXP_ZERO)
rem = sexp_make_flonum(ctx, sexp_fixnump(rem) ? sexp_unbox_fixnum(rem) : sexp_bignum_to_double(rem)); res = sexp_make_flonum(ctx, sexp_bignum_to_double(res));
rem = sexp_div(ctx, rem, a);
res = sexp_add(ctx, res, rem);
}
} }
sexp_gc_release3(ctx); sexp_gc_release4(ctx);
return sexp_bignum_normalize(res); return sexp_bignum_normalize(res);
} }

4
include/chibi/sexp.h
View file

@ -777,8 +777,8 @@ SEXP_API sexp sexp_make_unsigned_integer(sexp ctx, sexp_luint_t x);
#endif #endif
#define sexp_exact_negativep(x) (sexp_fixnump(x) ? (sexp_unbox_fixnum(x) < 0) \ #define sexp_exact_negativep(x) (sexp_fixnump(x) ? (sexp_unbox_fixnum(x) < 0) \
: (SEXP_USE_BIGNUMS && sexp_bignump(x)) \ : ((SEXP_USE_BIGNUMS && sexp_bignump(x)) \
&& (sexp_bignum_sign(x) < 0)) && (sexp_bignum_sign(x) < 0)))
#define sexp_negativep(x) (sexp_exact_negativep(x) || \ #define sexp_negativep(x) (sexp_exact_negativep(x) || \
(sexp_flonump(x) && sexp_flonum_value(x) < 0)) (sexp_flonump(x) && sexp_flonum_value(x) < 0))
#define sexp_positivep(x) (!(sexp_negativep(x))) #define sexp_positivep(x) (!(sexp_negativep(x)))