better inexact computation for ratios which overflow double (issue #671)

This commit is contained in:
Alex Shinn 2020-07-15 16:38:56 +09:00
parent 9104fcc44e
commit 983829cab1
7 changed files with 71 additions and 55 deletions

View file

@ -737,12 +737,25 @@ sexp sexp_bignum_sqrt (sexp ctx, sexp a, sexp* rem_out) {
#if SEXP_USE_RATIOS
double sexp_ratio_to_double (sexp rat) {
double sexp_ratio_to_double (sexp ctx, sexp rat) {
sexp_gc_var1(quot);
sexp num = sexp_ratio_numerator(rat), den = sexp_ratio_denominator(rat);
return (sexp_bignump(num) ? sexp_bignum_to_double(num)
double res = (sexp_bignump(num) ? sexp_bignum_to_double(num)
: sexp_fixnum_to_double(num))
/ (sexp_bignump(den) ? sexp_bignum_to_double(den)
: sexp_fixnum_to_double(den));
if (!isfinite(res)) {
sexp_gc_preserve1(ctx, quot);
if (sexp_unbox_fixnum(sexp_compare(ctx, sexp_ratio_numerator(rat), sexp_ratio_denominator(rat))) < 0) {
quot = sexp_quotient(ctx, sexp_ratio_denominator(rat), sexp_ratio_numerator(rat));
res = 1 / sexp_to_double(ctx, quot);
} else {
quot = sexp_quotient(ctx, sexp_ratio_numerator(rat), sexp_ratio_denominator(rat));
res = sexp_to_double(ctx, quot);
}
sexp_gc_release1(ctx);
}
return res;
}
sexp sexp_double_to_ratio (sexp ctx, double f) {
@ -895,7 +908,7 @@ sexp sexp_ratio_ceiling (sexp ctx, sexp a) {
#endif
double sexp_to_double (sexp x) {
double sexp_to_double (sexp ctx, sexp x) {
if (sexp_flonump(x))
return sexp_flonum_value(x);
else if (sexp_fixnump(x))
@ -904,7 +917,7 @@ double sexp_to_double (sexp x) {
return sexp_bignum_to_double(x);
#if SEXP_USE_RATIOS
else if (sexp_ratiop(x))
return sexp_ratio_to_double(x);
return sexp_ratio_to_double(ctx, x);
#endif
else
return 0.0;
@ -1001,7 +1014,7 @@ static sexp sexp_to_complex (sexp ctx, sexp x) {
} else if (sexp_ratiop(x)) {
sexp_gc_preserve1(ctx, tmp);
tmp = sexp_make_complex(ctx, SEXP_ZERO, SEXP_ZERO);
sexp_complex_real(tmp) = sexp_make_flonum(ctx, sexp_to_double(x));
sexp_complex_real(tmp) = sexp_make_flonum(ctx, sexp_to_double(ctx, x));
sexp_gc_release1(ctx);
return tmp;
#endif
@ -1011,8 +1024,8 @@ static sexp sexp_to_complex (sexp ctx, sexp x) {
}
sexp sexp_complex_exp (sexp ctx, sexp z) {
double e2x = exp(sexp_to_double(sexp_complex_real(z))),
y = sexp_to_double(sexp_complex_imag(z));
double e2x = exp(sexp_to_double(ctx, sexp_complex_real(z))),
y = sexp_to_double(ctx, sexp_complex_imag(z));
sexp_gc_var1(res);
sexp_gc_preserve1(ctx, res);
res = sexp_make_complex(ctx, SEXP_ZERO, SEXP_ZERO);
@ -1023,8 +1036,8 @@ sexp sexp_complex_exp (sexp ctx, sexp z) {
}
sexp sexp_complex_log (sexp ctx, sexp z) {
double x = sexp_to_double(sexp_complex_real(z)),
y = sexp_to_double(sexp_complex_imag(z));
double x = sexp_to_double(ctx, sexp_complex_real(z)),
y = sexp_to_double(ctx, sexp_complex_imag(z));
sexp_gc_var1(res);
sexp_gc_preserve1(ctx, res);
res = sexp_make_complex(ctx, SEXP_ZERO, SEXP_ZERO);
@ -1048,8 +1061,8 @@ sexp sexp_complex_expt (sexp ctx, sexp a, sexp b) {
#if SEXP_USE_MATH
sexp sexp_complex_sqrt (sexp ctx, sexp z) {
double x = sexp_to_double(sexp_complex_real(z)),
y = sexp_to_double(sexp_complex_imag(z)), r;
double x = sexp_to_double(ctx, sexp_complex_real(z)),
y = sexp_to_double(ctx, sexp_complex_imag(z)), r;
sexp_gc_var1(res);
sexp_gc_preserve1(ctx, res);
r = sqrt(x*x + y*y);
@ -1061,8 +1074,8 @@ sexp sexp_complex_sqrt (sexp ctx, sexp z) {
}
sexp sexp_complex_sin (sexp ctx, sexp z) {
double x = sexp_to_double(sexp_complex_real(z)),
y = sexp_to_double(sexp_complex_imag(z));
double x = sexp_to_double(ctx, sexp_complex_real(z)),
y = sexp_to_double(ctx, sexp_complex_imag(z));
sexp_gc_var1(res);
sexp_gc_preserve1(ctx, res);
res = sexp_make_complex(ctx, SEXP_ZERO, SEXP_ZERO);
@ -1073,8 +1086,8 @@ sexp sexp_complex_sin (sexp ctx, sexp z) {
}
sexp sexp_complex_cos (sexp ctx, sexp z) {
double x = sexp_to_double(sexp_complex_real(z)),
y = sexp_to_double(sexp_complex_imag(z));
double x = sexp_to_double(ctx, sexp_complex_real(z)),
y = sexp_to_double(ctx, sexp_complex_imag(z));
sexp_gc_var1(res);
sexp_gc_preserve1(ctx, res);
res = sexp_make_complex(ctx, SEXP_ZERO, SEXP_ZERO);
@ -1288,7 +1301,7 @@ sexp sexp_add (sexp ctx, sexp a, sexp b) {
break;
#if SEXP_USE_RATIOS
case SEXP_NUM_FLO_RAT:
r = sexp_make_flonum(ctx, sexp_flonum_value(a) + sexp_ratio_to_double(b));
r = sexp_make_flonum(ctx, sexp_flonum_value(a) + sexp_ratio_to_double(ctx, b));
break;
case SEXP_NUM_FIX_RAT:
case SEXP_NUM_BIG_RAT:
@ -1377,10 +1390,10 @@ sexp sexp_sub (sexp ctx, sexp a, sexp b) {
break;
#if SEXP_USE_RATIOS
case SEXP_NUM_FLO_RAT:
r = sexp_make_flonum(ctx, sexp_flonum_value(a) - sexp_ratio_to_double(b));
r = sexp_make_flonum(ctx, sexp_flonum_value(a) - sexp_ratio_to_double(ctx, b));
break;
case SEXP_NUM_RAT_FLO:
r = sexp_make_flonum(ctx, sexp_ratio_to_double(a) - sexp_flonum_value(b));
r = sexp_make_flonum(ctx, sexp_ratio_to_double(ctx, a) - sexp_flonum_value(b));
break;
case SEXP_NUM_RAT_FIX:
case SEXP_NUM_RAT_BIG:
@ -1409,10 +1422,10 @@ sexp sexp_sub (sexp ctx, sexp a, sexp b) {
#if SEXP_USE_COMPLEX
#if SEXP_USE_RATIOS
case SEXP_NUM_RAT_CPX:
a = tmp1 = sexp_make_flonum(ctx, sexp_ratio_to_double(a));
a = tmp1 = sexp_make_flonum(ctx, sexp_ratio_to_double(ctx, a));
goto complex_sub;
case SEXP_NUM_CPX_RAT:
b = tmp1 = sexp_make_flonum(ctx, sexp_ratio_to_double(b));
b = tmp1 = sexp_make_flonum(ctx, sexp_ratio_to_double(ctx, b));
/* ... FALLTHROUGH ... */
#endif
case SEXP_NUM_CPX_FLO:
@ -1488,7 +1501,7 @@ sexp sexp_mul (sexp ctx, sexp a, sexp b) {
break;
#if SEXP_USE_RATIOS
case SEXP_NUM_FLO_RAT:
r = sexp_make_flonum(ctx, sexp_flonum_value(a) * sexp_ratio_to_double(b));
r = sexp_make_flonum(ctx, sexp_flonum_value(a) * sexp_ratio_to_double(ctx, b));
break;
case SEXP_NUM_FIX_RAT:
case SEXP_NUM_BIG_RAT:
@ -1595,10 +1608,10 @@ sexp sexp_div (sexp ctx, sexp a, sexp b) {
break;
#if SEXP_USE_RATIOS
case SEXP_NUM_FLO_RAT:
r = sexp_make_flonum(ctx, sexp_flonum_value(a) / sexp_ratio_to_double(b));
r = sexp_make_flonum(ctx, sexp_flonum_value(a) / sexp_ratio_to_double(ctx, b));
break;
case SEXP_NUM_RAT_FLO:
r = sexp_make_flonum(ctx, sexp_ratio_to_double(a) / sexp_flonum_value(b));
r = sexp_make_flonum(ctx, sexp_ratio_to_double(ctx, a) / sexp_flonum_value(b));
break;
case SEXP_NUM_RAT_FIX:
case SEXP_NUM_RAT_BIG:
@ -1616,7 +1629,7 @@ sexp sexp_div (sexp ctx, sexp a, sexp b) {
#if SEXP_USE_COMPLEX
#if SEXP_USE_RATIOS
case SEXP_NUM_CPX_RAT:
b = tmp = sexp_make_flonum(ctx, sexp_ratio_to_double(b));
b = tmp = sexp_make_flonum(ctx, sexp_ratio_to_double(ctx, b));
/* ... FALLTHROUGH ... */
#endif
case SEXP_NUM_CPX_FLO:
@ -1627,7 +1640,7 @@ sexp sexp_div (sexp ctx, sexp a, sexp b) {
#if SEXP_USE_RATIOS
case SEXP_NUM_RAT_CPX:
if (sexp_ratiop(a))
a = tmp = sexp_make_flonum(ctx, sexp_ratio_to_double(a));
a = tmp = sexp_make_flonum(ctx, sexp_ratio_to_double(ctx, a));
/* ... FALLTHROUGH ... */
#endif
case SEXP_NUM_FLO_CPX:
@ -1878,7 +1891,7 @@ sexp sexp_compare (sexp ctx, sexp a, sexp b) {
} else if (isnan(f)) {
r = sexp_xtype_exception(ctx, NULL, "can't compare NaN", a);
} else {
g = sexp_ratio_to_double(b);
g = sexp_ratio_to_double(ctx, b);
r = sexp_make_fixnum(f < g ? -1 : f == g ? 0 : 1);
}
break;

20
eval.c
View file

@ -1486,10 +1486,10 @@ sexp sexp_register_optimization (sexp ctx, sexp self, sexp_sint_t n, sexp f, sex
#endif
#if SEXP_USE_RATIOS
#define maybe_convert_ratio(z) \
else if (sexp_ratiop(z)) d = sexp_ratio_to_double(z);
#define maybe_convert_ratio(ctx, z) \
else if (sexp_ratiop(z)) d = sexp_ratio_to_double(ctx, z);
#else
#define maybe_convert_ratio(z)
#define maybe_convert_ratio(ctx, z)
#endif
#if SEXP_USE_COMPLEX
@ -1507,7 +1507,7 @@ sexp sexp_register_optimization (sexp ctx, sexp self, sexp_sint_t n, sexp f, sex
d = sexp_flonum_value(z); \
else if (sexp_fixnump(z)) \
d = (double)sexp_unbox_fixnum(z); \
maybe_convert_ratio(z) \
maybe_convert_ratio(ctx, z) \
maybe_convert_bignum(z) \
maybe_convert_complex(z, f) \
else \
@ -1523,7 +1523,7 @@ sexp sexp_register_optimization (sexp ctx, sexp self, sexp_sint_t n, sexp f, sex
d = sexp_flonum_value(z); \
else if (sexp_fixnump(z)) \
d = (double)sexp_unbox_fixnum(z); \
maybe_convert_ratio(z) \
maybe_convert_ratio(ctx, z) \
maybe_convert_bignum(z) \
maybe_convert_complex(z, f) \
else \
@ -1586,7 +1586,7 @@ sexp sexp_log (sexp ctx, sexp self, sexp_sint_t n, sexp z) {
d = sexp_flonum_value(z);
else if (sexp_fixnump(z))
d = (double)sexp_unbox_fixnum(z);
maybe_convert_ratio(z)
maybe_convert_ratio(ctx, z)
maybe_convert_bignum(z)
else
return sexp_type_exception(ctx, self, SEXP_NUMBER, z);
@ -1613,7 +1613,7 @@ sexp sexp_inexact_sqrt (sexp ctx, sexp self, sexp_sint_t n, sexp z) {
d = sexp_flonum_value(z);
else if (sexp_fixnump(z))
d = (double)sexp_unbox_fixnum(z);
maybe_convert_ratio(z) /* XXXX add ratio sqrt */
maybe_convert_ratio(ctx, z) /* XXXX add ratio sqrt */
maybe_convert_complex(z, sexp_complex_sqrt)
else
return sexp_type_exception(ctx, self, SEXP_NUMBER, z);
@ -1730,7 +1730,7 @@ sexp sexp_expt_op (sexp ctx, sexp self, sexp_sint_t n, sexp x, sexp e) {
if (sexp_fixnump(e)) {
return sexp_generic_expt(ctx, x, sexp_unbox_fixnum(e));
} else {
x1 = sexp_ratio_to_double(x);
x1 = sexp_ratio_to_double(ctx, x);
}
}
#endif
@ -1742,7 +1742,7 @@ sexp sexp_expt_op (sexp ctx, sexp self, sexp_sint_t n, sexp x, sexp e) {
e1 = sexp_flonum_value(e);
#if SEXP_USE_RATIOS
else if (sexp_ratiop(e))
e1 = sexp_ratio_to_double(e);
e1 = sexp_ratio_to_double(ctx, e);
#endif
else
return sexp_type_exception(ctx, self, SEXP_FIXNUM, e);
@ -1804,7 +1804,7 @@ sexp sexp_exact_to_inexact (sexp ctx, sexp self, sexp_sint_t n, sexp i) {
#endif
#if SEXP_USE_RATIOS
else if (sexp_ratiop(i))
res = sexp_make_flonum(ctx, sexp_ratio_to_double(i));
res = sexp_make_flonum(ctx, sexp_ratio_to_double(ctx, i));
#endif
#if SEXP_USE_COMPLEX
else if (sexp_complexp(i)) {

View file

@ -1,5 +1,5 @@
/* bignum.h -- header for bignum utilities */
/* Copyright (c) 2009-2018 Alex Shinn. All rights reserved. */
/* Copyright (c) 2009-2020 Alex Shinn. All rights reserved. */
/* BSD-style license: http://synthcode.com/license.txt */
#ifndef SEXP_BIGNUM_H
@ -369,7 +369,7 @@ SEXP_API sexp_uint_t sexp_bignum_hi (sexp a);
SEXP_API sexp sexp_fixnum_to_bignum (sexp ctx, sexp a);
SEXP_API double sexp_bignum_to_double (sexp a);
SEXP_API sexp sexp_double_to_bignum (sexp ctx, double f);
SEXP_API double sexp_to_double (sexp x);
SEXP_API double sexp_to_double (sexp ctx, sexp x);
SEXP_API sexp sexp_bignum_fxadd (sexp ctx, sexp a, sexp_uint_t b);
SEXP_API sexp sexp_bignum_fxsub (sexp ctx, sexp a, sexp_uint_t b);
SEXP_API sexp sexp_bignum_fxmul (sexp ctx, sexp d, sexp a, sexp_uint_t b, int offset);
@ -389,7 +389,7 @@ 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 double sexp_ratio_to_double (sexp ctx, 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);
SEXP_API sexp sexp_ratio_round (sexp ctx, sexp a);

View file

@ -347,6 +347,9 @@
(let ((g 18446744073709551616/6148914691236517205))
(test 36893488147419103231/113427455640312821148309287786019553280
(- g (/ 9 g))))
(parameterize ((current-test-epsilon 1e-2))
(test 0.011459
(inexact 12017605191576564605523818479699032725649277767474470915174336170848732393575453450735711300329008209389459924841495413846951496975780183380894605321251029061627519460167494925074778514559449442501939557451941024018468063387066190930414864151964948746353911701589822769165039062500000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000/1048736796801977873293967427489899791695438298067922611215197702964142873750784860800028103241622682379099257583628328311596508052097172881170612593367688231078402334012314961284218456715145088673704198639789178608986339123909192734212369587774017431294143287678378871830938194004625162192783829580954040802482154686468390982830152927306297304746125796640118510663301998101921329795871645487301048602221374988871811684449297195367)))
(let ((r (/ (expt 2 61) 3)))
(test 0 (- r r))
(test 2305843009213693952/3 r))

View file

@ -98,9 +98,9 @@ sexp c64vector_ref(sexp ctx, float* uv, int i) {
sexp_gc_release3(ctx);
return res;
}
void c64vector_set(float* uv, int i, sexp v) {
uv[i*2] = sexp_to_double(sexp_real_part(v));
uv[i*2 + 1] = sexp_to_double(sexp_imag_part(v));
void c64vector_set(sexp ctx, float* uv, int i, sexp v) {
uv[i*2] = sexp_to_double(ctx, sexp_real_part(v));
uv[i*2 + 1] = sexp_to_double(ctx, sexp_imag_part(v));
}
sexp c128vector_ref(sexp ctx, double* uv, int i) {
@ -112,9 +112,9 @@ sexp c128vector_ref(sexp ctx, double* uv, int i) {
sexp_gc_release3(ctx);
return res;
}
void c128vector_set(double* uv, int i, sexp v) {
uv[i*2] = sexp_to_double(sexp_real_part(v));
uv[i*2 + 1] = sexp_to_double(sexp_imag_part(v));
void c128vector_set(sexp ctx, double* uv, int i, sexp v) {
uv[i*2] = sexp_to_double(ctx, sexp_real_part(v));
uv[i*2 + 1] = sexp_to_double(ctx, sexp_imag_part(v));
}
")
@ -164,7 +164,7 @@ void c128vector_set(double* uv, int i, sexp v) {
(define-c void (f64vector-set! "f64vector_set") (f64vector int double))
(define-c sexp c64vector-ref ((value ctx sexp) c64vector int))
(define-c void (c64vector-set! "c64vector_set") (c64vector int sexp))
(define-c void (c64vector-set! "c64vector_set") ((value ctx sexp) c64vector int sexp))
(define-c sexp c128vector-ref ((value ctx sexp) c128vector int))
(define-c void (c128vector-set! "c128vector_set") (c128vector int sexp))
(define-c void (c128vector-set! "c128vector_set") ((value ctx sexp) c128vector int sexp))

View file

@ -180,7 +180,7 @@ static void sexp_insert_timed (sexp ctx, sexp thread, sexp timeout) {
#endif
#if SEXP_USE_RATIOS
} else if (sexp_ratiop(timeout)) {
d = sexp_ratio_to_double(timeout);
d = sexp_ratio_to_double(ctx, timeout);
sexp_context_timeval(thread).tv_sec += trunc(d);
sexp_context_timeval(thread).tv_usec += (d-trunc(d))*1000000;
if (sexp_context_timeval(thread).tv_usec > 1000000) {

14
sexp.c
View file

@ -3073,22 +3073,22 @@ sexp sexp_list_to_uvector_op(sexp ctx, sexp self, sexp_sint_t n, sexp etype, sex
break;
#if SEXP_USE_FLONUMS
case SEXP_F32:
((float*)sexp_uvector_data(res))[i] = sexp_to_double(sexp_car(ls)); break;
((float*)sexp_uvector_data(res))[i] = sexp_to_double(ctx, sexp_car(ls)); break;
case SEXP_F64:
((double*)sexp_uvector_data(res))[i] = sexp_to_double(sexp_car(ls)); break;
((double*)sexp_uvector_data(res))[i] = sexp_to_double(ctx, sexp_car(ls)); break;
#endif
#if SEXP_USE_COMPLEX
case SEXP_C64:
((float*)sexp_uvector_data(res))[i*2] =
sexp_to_double(sexp_real_part(sexp_car(ls)));
sexp_to_double(ctx, sexp_real_part(sexp_car(ls)));
((float*)sexp_uvector_data(res))[i*2 + 1] =
sexp_to_double(sexp_imag_part(sexp_car(ls)));
sexp_to_double(ctx, sexp_imag_part(sexp_car(ls)));
break;
case SEXP_C128:
((double*)sexp_uvector_data(res))[i*2] =
sexp_to_double(sexp_real_part(sexp_car(ls)));
sexp_to_double(ctx, sexp_real_part(sexp_car(ls)));
((double*)sexp_uvector_data(res))[i*2 + 1] =
sexp_to_double(sexp_imag_part(sexp_car(ls)));
sexp_to_double(ctx, sexp_imag_part(sexp_car(ls)));
break;
#endif
}
@ -3293,7 +3293,7 @@ sexp sexp_read_raw (sexp ctx, sexp in, sexp *shares) {
res = sexp_make_flonum(ctx, sexp_unbox_fixnum(res));
#if SEXP_USE_RATIOS
else if (sexp_ratiop(res))
res = sexp_make_flonum(ctx, sexp_ratio_to_double(res));
res = sexp_make_flonum(ctx, sexp_ratio_to_double(ctx, res));
#endif
break;
case 'f': case 'F':