diff --git a/ChangeLog b/ChangeLog index 29c5dc28f0..a10bb8bc91 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,3 +1,17 @@ +2012-03-02 Joseph Myers + + [BZ #3976] + * sysdeps/ieee754/dbl-64/e_exp.c: Include . + (__ieee754_exp): Save and restore rounding mode and use + round-to-nearest for all computations. + * math/libm-test.inc (exp_test_tonearest): New function. + (exp_test_towardzero): Likewise. + (exp_test_downward): Likewise. + (exp_test_upward): Likewise. + (main): Call the new functions. + * sysdeps/i386/fpu/libm-test-ulps: Update. + * sysdeps/x86_64/fpu/libm-test-ulps: Likewise. + 2012-03-01 Chris Demetriou * sysdeps/gnu/errlist-compat.awk: Don't depend on AWK internals to diff --git a/math/libm-test.inc b/math/libm-test.inc index 9f7d4896d8..5bc0d40872 100644 --- a/math/libm-test.inc +++ b/math/libm-test.inc @@ -2531,6 +2531,114 @@ exp_test (void) } +static void +exp_test_tonearest (void) +{ + int save_round_mode; + errno = 0; + FUNC(exp) (0); + if (errno == ENOSYS) + /* Function not implemented. */ + return; + + START (exp_tonearest); + + save_round_mode = fegetround (); + + if (!fesetround (FE_TONEAREST)) + { + TEST_f_f (exp, 1, M_El); + TEST_f_f (exp, 2, M_E2l); + TEST_f_f (exp, 3, M_E3l); + } + + fesetround (save_round_mode); + + END (exp_tonearest); +} + + +static void +exp_test_towardzero (void) +{ + int save_round_mode; + errno = 0; + FUNC(exp) (0); + if (errno == ENOSYS) + /* Function not implemented. */ + return; + + START (exp_towardzero); + + save_round_mode = fegetround (); + + if (!fesetround (FE_TOWARDZERO)) + { + TEST_f_f (exp, 1, M_El); + TEST_f_f (exp, 2, M_E2l); + TEST_f_f (exp, 3, M_E3l); + } + + fesetround (save_round_mode); + + END (exp_towardzero); +} + + +static void +exp_test_downward (void) +{ + int save_round_mode; + errno = 0; + FUNC(exp) (0); + if (errno == ENOSYS) + /* Function not implemented. */ + return; + + START (exp_downward); + + save_round_mode = fegetround (); + + if (!fesetround (FE_DOWNWARD)) + { + TEST_f_f (exp, 1, M_El); + TEST_f_f (exp, 2, M_E2l); + TEST_f_f (exp, 3, M_E3l); + } + + fesetround (save_round_mode); + + END (exp_downward); +} + + +static void +exp_test_upward (void) +{ + int save_round_mode; + errno = 0; + FUNC(exp) (0); + if (errno == ENOSYS) + /* Function not implemented. */ + return; + + START (exp_upward); + + save_round_mode = fegetround (); + + if (!fesetround (FE_UPWARD)) + { + TEST_f_f (exp, 1, M_El); + TEST_f_f (exp, 2, M_E2l); + TEST_f_f (exp, 3, M_E3l); + } + + fesetround (save_round_mode); + + END (exp_upward); +} + + static void exp10_test (void) { @@ -6400,6 +6508,10 @@ main (int argc, char **argv) /* Exponential and logarithmic functions: */ exp_test (); + exp_test_tonearest (); + exp_test_towardzero (); + exp_test_downward (); + exp_test_upward (); exp10_test (); exp2_test (); expm1_test (); diff --git a/sysdeps/i386/fpu/libm-test-ulps b/sysdeps/i386/fpu/libm-test-ulps index 56293973dd..389253e1e3 100644 --- a/sysdeps/i386/fpu/libm-test-ulps +++ b/sysdeps/i386/fpu/libm-test-ulps @@ -460,6 +460,51 @@ Test "exp10 (3) == 1000": ildouble: 8 ldouble: 8 +# exp_downward +Test "exp_downward (1) == e": +ildouble: 1 +ldouble: 1 +Test "exp_downward (2) == e^2": +double: 1 +float: 1 +idouble: 1 +ifloat: 1 +ildouble: 2 +ldouble: 2 +Test "exp_downward (3) == e^3": +double: 1 +float: 1 +idouble: 1 +ifloat: 1 +ildouble: 1 +ldouble: 1 + +# exp_towardzero +Test "exp_towardzero (1) == e": +ildouble: 1 +ldouble: 1 +Test "exp_towardzero (2) == e^2": +double: 1 +float: 1 +idouble: 1 +ifloat: 1 +ildouble: 2 +ldouble: 2 +Test "exp_towardzero (3) == e^3": +double: 1 +float: 1 +idouble: 1 +ifloat: 1 +ildouble: 1 +ldouble: 1 + +# exp_upward +Test "exp_upward (1) == e": +double: 1 +float: 1 +idouble: 1 +ifloat: 1 + # expm1 Test "expm1 (1) == M_El - 1.0": ildouble: 1 @@ -1184,6 +1229,28 @@ Function: "exp10": ildouble: 8 ldouble: 8 +Function: "exp_downward": +double: 1 +float: 1 +idouble: 1 +ifloat: 1 +ildouble: 2 +ldouble: 2 + +Function: "exp_towardzero": +double: 1 +float: 1 +idouble: 1 +ifloat: 1 +ildouble: 2 +ldouble: 2 + +Function: "exp_upward": +double: 1 +float: 1 +idouble: 1 +ifloat: 1 + Function: "expm1": ildouble: 1 ldouble: 1 diff --git a/sysdeps/ieee754/dbl-64/e_exp.c b/sysdeps/ieee754/dbl-64/e_exp.c index cfdb8e2c7d..8231c5698c 100644 --- a/sysdeps/ieee754/dbl-64/e_exp.c +++ b/sysdeps/ieee754/dbl-64/e_exp.c @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * written by International Business Machines Corp. - * Copyright (C) 2001, 2011 Free Software Foundation + * Copyright (C) 2001-2012 Free Software Foundation * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by @@ -38,6 +38,7 @@ #include "MathLib.h" #include "uexp.tbl" #include "math_private.h" +#include #ifndef SECTION # define SECTION @@ -58,6 +59,10 @@ __ieee754_exp(double x) { int4 k; #endif int4 i,j,m,n,ex; + fenv_t env; + double retval; + + libc_feholdexcept_setround (&env, FE_TONEAREST); junk1.x = x; m = junk1.i[HIGH_HALF]; @@ -90,18 +95,19 @@ __ieee754_exp(double x) { rem=(bet + bet*eps)+al*eps; res = al + rem; cor = (al - res) + rem; - if (res == (res+cor*err_0)) return res*binexp.x; - else return __slowexp(x); /*if error is over bound */ + if (res == (res+cor*err_0)) { retval = res*binexp.x; goto ret; } + else { retval = __slowexp(x); goto ret; } /*if error is over bound */ } - if (n <= smallint) return 1.0; + if (n <= smallint) { retval = 1.0; goto ret; } if (n >= badint) { - if (n > infint) return(x+x); /* x is NaN */ - if (n < infint) return ( (x>0) ? (hhuge*hhuge) : (tiny*tiny) ); + if (n > infint) { retval = x+x; goto ret; } /* x is NaN */ + if (n < infint) { retval = (x>0) ? (hhuge*hhuge) : (tiny*tiny); goto ret; } /* x is finite, cause either overflow or underflow */ - if (junk1.i[LOW_HALF] != 0) return (x+x); /* x is NaN */ - return ((x>0)?inf.x:zero ); /* |x| = inf; return either inf or 0 */ + if (junk1.i[LOW_HALF] != 0) { retval = x+x; goto ret; } /* x is NaN */ + retval = (x>0)?inf.x:zero; /* |x| = inf; return either inf or 0 */ + goto ret; } y = x*log2e.x + three51.x; @@ -126,8 +132,8 @@ __ieee754_exp(double x) { if (res < 1.0) {res+=res; cor+=cor; ex-=1;} if (ex >=-1022) { binexp.i[HIGH_HALF] = (1023+ex)<<20; - if (res == (res+cor*err_0)) return res*binexp.x; - else return __slowexp(x); /*if error is over bound */ + if (res == (res+cor*err_0)) { retval = res*binexp.x; goto ret; } + else { retval = __slowexp(x); goto ret; } /*if error is over bound */ } ex = -(1022+ex); binexp.i[HIGH_HALF] = (1023-ex)<<20; @@ -140,15 +146,19 @@ __ieee754_exp(double x) { cor = (t-res)+y; if (res == (res + eps*cor)) { binexp.i[HIGH_HALF] = 0x00100000; - return (res-1.0)*binexp.x; + retval = (res-1.0)*binexp.x; + goto ret; } - else return __slowexp(x); /* if error is over bound */ + else { retval = __slowexp(x); goto ret; } /* if error is over bound */ } else { binexp.i[HIGH_HALF] =(junk1.i[LOW_HALF]+767)<<20; - if (res == (res+cor*err_0)) return res*binexp.x*t256.x; - else return __slowexp(x); + if (res == (res+cor*err_0)) { retval = res*binexp.x*t256.x; goto ret; } + else { retval = __slowexp(x); goto ret; } } + ret: + libc_feupdateenv (&env); + return retval; } #ifndef __ieee754_exp strong_alias (__ieee754_exp, __exp_finite) diff --git a/sysdeps/x86_64/fpu/libm-test-ulps b/sysdeps/x86_64/fpu/libm-test-ulps index dac90babad..74df77ce2a 100644 --- a/sysdeps/x86_64/fpu/libm-test-ulps +++ b/sysdeps/x86_64/fpu/libm-test-ulps @@ -512,6 +512,41 @@ ifloat: 2 ildouble: 8 ldouble: 8 +# exp_downward +Test "exp_downward (1) == e": +ildouble: 1 +ldouble: 1 +Test "exp_downward (2) == e^2": +float: 1 +ifloat: 1 +ildouble: 2 +ldouble: 2 +Test "exp_downward (3) == e^3": +float: 1 +ifloat: 1 +ildouble: 1 +ldouble: 1 + +# exp_towardzero +Test "exp_towardzero (1) == e": +ildouble: 1 +ldouble: 1 +Test "exp_towardzero (2) == e^2": +float: 1 +ifloat: 1 +ildouble: 2 +ldouble: 2 +Test "exp_towardzero (3) == e^3": +float: 1 +ifloat: 1 +ildouble: 1 +ldouble: 1 + +# exp_upward +Test "exp_upward (1) == e": +float: 1 +ifloat: 1 + # expm1 Test "expm1 (0.75) == 1.11700001661267466854536981983709561": double: 1 @@ -1275,6 +1310,22 @@ ifloat: 2 ildouble: 8 ldouble: 8 +Function: "exp_downward": +float: 1 +ifloat: 1 +ildouble: 2 +ldouble: 2 + +Function: "exp_towardzero": +float: 1 +ifloat: 1 +ildouble: 2 +ldouble: 2 + +Function: "exp_upward": +float: 1 +ifloat: 1 + Function: "expm1": double: 1 float: 1