From b25e46b11bbff866962243a3d4e72e807314266d Mon Sep 17 00:00:00 2001 From: Jim Rees Date: Thu, 22 Mar 2018 22:19:39 -0400 Subject: [PATCH] Introduced a second version of sexp_double_to_ratio, named sexp_double_to_ratio_2, which converts without introducing round-off errors the way sexp_double_to_ratio does when it multiplies by 10. Changed sexp_inexact_to_exact to use this new function when a non-zero fractional part of the input exists. --- bignum.c | 35 +++++++++++++++++++++++++++++++++++ eval.c | 2 +- include/chibi/bignum.h | 1 + 3 files changed, 37 insertions(+), 1 deletion(-) diff --git a/bignum.c b/bignum.c index 757bc5de..b2e56a63 100644 --- a/bignum.c +++ b/bignum.c @@ -730,6 +730,41 @@ sexp sexp_double_to_ratio (sexp ctx, double f) { return res; } +// +// For conversion that does not introduce round-off error, +// no matter what FLT_RADIX is. +// +sexp sexp_double_to_ratio_2 (sexp ctx, double f) { + int sign,i; + sexp_gc_var3(res, whole, scale); + if (f == trunc(f)) + return sexp_bignum_normalize(sexp_double_to_bignum(ctx, f)); + sexp_gc_preserve3(ctx, res, whole, scale); + whole = sexp_double_to_bignum(ctx, trunc(f)); + res = sexp_fixnum_to_bignum(ctx, SEXP_ZERO); + scale = SEXP_ONE; + sign = (f < 0 ? -1 : 1); + f = fabs(f-trunc(f)); + while(f) { + res = sexp_bignum_fxmul(ctx, NULL, res, FLT_RADIX, 0); + scale = sexp_mul(ctx, scale, sexp_make_fixnum(FLT_RADIX)); + f *= FLT_RADIX; + i = trunc(f); + if (i) { + f -= i; + res = sexp_bignum_fxadd(ctx, res, i); + } + } + sexp_bignum_sign(res) = sign; + res = sexp_bignum_normalize(res); + scale = sexp_bignum_normalize(scale); + res = sexp_make_ratio(ctx, res, scale); + res = sexp_ratio_normalize(ctx, res, SEXP_FALSE); + res = sexp_add(ctx, res, whole); + sexp_gc_release3(ctx); + return res; +} + sexp sexp_ratio_add (sexp ctx, sexp a, sexp b) { sexp_gc_var3(res, num, den); sexp_gc_preserve3(ctx, res, num, den); diff --git a/eval.c b/eval.c index 6d6fccca..fada559f 100644 --- a/eval.c +++ b/eval.c @@ -1811,7 +1811,7 @@ sexp sexp_inexact_to_exact (sexp ctx, sexp self, sexp_sint_t n, sexp z) { res = sexp_xtype_exception(ctx, self, "exact: not a finite number", z); } else if (sexp_flonum_value(z) != trunc(sexp_flonum_value(z))) { #if SEXP_USE_RATIOS - res = sexp_double_to_ratio(ctx, sexp_flonum_value(z)); + res = sexp_double_to_ratio_2(ctx, sexp_flonum_value(z)); #else res = sexp_xtype_exception(ctx, self, "exact: not an integer", z); #endif diff --git a/include/chibi/bignum.h b/include/chibi/bignum.h index f554f55b..bd5b9791 100644 --- a/include/chibi/bignum.h +++ b/include/chibi/bignum.h @@ -44,6 +44,7 @@ SEXP_API sexp sexp_quotient (sexp ctx, sexp a, sexp b); SEXP_API sexp sexp_remainder (sexp ctx, sexp a, sexp b); #if SEXP_USE_RATIOS SEXP_API sexp sexp_double_to_ratio (sexp ctx, double f); +SEXP_API sexp sexp_double_to_ratio_2 (sexp ctx, double f); SEXP_API double sexp_ratio_to_double (sexp rat); SEXP_API sexp sexp_make_ratio (sexp ctx, sexp num, sexp den); SEXP_API sexp sexp_ratio_normalize (sexp ctx, sexp rat, sexp in);