This is the mail archive of the
libc-alpha@sourceware.org
mailing list for the glibc project.
Fix exp in non-default rounding modes (bug 3976)
- From: "Joseph S. Myers" <joseph at codesourcery dot com>
- To: libc-alpha at sourceware dot org
- Date: Wed, 29 Feb 2012 01:10:01 +0000 (UTC)
- Subject: Fix exp in non-default rounding modes (bug 3976)
Bug 3976 discusses how some libm functions produce wildly incorrect
results in rounding modes other than round-to-nearest (they don't need
to be correctly rounded for each mode, but they should still produce
some reasonable approximation to the correct value of the function).
I propose this patch which fixes exp by saving and restoring the
rounding mode around the implementation which depends on
round-to-nearest, and adds associated tests. (Each affected function
would be fixed separately.)
Tested x86_64 and i686 and ULPs updated for those platforms based on
those tests.
2012-02-29 Joseph Myers <joseph@codesourcery.com>
[BZ #3976]
* sysdeps/ieee754/dbl-64/e_exp.c: Include <fenv.h>.
(__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.
diff --git a/math/libm-test.inc b/math/libm-test.inc
index 1016753..f9c2c6d 100644
--- a/math/libm-test.inc
+++ b/math/libm-test.inc
@@ -2526,6 +2526,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)
{
errno = 0;
@@ -6389,6 +6497,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 83a68af..b0bef76 100644
--- a/sysdeps/i386/fpu/libm-test-ulps
+++ b/sysdeps/i386/fpu/libm-test-ulps
@@ -453,6 +453,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
@@ -1182,6 +1227,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
diff --git a/sysdeps/ieee754/dbl-64/e_exp.c b/sysdeps/ieee754/dbl-64/e_exp.c
index cfdb8e2..8231c56 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 <fenv.h>
#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 a2a82e6..5a069b1 100644
--- a/sysdeps/x86_64/fpu/libm-test-ulps
+++ b/sysdeps/x86_64/fpu/libm-test-ulps
@@ -502,6 +502,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
@@ -1260,6 +1295,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
--
Joseph S. Myers
joseph@codesourcery.com