From 28afd92dbdb4fef4358051aad5cb944a9527a4b5 Mon Sep 17 00:00:00 2001 From: Joseph Myers Date: Fri, 2 Mar 2012 15:12:53 +0000 Subject: Fix exp in non-default rounding modes (bug 3976). --- sysdeps/i386/fpu/libm-test-ulps | 67 +++++++++++++++++++++++++++++++++++++++ sysdeps/ieee754/dbl-64/e_exp.c | 38 ++++++++++++++-------- sysdeps/x86_64/fpu/libm-test-ulps | 51 +++++++++++++++++++++++++++++ 3 files changed, 142 insertions(+), 14 deletions(-) (limited to 'sysdeps') 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 -- cgit 1.4.1