From a31b91ac4235e4fdfa97ffb4fd2eda1c9fc4e158 Mon Sep 17 00:00:00 2001 From: Juan Jose Garcia Ripoll Date: Sun, 15 Nov 2009 23:42:37 +0100 Subject: [PATCH] Simplify LOG/LOG1P so that when passed a NaN they do not enter an infinite recursion. --- src/CHANGELOG | 2 + src/c/num_sfun.d | 137 +++++++++++++++++++++++++++++------------------ 2 files changed, 88 insertions(+), 51 deletions(-) diff --git a/src/CHANGELOG b/src/CHANGELOG index c0e48a473..ccb6c2086 100755 --- a/src/CHANGELOG +++ b/src/CHANGELOG @@ -61,6 +61,8 @@ ECL 9.11.1: - Solved a hard to hit bug in DEFCLASS's routine for detecting collisions in slot names + - LOG and LOG1P did not work properly with NaNs under linux. + * Sockets: - The socket option TCP_NODELAY option has been fixed: it was improperly using diff --git a/src/c/num_sfun.d b/src/c/num_sfun.d index d07c93064..a5e71b3f5 100644 --- a/src/c/num_sfun.d +++ b/src/c/num_sfun.d @@ -253,49 +253,67 @@ cl_object ecl_log1(cl_object x) { cl_object output; - cl_type tx; ECL_MATHERR_CLEAR; AGAIN: - tx = type_of(x); - if (!ECL_NUMBER_TYPE_P(tx)) { - x = ecl_type_error(@'log',"argument",x,@'number'); - goto AGAIN; - } - if (tx == t_complex) { - output = ecl_log1_complex(x->complex.real, x->complex.imag); - } else if (ecl_minusp(x)) { - output = ecl_log1_complex(x, MAKE_FIXNUM(0)); - } else switch (tx) { + switch (type_of(x)) { case t_fixnum: - case t_ratio: + case t_ratio: { + float f = number_to_float(x); + if (f < 0) goto COMPLEX; output = ecl_make_singlefloat(logf(number_to_float(x))); break; + } case t_bignum: { cl_fixnum l = fix(cl_integer_length(x)) - 1; cl_object r = ecl_make_ratio(x, ecl_ash(MAKE_FIXNUM(1), l)); float d = logf(number_to_float(r)) + l * logf(2.0); + if (d < 0) goto COMPLEX; output = ecl_make_singlefloat(d); break; } #ifdef ECL_SHORT_FLOAT - case t_shortfloat: + case t_shortfloat: { + float f = ecl_short_float(d); + if (f < 0) goto COMPLEX; output = make_shortfloat(logf(ecl_short_float(x))); break; + } #endif - case t_singlefloat: - output = ecl_make_singlefloat(logf(sf(x))); + case t_singlefloat: { + float f = sf(x); + if (isnan(f)) goto ISNAN; + if (f < 0) goto COMPLEX; + output = ecl_make_singlefloat(logf(f)); break; - case t_doublefloat: - output = ecl_make_doublefloat(log(df(x))); + } + case t_doublefloat: { + double f = df(x); + if (isnan(f)) goto ISNAN; + if (f < 0) goto COMPLEX; + output = ecl_make_doublefloat(log(f)); break; + } #ifdef ECL_LONG_FLOAT - case t_longfloat: - output = ecl_make_longfloat(logl(ecl_long_float(x))); + case t_longfloat: { + long double f = ecl_long_float(x); + if (isnan(f)) goto ISNAN; + if (f < 0) goto COMPLEX; + output = ecl_make_longfloat(logl(f)); break; + } #endif + case t_complex: + output = ecl_log1_complex(x->complex.real, x->complex.imag); + break; + ISNAN: + output = x; + break; + COMPLEX: + output = ecl_log1_complex(x, MAKE_FIXNUM(0)); + break; default: - /* We do not reach here */ - (void)0; + x = ecl_type_error(@'log',"argument",x,@'number'); + goto AGAIN; } ECL_MATHERR_TEST; return output; @@ -311,45 +329,62 @@ cl_object ecl_log1p(cl_object x) { cl_object output; - cl_type tx; ECL_MATHERR_CLEAR; AGAIN: - tx = type_of(x); - if (!ECL_NUMBER_TYPE_P(tx)) { - x = ecl_type_error(@'log',"argument",x,@'number'); - goto AGAIN; - } - if (tx == t_complex) { - output = ecl_log1(ecl_plus(MAKE_FIXNUM(1), x)); - } else if (ecl_number_compare(x, MAKE_FIXNUM(-1)) < 0) { - output = ecl_log1p(ecl_make_complex(x, MAKE_FIXNUM(0))); - } else { - switch (tx) { + switch (type_of(x)) { case t_fixnum: case t_bignum: - case t_ratio: + case t_ratio: { + float f = number_to_float(x); + if (f < -1) goto COMPLEX; output = ecl_make_singlefloat(log1pf(number_to_float(x))); break; + } #ifdef ECL_SHORT_FLOAT - case t_shortfloat: + case t_shortfloat: { + float f = ecl_short_float(x); + if (isnan(f)) goto ISNAN; + if (f < -1) goto COMPLEX; output = make_shortfloat(log1pf(ecl_short_float(x))); break; -#endif - case t_singlefloat: - output = ecl_make_singlefloat(log1pf(sf(x))); - break; - case t_doublefloat: - output = ecl_make_doublefloat(log1p(df(x))); - break; -#ifdef ECL_LONG_FLOAT - case t_longfloat: - output = ecl_make_longfloat(log1pl(ecl_long_float(x))); - break; -#endif - default: - /* We do not reach here */ - (void)0; } +#endif + case t_singlefloat: { + float f = sf(x); + if (isnan(f)) goto ISNAN; + if (f < -1) goto COMPLEX; + output = ecl_make_singlefloat(log1pf(f)); + break; + } + case t_doublefloat: { + double f = df(x); + if (isnan(f)) goto ISNAN; + if (f < -1) goto COMPLEX; + output = ecl_make_doublefloat(log1p(f)); + break; + } +#ifdef ECL_LONG_FLOAT + case t_longfloat: { + long double f = ecl_long_float(x); + if (isnan(f)) goto ISNAN; + if (f < -1) goto COMPLEX; + output = ecl_make_longfloat(log1pl(f)); + break; + } +#endif + case t_complex: + output = ecl_log1(ecl_plus(MAKE_FIXNUM(1), x)); + break; + ISNAN: + output = x; + break; + COMPLEX: + output = ecl_log1_complex(ecl_plus(x, MAKE_FIXNUM(1)), + MAKE_FIXNUM(0)); + break; + default: + x = ecl_type_error(@'log',"argument",x,@'number'); + goto AGAIN; } ECL_MATHERR_TEST; return output; @@ -509,7 +544,7 @@ cl_object ecl_atan1(cl_object y) { if (type_of(y) == t_complex) { -#if 1 /* ANSI states it should be this first part */ +#if 0 /* ANSI states it should be this first part */ cl_object z = ecl_times(cl_core.imag_unit, y); z = ecl_plus(ecl_log1(ecl_one_plus(z)), ecl_log1(ecl_minus(MAKE_FIXNUM(1), z)));