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.
This commit is contained in:
Jim Rees 2018-03-22 22:19:39 -04:00
parent 406aacf4dd
commit b25e46b11b
3 changed files with 37 additions and 1 deletions

View file

@ -730,6 +730,41 @@ sexp sexp_double_to_ratio (sexp ctx, double f) {
return res; 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 sexp_ratio_add (sexp ctx, sexp a, sexp b) {
sexp_gc_var3(res, num, den); sexp_gc_var3(res, num, den);
sexp_gc_preserve3(ctx, res, num, den); sexp_gc_preserve3(ctx, res, num, den);

2
eval.c
View file

@ -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); res = sexp_xtype_exception(ctx, self, "exact: not a finite number", z);
} else if (sexp_flonum_value(z) != trunc(sexp_flonum_value(z))) { } else if (sexp_flonum_value(z) != trunc(sexp_flonum_value(z))) {
#if SEXP_USE_RATIOS #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 #else
res = sexp_xtype_exception(ctx, self, "exact: not an integer", z); res = sexp_xtype_exception(ctx, self, "exact: not an integer", z);
#endif #endif

View file

@ -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); SEXP_API sexp sexp_remainder (sexp ctx, sexp a, sexp b);
#if SEXP_USE_RATIOS #if SEXP_USE_RATIOS
SEXP_API sexp sexp_double_to_ratio (sexp ctx, double f); 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 double sexp_ratio_to_double (sexp rat);
SEXP_API sexp sexp_make_ratio (sexp ctx, sexp num, sexp den); SEXP_API sexp sexp_make_ratio (sexp ctx, sexp num, sexp den);
SEXP_API sexp sexp_ratio_normalize (sexp ctx, sexp rat, sexp in); SEXP_API sexp sexp_ratio_normalize (sexp ctx, sexp rat, sexp in);