diff --git a/src/CHANGELOG b/src/CHANGELOG index 4d3237602..fcb584f8d 100644 --- a/src/CHANGELOG +++ b/src/CHANGELOG @@ -28,6 +28,10 @@ ECL 0.9f - Where available, ECL now uses _setjmp/_longjmp for control structures. These functions are faster, as they do not save signals. + - (EXPT 10.0l0 308) failed because the routine EXPT computed too many powers + of 10.0l0, some of which (in particular 10^512) were not required and + overflowed the machine accuracy. + * Foreign function interface (FFI): - ext:c-uint-max and ext:c-ulong-max did not have the right bignum value. diff --git a/src/c/num_sfun.d b/src/c/num_sfun.d index acf643bb5..b54018911 100644 --- a/src/c/num_sfun.d +++ b/src/c/num_sfun.d @@ -138,28 +138,26 @@ cl_expt(cl_object x, cl_object y) if (number_zerop(x)) { if (!number_plusp(ty==t_complex?y->complex.real:y)) FEerror("Cannot raise zero to the power ~S.", 1, y); - return1(number_times(x, y)); - } - if (ty == t_fixnum || ty == t_bignum) { - if (number_minusp(y)) { - z = number_negate(y); - z = cl_expt(x, z); - z = number_divide(MAKE_FIXNUM(1), z); - return1(z); - } + z = number_times(x, y); + } else if (ty != t_fixnum && ty != t_bignum) { + z = cl_log1(x); + z = number_times(z, y); + z = cl_exp(z); + } else if (number_minusp(y)) { + z = number_negate(y); + z = cl_expt(x, z); + z = number_divide(MAKE_FIXNUM(1), z); + } else { z = MAKE_FIXNUM(1); do { /* INV: integer_divide outputs an integer */ if (!number_evenp(y)) z = number_times(z, x); - x = number_times(x, x); y = integer_divide(y, MAKE_FIXNUM(2)); - } while (number_plusp(y)); - return1(z); + if (number_zerop(y)) break; + x = number_times(x, x); + } while (1); } - z = cl_log1(x); - z = number_times(z, y); - z = cl_exp(z); return1(z); }