From 3afa3e1a421463b0e0c4e032a70b96e0feabda3f Mon Sep 17 00:00:00 2001 From: Marius Gerbershagen Date: Sun, 16 Feb 2020 20:57:59 +0100 Subject: [PATCH] ieee_fp: move fetestexcept checks directly after floating point calculations The C standard allows compilers to ignore side effects of floating point operations unless the STDC FENV_ACCESS pragma is in effect. This could lead to spurious floating point exceptions being signaled because the C compiler optimized a calculation that would ordinarily not set fpe bits to one that did (observed with clang for a cast double to cl_fixnum). To solve this, we resurrect the ECL_MATHERR macro which was removed in commit cb03494a6d7f33bb7b8d0f40da649b25934ba300. We handle fpe bits only around calculations which are known to be "unsafe" in the sense that they can produce floating point exceptions. Before the calculation, all bits are unset and afterwards, we test for exceptions. The disadvantage of this method is that optimizations by the C compiler resulting in unboxed float arithmetic might lead to missing exceptions, because our native C compiler does not insert ECL_MATHERR statements around these calculations. --- src/c/number.d | 10 ++-- src/c/numbers/atan.d | 32 +++++++------ src/c/numbers/ceiling.d | 10 ++++ src/c/numbers/divide.d | 90 +++++++++++++++++++++++------------ src/c/numbers/expt.d | 28 ++++++++--- src/c/numbers/floor.d | 10 ++++ src/c/numbers/minus.d | 94 +++++++++++++++++++++++------------- src/c/numbers/one_minus.d | 5 ++ src/c/numbers/one_plus.d | 5 ++ src/c/numbers/plus.d | 94 +++++++++++++++++++++++------------- src/c/numbers/round.d | 6 +++ src/c/numbers/times.d | 97 +++++++++++++++++++++++++------------- src/c/numbers/truncate.d | 5 ++ src/h/impl/math_dispatch.h | 2 + src/h/impl/math_fenv.h | 17 +++++++ 15 files changed, 352 insertions(+), 153 deletions(-) diff --git a/src/c/number.d b/src/c/number.d index fdedceee2..ba28cfe57 100644 --- a/src/c/number.d +++ b/src/c/number.d @@ -38,13 +38,11 @@ # define DO_DETECT_FPE2(f1,f2) DO_DETECT_FPE(f1) # else /* - * We need explicit checks for floating point exception bits being set + * We need explicit checks for floating point exception bits being + * set; ECL_MATHERR_TEST handles this for us, so nothing to do here. */ -# define DO_DETECT_FPE(f) do { \ - int status = fetestexcept(ecl_process_env()->trap_fpe_bits); \ - unlikely_if (status) ecl_deliver_fpe(status); \ - } while (0) -# define DO_DETECT_FPE2(f1,f2) DO_DETECT_FPE(f1) +# define DO_DETECT_FPE(f) +# define DO_DETECT_FPE2(f1,f2) # endif #else /* diff --git a/src/c/numbers/atan.d b/src/c/numbers/atan.d index e4c3b3269..e63bc8a5e 100644 --- a/src/c/numbers/atan.d +++ b/src/c/numbers/atan.d @@ -24,23 +24,27 @@ cl_object ecl_atan2(cl_object y, cl_object x) { cl_object output; - int tx = ecl_t_of(x); - int ty = ecl_t_of(y); - if (tx < ty) - tx = ty; - if (tx == t_longfloat) { - long double d = atan2l(ecl_to_long_double(y), ecl_to_long_double(x)); - output = ecl_make_long_float(d); - } else { - double dx = ecl_to_double(x); - double dy = ecl_to_double(y); - double dz = atan2(dy, dx); - if (tx == t_doublefloat) { - output = ecl_make_double_float(dz); + ECL_MATHERR_CLEAR; + { + int tx = ecl_t_of(x); + int ty = ecl_t_of(y); + if (tx < ty) + tx = ty; + if (tx == t_longfloat) { + long double d = atan2l(ecl_to_long_double(y), ecl_to_long_double(x)); + output = ecl_make_long_float(d); } else { - output = ecl_make_single_float(dz); + double dx = ecl_to_double(x); + double dy = ecl_to_double(y); + double dz = atan2(dy, dx); + if (tx == t_doublefloat) { + output = ecl_make_double_float(dz); + } else { + output = ecl_make_single_float(dz); + } } } + ECL_MATHERR_TEST; return output; } diff --git a/src/c/numbers/ceiling.d b/src/c/numbers/ceiling.d index 4599cca33..2cea13f36 100644 --- a/src/c/numbers/ceiling.d +++ b/src/c/numbers/ceiling.d @@ -18,6 +18,8 @@ #endif #include +#pragma STDC FENV_ACCESS ON + @(defun ceiling (x &optional (y OBJNULL)) @ if (narg == 1) @@ -30,6 +32,8 @@ cl_object ecl_ceiling1(cl_object x) { cl_object v0, v1; + ECL_MATHERR_CLEAR; + switch (ecl_t_of(x)) { case t_fixnum: case t_bignum: @@ -66,6 +70,8 @@ ecl_ceiling1(cl_object x) default: FEwrong_type_nth_arg(@[ceiling],1,x,@[real]); } + + ECL_MATHERR_TEST; @(return v0 v1); } @@ -75,6 +81,8 @@ ecl_ceiling2(cl_object x, cl_object y) const cl_env_ptr the_env = ecl_process_env(); cl_object v0, v1; cl_type ty; + ECL_MATHERR_CLEAR; + v0 = v1 = ECL_NIL; ty = ecl_t_of(y); if (ecl_unlikely(!ECL_REAL_TYPE_P(ty))) { @@ -221,5 +229,7 @@ ecl_ceiling2(cl_object x, cl_object y) default: FEwrong_type_nth_arg(@[ceiling], 1, x, @[real]); } + + ECL_MATHERR_TEST; ecl_return2(the_env, v0, v1); } diff --git a/src/c/numbers/divide.d b/src/c/numbers/divide.d index aadff99f1..2d8393963 100644 --- a/src/c/numbers/divide.d +++ b/src/c/numbers/divide.d @@ -12,9 +12,12 @@ * */ +#define ECL_INCLUDE_MATH_H #include #include +#pragma STDC FENV_ACCESS ON + @(defun / (num &rest nums) @ /* INV: type check is in ecl_divide() */ @@ -41,6 +44,9 @@ complex_divide(cl_object ar, cl_object ai, cl_object br, cl_object bi) cl_object ecl_divide(cl_object x, cl_object y) { + cl_object ret; + ECL_MATHERR_CLEAR; + MATH_DISPATCH2_BEGIN(x,y) { CASE_FIXNUM_FIXNUM; @@ -58,18 +64,22 @@ ecl_divide(cl_object x, cl_object y) y->ratio.num); } CASE_FIXNUM_SINGLE_FLOAT { - return ecl_make_single_float(ecl_fixnum(x) / ecl_single_float(y)); + ret = ecl_make_single_float(ecl_fixnum(x) / ecl_single_float(y)); + break; } CASE_FIXNUM_DOUBLE_FLOAT { - return ecl_make_double_float(ecl_fixnum(x) / ecl_double_float(y)); + ret = ecl_make_double_float(ecl_fixnum(x) / ecl_double_float(y)); + break; } CASE_BIGNUM_SINGLE_FLOAT; CASE_RATIO_SINGLE_FLOAT { - return ecl_make_single_float(ecl_to_float(x) / ecl_single_float(y)); + ret = ecl_make_single_float(ecl_to_float(x) / ecl_single_float(y)); + break; } CASE_BIGNUM_DOUBLE_FLOAT; CASE_RATIO_DOUBLE_FLOAT { - return ecl_make_double_float(ecl_to_double(x) / ecl_double_float(y)); + ret = ecl_make_double_float(ecl_to_double(x) / ecl_double_float(y)); + break; } CASE_RATIO_FIXNUM { if (y == ecl_make_fixnum(0)) { @@ -86,59 +96,76 @@ ecl_divide(cl_object x, cl_object y) return ecl_make_ratio(num, den); } CASE_SINGLE_FLOAT_FIXNUM { - return ecl_make_single_float(ecl_single_float(x) / ecl_fixnum(y)); + ret = ecl_make_single_float(ecl_single_float(x) / ecl_fixnum(y)); + break; } CASE_SINGLE_FLOAT_BIGNUM; CASE_SINGLE_FLOAT_RATIO { - return ecl_make_single_float(ecl_single_float(x) / ecl_to_float(y)); + ret = ecl_make_single_float(ecl_single_float(x) / ecl_to_float(y)); + break; } CASE_SINGLE_FLOAT_SINGLE_FLOAT { - return ecl_make_single_float(ecl_single_float(x) / ecl_single_float(y)); + ret = ecl_make_single_float(ecl_single_float(x) / ecl_single_float(y)); + break; } CASE_SINGLE_FLOAT_DOUBLE_FLOAT { - return ecl_make_double_float(ecl_single_float(x) / ecl_double_float(y)); + ret = ecl_make_double_float(ecl_single_float(x) / ecl_double_float(y)); + break; } CASE_DOUBLE_FLOAT_FIXNUM { - return ecl_make_double_float(ecl_double_float(x) / ecl_fixnum(y)); + ret = ecl_make_double_float(ecl_double_float(x) / ecl_fixnum(y)); + break; } CASE_DOUBLE_FLOAT_BIGNUM; CASE_DOUBLE_FLOAT_RATIO { - return ecl_make_double_float(ecl_double_float(x) / ecl_to_double(y)); + ret = ecl_make_double_float(ecl_double_float(x) / ecl_to_double(y)); + break; } CASE_DOUBLE_FLOAT_SINGLE_FLOAT { - return ecl_make_double_float(ecl_double_float(x) / ecl_single_float(y)); + ret = ecl_make_double_float(ecl_double_float(x) / ecl_single_float(y)); + break; } CASE_DOUBLE_FLOAT_DOUBLE_FLOAT { - return ecl_make_double_float(ecl_double_float(x) / ecl_double_float(y)); + ret = ecl_make_double_float(ecl_double_float(x) / ecl_double_float(y)); + break; } CASE_FIXNUM_LONG_FLOAT { - return ecl_make_long_float(ecl_fixnum(x) / ecl_long_float(y)); + ret = ecl_make_long_float(ecl_fixnum(x) / ecl_long_float(y)); + break; } CASE_BIGNUM_LONG_FLOAT; CASE_RATIO_LONG_FLOAT { - return ecl_make_long_float(ecl_to_long_double(x) / ecl_long_float(y)); + ret = ecl_make_long_float(ecl_to_long_double(x) / ecl_long_float(y)); + break; } CASE_SINGLE_FLOAT_LONG_FLOAT { - return ecl_make_long_float(ecl_single_float(x) / ecl_long_float(y)); + ret = ecl_make_long_float(ecl_single_float(x) / ecl_long_float(y)); + break; } CASE_DOUBLE_FLOAT_LONG_FLOAT { - return ecl_make_long_float(ecl_double_float(x) / ecl_long_float(y)); + ret = ecl_make_long_float(ecl_double_float(x) / ecl_long_float(y)); + break; } CASE_LONG_FLOAT_FIXNUM { - return ecl_make_long_float(ecl_long_float(x) / ecl_fixnum(y)); + ret = ecl_make_long_float(ecl_long_float(x) / ecl_fixnum(y)); + break; } CASE_LONG_FLOAT_BIGNUM; CASE_LONG_FLOAT_RATIO { - return ecl_make_long_float(ecl_long_float(x) / ecl_to_long_double(y)); + ret = ecl_make_long_float(ecl_long_float(x) / ecl_to_long_double(y)); + break; } CASE_LONG_FLOAT_SINGLE_FLOAT { - return ecl_make_long_float(ecl_long_float(x) / ecl_single_float(y)); + ret = ecl_make_long_float(ecl_long_float(x) / ecl_single_float(y)); + break; } CASE_LONG_FLOAT_DOUBLE_FLOAT { - return ecl_make_long_float(ecl_long_float(x) / ecl_double_float(y)); + ret = ecl_make_long_float(ecl_long_float(x) / ecl_double_float(y)); + break; } CASE_LONG_FLOAT_LONG_FLOAT { - return ecl_make_long_float(ecl_long_float(x) / ecl_long_float(y)); + ret = ecl_make_long_float(ecl_long_float(x) / ecl_long_float(y)); + break; } CASE_LONG_FLOAT_COMPLEX { goto COMPLEX_Y; @@ -180,9 +207,9 @@ ecl_divide(cl_object x, cl_object y) CASE_SINGLE_FLOAT_CSFLOAT; CASE_COMPLEX_CSFLOAT; CASE_CSFLOAT_CSFLOAT { - cl_object aux = ecl_alloc_object(t_csfloat); - ecl_csfloat(aux) = ecl_to_csfloat(x) / ecl_to_csfloat(y); - return aux; + ret = ecl_alloc_object(t_csfloat); + ecl_csfloat(ret) = ecl_to_csfloat(x) / ecl_to_csfloat(y); + break; } /* upgraded type cdfloat */ CASE_CSFLOAT_DOUBLE_FLOAT; @@ -204,9 +231,9 @@ ecl_divide(cl_object x, cl_object y) CASE_COMPLEX_CDFLOAT; CASE_CSFLOAT_CDFLOAT; CASE_CDFLOAT_CDFLOAT { - cl_object aux = ecl_alloc_object(t_cdfloat); - ecl_cdfloat(aux) = ecl_to_cdfloat(x) / ecl_to_cdfloat(y); - return aux; + ret = ecl_alloc_object(t_cdfloat); + ecl_cdfloat(ret) = ecl_to_cdfloat(x) / ecl_to_cdfloat(y); + break; } /* upgraded type clfloat */ CASE_CSFLOAT_LONG_FLOAT; @@ -234,12 +261,15 @@ ecl_divide(cl_object x, cl_object y) CASE_CSFLOAT_CLFLOAT; CASE_CDFLOAT_CLFLOAT; CASE_CLFLOAT_CLFLOAT { - cl_object aux = ecl_alloc_object(t_clfloat); - ecl_clfloat(aux) = ecl_to_clfloat(x) / ecl_to_clfloat(y); - return aux; + ret = ecl_alloc_object(t_clfloat); + ecl_clfloat(ret) = ecl_to_clfloat(x) / ecl_to_clfloat(y); + break; } #endif CASE_UNKNOWN(@[/],x,y,@[number]); } MATH_DISPATCH2_END; + + ECL_MATHERR_TEST; + return ret; } diff --git a/src/c/numbers/expt.d b/src/c/numbers/expt.d index cb83ada17..66b16276e 100644 --- a/src/c/numbers/expt.d +++ b/src/c/numbers/expt.d @@ -119,17 +119,25 @@ ecl_expt_float(cl_object x, cl_object y) { cl_type tx = ecl_t_of(x), ty = ecl_t_of(y); + cl_object ret; + + ECL_MATHERR_CLEAR; switch((ty > tx) ? ty : tx) { case t_longfloat: - return ecl_make_long_float + ret = ecl_make_long_float (powl(ecl_to_long_double(x), ecl_to_long_double(y))); + break; case t_doublefloat: - return ecl_make_double_float + ret = ecl_make_double_float (pow(ecl_to_double(x), ecl_to_double(y))); + break; default: - return ecl_make_single_float + ret = ecl_make_single_float (powf(ecl_to_float(x), ecl_to_float(y))); + break; } + ECL_MATHERR_TEST; + return ret; } #ifdef ECL_COMPLEX_FLOAT @@ -138,19 +146,27 @@ ecl_expt_complex_float(cl_object x, cl_object y) { cl_type tx = ecl_t_of(x), ty = ecl_t_of(y); + cl_object ret; + + ECL_MATHERR_CLEAR; switch ((ty > tx)? ty : tx) { case t_clfloat: case t_longfloat: - return ecl_make_clfloat + ret = ecl_make_clfloat (cpowl(ecl_to_clfloat(x), ecl_to_clfloat(y))); + break; case t_cdfloat: case t_doublefloat: - return ecl_make_cdfloat + ret = ecl_make_cdfloat (cpow (ecl_to_cdfloat(x), ecl_to_cdfloat(y))); + break; default: - return ecl_make_csfloat + ret = ecl_make_csfloat (cpowf(ecl_to_csfloat(x), ecl_to_csfloat(y))); + break; } + ECL_MATHERR_TEST; + return ret; } #endif diff --git a/src/c/numbers/floor.d b/src/c/numbers/floor.d index 7648473cf..57d7ba035 100644 --- a/src/c/numbers/floor.d +++ b/src/c/numbers/floor.d @@ -21,6 +21,8 @@ #include #include +#pragma STDC FENV_ACCESS ON + @(defun floor (x &optional (y OBJNULL)) @ if (narg == 1) @@ -34,6 +36,8 @@ ecl_floor1(cl_object x) { const cl_env_ptr the_env = ecl_process_env(); cl_object v0, v1; + ECL_MATHERR_CLEAR; + switch (ecl_t_of(x)) { case t_fixnum: case t_bignum: @@ -68,6 +72,8 @@ ecl_floor1(cl_object x) default: FEwrong_type_nth_arg(@[floor],1,x,@[real]); } + + ECL_MATHERR_TEST; ecl_return2(the_env, v0, v1); } @@ -76,6 +82,8 @@ ecl_floor2(cl_object x, cl_object y) { const cl_env_ptr the_env = ecl_process_env(); cl_object v0, v1; + ECL_MATHERR_CLEAR; + MATH_DISPATCH2_BEGIN(x,y) { CASE_FIXNUM_FIXNUM { @@ -235,6 +243,8 @@ ecl_floor2(cl_object x, cl_object y) } } MATH_DISPATCH2_END; + + ECL_MATHERR_TEST; ecl_return2(the_env, v0, v1); } diff --git a/src/c/numbers/minus.d b/src/c/numbers/minus.d index 1d8b4d051..11e276572 100644 --- a/src/c/numbers/minus.d +++ b/src/c/numbers/minus.d @@ -12,10 +12,12 @@ * */ - +#define ECL_INCLUDE_MATH_H #include #include +#pragma STDC FENV_ACCESS ON + @(defun - (num &rest nums) cl_object diff; @ @@ -31,6 +33,9 @@ cl_object ecl_minus(cl_object x, cl_object y) { + cl_object ret; + ECL_MATHERR_CLEAR; + MATH_DISPATCH2_BEGIN(x,y) { CASE_FIXNUM_FIXNUM { @@ -46,10 +51,12 @@ ecl_minus(cl_object x, cl_object y) return ecl_make_ratio(z, y->ratio.den); } CASE_FIXNUM_SINGLE_FLOAT { - return ecl_make_single_float(ecl_fixnum(x) - ecl_single_float(y)); + ret = ecl_make_single_float(ecl_fixnum(x) - ecl_single_float(y)); + break; } CASE_FIXNUM_DOUBLE_FLOAT { - return ecl_make_double_float(ecl_fixnum(x) - ecl_double_float(y)); + ret = ecl_make_double_float(ecl_fixnum(x) - ecl_double_float(y)); + break; } CASE_BIGNUM_FIXNUM { return _ecl_big_plus_fix(x, -ecl_fixnum(y)); @@ -59,11 +66,13 @@ ecl_minus(cl_object x, cl_object y) } CASE_BIGNUM_SINGLE_FLOAT; CASE_RATIO_SINGLE_FLOAT { - return ecl_make_single_float(ecl_to_float(x) - ecl_single_float(y)); + ret = ecl_make_single_float(ecl_to_float(x) - ecl_single_float(y)); + break; } CASE_BIGNUM_DOUBLE_FLOAT; CASE_RATIO_DOUBLE_FLOAT { - return ecl_make_double_float(ecl_to_double(x) - ecl_double_float(y)); + ret = ecl_make_double_float(ecl_to_double(x) - ecl_double_float(y)); + break; } CASE_RATIO_FIXNUM; /* fallthrough */ @@ -80,61 +89,79 @@ ecl_minus(cl_object x, cl_object y) return ecl_make_ratio(z, z1); } CASE_SINGLE_FLOAT_FIXNUM { - return ecl_make_single_float(ecl_single_float(x) - ecl_fixnum(y)); + ret = ecl_make_single_float(ecl_single_float(x) - ecl_fixnum(y)); + break; } CASE_SINGLE_FLOAT_BIGNUM; CASE_SINGLE_FLOAT_RATIO { - return ecl_make_single_float(ecl_single_float(x) - ecl_to_float(y)); + ret = ecl_make_single_float(ecl_single_float(x) - ecl_to_float(y)); + break; } CASE_SINGLE_FLOAT_SINGLE_FLOAT { - return ecl_make_single_float(ecl_single_float(x) - ecl_single_float(y)); + ret = ecl_make_single_float(ecl_single_float(x) - ecl_single_float(y)); + break; } CASE_SINGLE_FLOAT_DOUBLE_FLOAT { - return ecl_make_double_float(ecl_single_float(x) - ecl_double_float(y)); + ret = ecl_make_double_float(ecl_single_float(x) - ecl_double_float(y)); + break; } CASE_DOUBLE_FLOAT_FIXNUM { - return ecl_make_double_float(ecl_double_float(x) - ecl_fixnum(y)); + ret = ecl_make_double_float(ecl_double_float(x) - ecl_fixnum(y)); + break; } CASE_DOUBLE_FLOAT_BIGNUM; CASE_DOUBLE_FLOAT_RATIO { - return ecl_make_double_float(ecl_double_float(x) - ecl_to_double(y)); + ret = ecl_make_double_float(ecl_double_float(x) - ecl_to_double(y)); + break; } CASE_DOUBLE_FLOAT_SINGLE_FLOAT { - return ecl_make_double_float(ecl_double_float(x) - ecl_single_float(y)); + ret = ecl_make_double_float(ecl_double_float(x) - ecl_single_float(y)); + break; } CASE_DOUBLE_FLOAT_DOUBLE_FLOAT { - return ecl_make_double_float(ecl_double_float(x) - ecl_double_float(y)); + ret = ecl_make_double_float(ecl_double_float(x) - ecl_double_float(y)); + break; } CASE_FIXNUM_LONG_FLOAT { - return ecl_make_long_float(ecl_fixnum(x) - ecl_long_float(y)); + ret = ecl_make_long_float(ecl_fixnum(x) - ecl_long_float(y)); + break; } CASE_BIGNUM_LONG_FLOAT { - return ecl_make_long_float(ecl_to_long_double(x) - ecl_long_float(y)); + ret = ecl_make_long_float(ecl_to_long_double(x) - ecl_long_float(y)); + break; } CASE_RATIO_LONG_FLOAT { - return ecl_make_long_float(ecl_to_long_double(x) - ecl_long_float(y)); + ret = ecl_make_long_float(ecl_to_long_double(x) - ecl_long_float(y)); + break; } CASE_SINGLE_FLOAT_LONG_FLOAT { - return ecl_make_long_float(ecl_single_float(x) - ecl_long_float(y)); + ret = ecl_make_long_float(ecl_single_float(x) - ecl_long_float(y)); + break; } CASE_DOUBLE_FLOAT_LONG_FLOAT { - return ecl_make_long_float(ecl_double_float(x) - ecl_long_float(y)); + ret = ecl_make_long_float(ecl_double_float(x) - ecl_long_float(y)); + break; } CASE_LONG_FLOAT_FIXNUM { - return ecl_make_long_float(ecl_long_float(x) - ecl_fixnum(y)); + ret = ecl_make_long_float(ecl_long_float(x) - ecl_fixnum(y)); + break; } CASE_LONG_FLOAT_BIGNUM; CASE_LONG_FLOAT_RATIO { - return ecl_make_long_float(ecl_long_float(x) - ecl_to_long_double(y)); + ret = ecl_make_long_float(ecl_long_float(x) - ecl_to_long_double(y)); + break; } CASE_LONG_FLOAT_SINGLE_FLOAT { - return ecl_make_long_float(ecl_long_float(x) - ecl_single_float(y)); + ret = ecl_make_long_float(ecl_long_float(x) - ecl_single_float(y)); + break; } CASE_LONG_FLOAT_DOUBLE_FLOAT { - return ecl_make_long_float(ecl_long_float(x) - ecl_double_float(y)); + ret = ecl_make_long_float(ecl_long_float(x) - ecl_double_float(y)); + break; } CASE_LONG_FLOAT_LONG_FLOAT { - return ecl_make_long_float(ecl_long_float(x) - ecl_long_float(y)); + ret = ecl_make_long_float(ecl_long_float(x) - ecl_long_float(y)); + break; } CASE_LONG_FLOAT_COMPLEX { goto COMPLEX_Y; @@ -179,9 +206,9 @@ ecl_minus(cl_object x, cl_object y) CASE_SINGLE_FLOAT_CSFLOAT; CASE_COMPLEX_CSFLOAT; CASE_CSFLOAT_CSFLOAT { - cl_object aux = ecl_alloc_object(t_csfloat); - ecl_csfloat(aux) = ecl_to_csfloat(x) - ecl_to_csfloat(y); - return aux; + ret = ecl_alloc_object(t_csfloat); + ecl_csfloat(ret) = ecl_to_csfloat(x) - ecl_to_csfloat(y); + break; } /* upgraded type cdfloat */ CASE_CSFLOAT_DOUBLE_FLOAT; @@ -203,9 +230,9 @@ ecl_minus(cl_object x, cl_object y) CASE_COMPLEX_CDFLOAT; CASE_CSFLOAT_CDFLOAT; CASE_CDFLOAT_CDFLOAT { - cl_object aux = ecl_alloc_object(t_cdfloat); - ecl_cdfloat(aux) = ecl_to_cdfloat(x) - ecl_to_cdfloat(y); - return aux; + ret = ecl_alloc_object(t_cdfloat); + ecl_cdfloat(ret) = ecl_to_cdfloat(x) - ecl_to_cdfloat(y); + break; } /* upgraded type clfloat */ CASE_CSFLOAT_LONG_FLOAT; @@ -234,12 +261,15 @@ ecl_minus(cl_object x, cl_object y) CASE_CDFLOAT_CLFLOAT; CASE_CLFLOAT_CLFLOAT { - cl_object aux = ecl_alloc_object(t_clfloat); - ecl_clfloat(aux) = ecl_to_clfloat(x) - ecl_to_clfloat(y); - return aux; + ret = ecl_alloc_object(t_clfloat); + ecl_clfloat(ret) = ecl_to_clfloat(x) - ecl_to_clfloat(y); + break; } #endif CASE_UNKNOWN(@[-],x,y,@[number]); } MATH_DISPATCH2_END; + + ECL_MATHERR_TEST; + return ret; } diff --git a/src/c/numbers/one_minus.d b/src/c/numbers/one_minus.d index 7f655abdc..1a7c12899 100644 --- a/src/c/numbers/one_minus.d +++ b/src/c/numbers/one_minus.d @@ -16,6 +16,11 @@ #include #include +/* INV: FLT_MIN - 1 == FLT_MIN + * DBL_MIN - 1 == DBL_MIN + * LDBL_MIN - 1 == LDBL_MIN + * (no ECL_MATHERR_TEST needed) */ + static cl_object ecl_one_minus_fix(cl_object x) { diff --git a/src/c/numbers/one_plus.d b/src/c/numbers/one_plus.d index 7d56af8ea..5add18b94 100644 --- a/src/c/numbers/one_plus.d +++ b/src/c/numbers/one_plus.d @@ -16,6 +16,11 @@ #include #include +/* INV: FLT_MAX + 1 == FLT_MAX + * DBL_MAX + 1 == DBL_MAX + * LDBL_MAX + 1 == LDBL_MAX + * (no ECL_MATHERR_TEST needed) */ + static cl_object ecl_one_plus_fix(cl_object x) { diff --git a/src/c/numbers/plus.d b/src/c/numbers/plus.d index 02b68dd9c..aceb34b41 100644 --- a/src/c/numbers/plus.d +++ b/src/c/numbers/plus.d @@ -12,10 +12,12 @@ * */ - +#define ECL_INCLUDE_MATH_H #include #include +#pragma STDC FENV_ACCESS ON + @(defun + (&rest nums) cl_object sum = ecl_make_fixnum(0); @ @@ -27,6 +29,9 @@ cl_object ecl_plus(cl_object x, cl_object y) { + cl_object ret; + ECL_MATHERR_CLEAR; + MATH_DISPATCH2_BEGIN(x,y) { CASE_FIXNUM_FIXNUM { @@ -42,10 +47,12 @@ ecl_plus(cl_object x, cl_object y) { return ecl_make_ratio(z, y->ratio.den); } CASE_FIXNUM_SINGLE_FLOAT { - return ecl_make_single_float(ecl_fixnum(x) + ecl_single_float(y)); + ret = ecl_make_single_float(ecl_fixnum(x) + ecl_single_float(y)); + break; } CASE_FIXNUM_DOUBLE_FLOAT { - return ecl_make_double_float(ecl_fixnum(x) + ecl_double_float(y)); + ret = ecl_make_double_float(ecl_fixnum(x) + ecl_double_float(y)); + break; } CASE_BIGNUM_FIXNUM { return _ecl_big_plus_fix(x, ecl_fixnum(y)); @@ -55,11 +62,13 @@ ecl_plus(cl_object x, cl_object y) { } CASE_BIGNUM_SINGLE_FLOAT; CASE_RATIO_SINGLE_FLOAT { - return ecl_make_single_float(ecl_to_float(x) + ecl_single_float(y)); + ret = ecl_make_single_float(ecl_to_float(x) + ecl_single_float(y)); + break; } CASE_BIGNUM_DOUBLE_FLOAT; CASE_RATIO_DOUBLE_FLOAT { - return ecl_make_double_float(ecl_to_double(x) + ecl_double_float(y)); + ret = ecl_make_double_float(ecl_to_double(x) + ecl_double_float(y)); + break; } CASE_RATIO_FIXNUM; CASE_RATIO_BIGNUM { @@ -75,61 +84,79 @@ ecl_plus(cl_object x, cl_object y) { return ecl_make_ratio(z, z1); } CASE_SINGLE_FLOAT_FIXNUM { - return ecl_make_single_float(ecl_single_float(x) + ecl_fixnum(y)); + ret = ecl_make_single_float(ecl_single_float(x) + ecl_fixnum(y)); + break; } CASE_SINGLE_FLOAT_BIGNUM; CASE_SINGLE_FLOAT_RATIO { - return ecl_make_single_float(ecl_single_float(x) + ecl_to_float(y)); + ret = ecl_make_single_float(ecl_single_float(x) + ecl_to_float(y)); + break; } CASE_SINGLE_FLOAT_SINGLE_FLOAT { - return ecl_make_single_float(ecl_single_float(x) + ecl_single_float(y)); + ret = ecl_make_single_float(ecl_single_float(x) + ecl_single_float(y)); + break; } CASE_SINGLE_FLOAT_DOUBLE_FLOAT { - return ecl_make_double_float(ecl_single_float(x) + ecl_double_float(y)); + ret = ecl_make_double_float(ecl_single_float(x) + ecl_double_float(y)); + break; } CASE_DOUBLE_FLOAT_FIXNUM { - return ecl_make_double_float(ecl_double_float(x) + ecl_fixnum(y)); + ret = ecl_make_double_float(ecl_double_float(x) + ecl_fixnum(y)); + break; } CASE_DOUBLE_FLOAT_BIGNUM; CASE_DOUBLE_FLOAT_RATIO { - return ecl_make_double_float(ecl_double_float(x) + ecl_to_double(y)); + ret = ecl_make_double_float(ecl_double_float(x) + ecl_to_double(y)); + break; } CASE_DOUBLE_FLOAT_SINGLE_FLOAT { - return ecl_make_double_float(ecl_double_float(x) + ecl_single_float(y)); + ret = ecl_make_double_float(ecl_double_float(x) + ecl_single_float(y)); + break; } CASE_DOUBLE_FLOAT_DOUBLE_FLOAT { - return ecl_make_double_float(ecl_double_float(x) + ecl_double_float(y)); + ret = ecl_make_double_float(ecl_double_float(x) + ecl_double_float(y)); + break; } CASE_FIXNUM_LONG_FLOAT { - return ecl_make_long_float(ecl_fixnum(x) + ecl_long_float(y)); + ret = ecl_make_long_float(ecl_fixnum(x) + ecl_long_float(y)); + break; } CASE_BIGNUM_LONG_FLOAT { - return ecl_make_long_float(ecl_to_long_double(x) + ecl_long_float(y)); + ret = ecl_make_long_float(ecl_to_long_double(x) + ecl_long_float(y)); + break; } CASE_RATIO_LONG_FLOAT { - return ecl_make_long_float(ecl_to_long_double(x) + ecl_long_float(y)); + ret = ecl_make_long_float(ecl_to_long_double(x) + ecl_long_float(y)); + break; } CASE_SINGLE_FLOAT_LONG_FLOAT { - return ecl_make_long_float(ecl_single_float(x) + ecl_long_float(y)); + ret = ecl_make_long_float(ecl_single_float(x) + ecl_long_float(y)); + break; } CASE_DOUBLE_FLOAT_LONG_FLOAT { - return ecl_make_long_float(ecl_double_float(x) + ecl_long_float(y)); + ret = ecl_make_long_float(ecl_double_float(x) + ecl_long_float(y)); + break; } CASE_LONG_FLOAT_FIXNUM { - return ecl_make_long_float(ecl_long_float(x) + ecl_fixnum(y)); + ret = ecl_make_long_float(ecl_long_float(x) + ecl_fixnum(y)); + break; } CASE_LONG_FLOAT_BIGNUM; CASE_LONG_FLOAT_RATIO { - return ecl_make_long_float(ecl_long_float(x) + ecl_to_long_double(y)); + ret = ecl_make_long_float(ecl_long_float(x) + ecl_to_long_double(y)); + break; } CASE_LONG_FLOAT_SINGLE_FLOAT { - return ecl_make_long_float(ecl_long_float(x) + ecl_single_float(y)); + ret = ecl_make_long_float(ecl_long_float(x) + ecl_single_float(y)); + break; } CASE_LONG_FLOAT_DOUBLE_FLOAT { - return ecl_make_long_float(ecl_long_float(x) + ecl_double_float(y)); + ret = ecl_make_long_float(ecl_long_float(x) + ecl_double_float(y)); + break; } CASE_LONG_FLOAT_LONG_FLOAT { - return ecl_make_long_float(ecl_long_float(x) + ecl_long_float(y)); + ret = ecl_make_long_float(ecl_long_float(x) + ecl_long_float(y)); + break; } CASE_LONG_FLOAT_COMPLEX { goto COMPLEX_Y; @@ -174,9 +201,9 @@ ecl_plus(cl_object x, cl_object y) { CASE_SINGLE_FLOAT_CSFLOAT; CASE_COMPLEX_CSFLOAT; CASE_CSFLOAT_CSFLOAT { - cl_object aux = ecl_alloc_object(t_csfloat); - ecl_csfloat(aux) = ecl_to_csfloat(x) + ecl_to_csfloat(y); - return aux; + ret = ecl_alloc_object(t_csfloat); + ecl_csfloat(ret) = ecl_to_csfloat(x) + ecl_to_csfloat(y); + break; } /* upgraded type cdfloat */ CASE_CSFLOAT_DOUBLE_FLOAT; @@ -198,9 +225,9 @@ ecl_plus(cl_object x, cl_object y) { CASE_COMPLEX_CDFLOAT; CASE_CSFLOAT_CDFLOAT; CASE_CDFLOAT_CDFLOAT { - cl_object aux = ecl_alloc_object(t_cdfloat); - ecl_cdfloat(aux) = ecl_to_cdfloat(x) + ecl_to_cdfloat(y); - return aux; + ret = ecl_alloc_object(t_cdfloat); + ecl_cdfloat(ret) = ecl_to_cdfloat(x) + ecl_to_cdfloat(y); + break; } /* upgraded type clfloat */ CASE_CSFLOAT_LONG_FLOAT; @@ -229,12 +256,15 @@ ecl_plus(cl_object x, cl_object y) { CASE_CDFLOAT_CLFLOAT; CASE_CLFLOAT_CLFLOAT { - cl_object aux = ecl_alloc_object(t_clfloat); - ecl_clfloat(aux) = ecl_to_clfloat(x) + ecl_to_clfloat(y); - return aux; + ret = ecl_alloc_object(t_clfloat); + ecl_clfloat(ret) = ecl_to_clfloat(x) + ecl_to_clfloat(y); + break; } #endif CASE_UNKNOWN(@[+],x,y,@[number]); } MATH_DISPATCH2_END; + + ECL_MATHERR_TEST; + return ret; } diff --git a/src/c/numbers/round.d b/src/c/numbers/round.d index d46eba28d..b1b97060f 100644 --- a/src/c/numbers/round.d +++ b/src/c/numbers/round.d @@ -23,6 +23,8 @@ #endif #include +#pragma STDC FENV_ACCESS ON + @(defun round (x &optional (y OBJNULL)) @ if (narg == 1) @@ -104,6 +106,8 @@ ecl_round1(cl_object x) { const cl_env_ptr the_env = ecl_process_env(); cl_object v0, v1; + ECL_MATHERR_CLEAR; + switch (ecl_t_of(x)) { case t_fixnum: case t_bignum: @@ -138,6 +142,8 @@ ecl_round1(cl_object x) default: FEwrong_type_nth_arg(@[round],1,x,@[real]); } + + ECL_MATHERR_TEST; ecl_return2(the_env, v0, v1); } diff --git a/src/c/numbers/times.d b/src/c/numbers/times.d index db3bfc52c..44b5b0833 100644 --- a/src/c/numbers/times.d +++ b/src/c/numbers/times.d @@ -12,10 +12,12 @@ * */ - +#define ECL_INCLUDE_MATH_H #include #include +#pragma STDC FENV_ACCESS ON + @(defun * (&rest nums) cl_object prod = ecl_make_fixnum(1); @ @@ -28,6 +30,9 @@ cl_object ecl_times(cl_object x, cl_object y) { + cl_object ret; + ECL_MATHERR_CLEAR; + MATH_DISPATCH2_BEGIN(x,y) { CASE_FIXNUM_FIXNUM { @@ -42,10 +47,12 @@ ecl_times(cl_object x, cl_object y) y->ratio.den); } CASE_FIXNUM_SINGLE_FLOAT { - return ecl_make_single_float(ecl_fixnum(x) * ecl_single_float(y)); + ret = ecl_make_single_float(ecl_fixnum(x) * ecl_single_float(y)); + break; } CASE_FIXNUM_DOUBLE_FLOAT { - return ecl_make_double_float(ecl_fixnum(x) * ecl_double_float(y)); + ret = ecl_make_double_float(ecl_fixnum(x) * ecl_double_float(y)); + break; } CASE_BIGNUM_FIXNUM { return _ecl_big_times_fix(x, ecl_fixnum(y)); @@ -54,10 +61,12 @@ ecl_times(cl_object x, cl_object y) return _ecl_big_times_big(x, y); } CASE_BIGNUM_SINGLE_FLOAT { - return ecl_make_single_float(ecl_to_float(x) * ecl_single_float(y)); + ret = ecl_make_single_float(ecl_to_float(x) * ecl_single_float(y)); + break; } CASE_BIGNUM_DOUBLE_FLOAT { - return ecl_make_double_float(ecl_to_double(x) * ecl_double_float(y)); + ret = ecl_make_double_float(ecl_to_double(x) * ecl_double_float(y)); + break; } CASE_RATIO_FIXNUM; CASE_RATIO_BIGNUM { @@ -70,65 +79,84 @@ ecl_times(cl_object x, cl_object y) return ecl_make_ratio(num, den); } CASE_RATIO_SINGLE_FLOAT { - return ecl_make_single_float(ecl_to_float(x) * ecl_single_float(y)); + ret = ecl_make_single_float(ecl_to_float(x) * ecl_single_float(y)); + break; } CASE_RATIO_DOUBLE_FLOAT { - return ecl_make_double_float(ecl_to_double(x) * ecl_double_float(y)); + ret = ecl_make_double_float(ecl_to_double(x) * ecl_double_float(y)); + break; } CASE_SINGLE_FLOAT_FIXNUM { - return ecl_make_single_float(ecl_single_float(x) * ecl_fixnum(y)); + ret = ecl_make_single_float(ecl_single_float(x) * ecl_fixnum(y)); + break; } CASE_SINGLE_FLOAT_BIGNUM; CASE_SINGLE_FLOAT_RATIO { - return ecl_make_single_float(ecl_single_float(x) * ecl_to_float(y)); + ret = ecl_make_single_float(ecl_single_float(x) * ecl_to_float(y)); + break; } CASE_SINGLE_FLOAT_SINGLE_FLOAT { - return ecl_make_single_float(ecl_single_float(x) * ecl_single_float(y)); + ret = ecl_make_single_float(ecl_single_float(x) * ecl_single_float(y)); + break; } CASE_SINGLE_FLOAT_DOUBLE_FLOAT { - return ecl_make_double_float(ecl_single_float(x) * ecl_double_float(y)); + ret = ecl_make_double_float(ecl_single_float(x) * ecl_double_float(y)); + break; } CASE_DOUBLE_FLOAT_FIXNUM { - return ecl_make_double_float(ecl_double_float(x) * ecl_fixnum(y)); + ret = ecl_make_double_float(ecl_double_float(x) * ecl_fixnum(y)); + break; } CASE_DOUBLE_FLOAT_BIGNUM; CASE_DOUBLE_FLOAT_RATIO { - return ecl_make_double_float(ecl_double_float(x) * ecl_to_double(y)); + ret = ecl_make_double_float(ecl_double_float(x) * ecl_to_double(y)); + break; } CASE_DOUBLE_FLOAT_SINGLE_FLOAT { - return ecl_make_double_float(ecl_double_float(x) * ecl_single_float(y)); + ret = ecl_make_double_float(ecl_double_float(x) * ecl_single_float(y)); + break; } CASE_DOUBLE_FLOAT_DOUBLE_FLOAT { - return ecl_make_double_float(ecl_double_float(x) * ecl_double_float(y)); + ret = ecl_make_double_float(ecl_double_float(x) * ecl_double_float(y)); + break; } CASE_FIXNUM_LONG_FLOAT { - return ecl_make_long_float(ecl_fixnum(x) * ecl_long_float(y)); + ret = ecl_make_long_float(ecl_fixnum(x) * ecl_long_float(y)); + break; } CASE_BIGNUM_LONG_FLOAT; CASE_RATIO_LONG_FLOAT { - return ecl_make_long_float(ecl_to_long_double(x) * ecl_long_float(y)); + ret = ecl_make_long_float(ecl_to_long_double(x) * ecl_long_float(y)); + break; } CASE_SINGLE_FLOAT_LONG_FLOAT { - return ecl_make_long_float(ecl_single_float(x) * ecl_long_float(y)); + ret = ecl_make_long_float(ecl_single_float(x) * ecl_long_float(y)); + break; } CASE_DOUBLE_FLOAT_LONG_FLOAT { - return ecl_make_long_float(ecl_double_float(x) * ecl_long_float(y)); + ret = ecl_make_long_float(ecl_double_float(x) * ecl_long_float(y)); + break; } CASE_LONG_FLOAT_FIXNUM { - return ecl_make_long_float(ecl_long_float(x) * ecl_fixnum(y)); + ret = ecl_make_long_float(ecl_long_float(x) * ecl_fixnum(y)); + break; } CASE_LONG_FLOAT_BIGNUM; CASE_LONG_FLOAT_RATIO { - return ecl_make_long_float(ecl_long_float(x) * ecl_to_long_double(y)); + ret = ecl_make_long_float(ecl_long_float(x) * ecl_to_long_double(y)); + break; } CASE_LONG_FLOAT_SINGLE_FLOAT { - return ecl_make_long_float(ecl_long_float(x) * ecl_single_float(y)); + ret = ecl_make_long_float(ecl_long_float(x) * ecl_single_float(y)); + break; } CASE_LONG_FLOAT_DOUBLE_FLOAT { - return ecl_make_long_float(ecl_long_float(x) * ecl_double_float(y)); + ret = ecl_make_long_float(ecl_long_float(x) * ecl_double_float(y)); + break; } CASE_LONG_FLOAT_LONG_FLOAT { - return ecl_make_long_float(ecl_long_float(x) * ecl_long_float(y)); + ret = ecl_make_long_float(ecl_long_float(x) * ecl_long_float(y)); + break; } CASE_LONG_FLOAT_COMPLEX { goto COMPLEX_Y; @@ -179,9 +207,9 @@ ecl_times(cl_object x, cl_object y) CASE_SINGLE_FLOAT_CSFLOAT; CASE_COMPLEX_CSFLOAT; CASE_CSFLOAT_CSFLOAT { - cl_object aux = ecl_alloc_object(t_csfloat); - ecl_csfloat(aux) = ecl_to_csfloat(x) * ecl_to_csfloat(y); - return aux; + ret = ecl_alloc_object(t_csfloat); + ecl_csfloat(ret) = ecl_to_csfloat(x) * ecl_to_csfloat(y); + break; } /* upgraded type cdfloat */ CASE_CSFLOAT_DOUBLE_FLOAT; @@ -203,9 +231,9 @@ ecl_times(cl_object x, cl_object y) CASE_COMPLEX_CDFLOAT; CASE_CSFLOAT_CDFLOAT; CASE_CDFLOAT_CDFLOAT { - cl_object aux = ecl_alloc_object(t_cdfloat); - ecl_cdfloat(aux) = ecl_to_cdfloat(x) * ecl_to_cdfloat(y); - return aux; + ret = ecl_alloc_object(t_cdfloat); + ecl_cdfloat(ret) = ecl_to_cdfloat(x) * ecl_to_cdfloat(y); + break; } /* upgraded type clfloat */ CASE_CSFLOAT_LONG_FLOAT; @@ -234,12 +262,15 @@ ecl_times(cl_object x, cl_object y) CASE_CDFLOAT_CLFLOAT; CASE_CLFLOAT_CLFLOAT { - cl_object aux = ecl_alloc_object(t_clfloat); - ecl_clfloat(aux) = ecl_to_clfloat(x) * ecl_to_clfloat(y); - return aux; + ret = ecl_alloc_object(t_clfloat); + ecl_clfloat(ret) = ecl_to_clfloat(x) * ecl_to_clfloat(y); + break; } #endif CASE_UNKNOWN(@[*],x,y,@[number]); } MATH_DISPATCH2_END; + + ECL_MATHERR_TEST; + return ret; } diff --git a/src/c/numbers/truncate.d b/src/c/numbers/truncate.d index d3a705db9..8d23ecfdf 100644 --- a/src/c/numbers/truncate.d +++ b/src/c/numbers/truncate.d @@ -22,10 +22,14 @@ #endif #include +#pragma STDC FENV_ACCESS ON + cl_object ecl_truncate1(cl_object x) { cl_object v0, v1; + ECL_MATHERR_CLEAR; + switch (ecl_t_of(x)) { case t_fixnum: case t_bignum: @@ -61,6 +65,7 @@ ecl_truncate1(cl_object x) default: FEwrong_type_nth_arg(@[truncate],1,x,@[real]); } + ECL_MATHERR_TEST; { const cl_env_ptr the_env = ecl_process_env(); ecl_return2(the_env, v0, v1); diff --git a/src/h/impl/math_dispatch.h b/src/h/impl/math_dispatch.h index a1ee0fca2..6acf3cffe 100644 --- a/src/h/impl/math_dispatch.h +++ b/src/h/impl/math_dispatch.h @@ -63,7 +63,9 @@ typedef cl_object (*math_one_arg_fn)(cl_object); cl_object ecl_##name(cl_object arg) \ { \ cl_object out; \ + ECL_MATHERR_CLEAR; \ out = ecl_##name##_ne(arg); \ + ECL_MATHERR_TEST; \ return out; \ } diff --git a/src/h/impl/math_fenv.h b/src/h/impl/math_fenv.h index f1c69ea6f..0a93c8e0a 100644 --- a/src/h/impl/math_fenv.h +++ b/src/h/impl/math_fenv.h @@ -95,4 +95,21 @@ # define ECL_WITH_LISP_FPE_END } while (0) #endif +#if defined(HAVE_FENV_H) && defined(ECL_IEEE_FP) && !defined(HAVE_FEENABLEEXCEPT) && !defined(ECL_AVOID_FPE_H) +# define ECL_USED_EXCEPTIONS (FE_DIVBYZERO|FE_INVALID|FE_OVERFLOW|FE_UNDERFLOW) +# define ECL_MATHERR_CLEAR feclearexcept(FE_ALL_EXCEPT) +# define ECL_MATHERR_TEST do { \ + int bits = fetestexcept(ECL_USED_EXCEPTIONS); \ + unlikely_if (bits) { \ + bits &= ecl_process_env()->trap_fpe_bits; \ + if (bits) ecl_deliver_fpe(bits); \ + } \ + } while(0) +#else +# define ECL_MATHERR_CLEAR +# define ECL_MATHERR_TEST +#endif + +extern void ecl_deliver_fpe(int flags); + #endif /* !ECL_MATH_FENV_H */