Bug#153022: [powerpc libm] exp() in directed rounding modes gives wrong results
tags 153022 + upstream patch moreinfo
quit
On 2013-03-03 16:33:45 -0800, Jonathan Nieder wrote:
> If someone prepares a backport of the fix for 2.13.y, would you
> like to test it?
debdiff attached. Completely untested.
diff -u eglibc-2.13/debian/changelog eglibc-2.13/debian/changelog
--- eglibc-2.13/debian/changelog
+++ eglibc-2.13/debian/changelog
@@ -1,3 +1,11 @@
+eglibc (2.13-38+test) local; urgency=low
+
+ * patches/any/cvs-exp-rounding-mode.diff: __ieee754_exp: save and
+ restore rounding mode and use round-to-nearest for all computations.
+ Closes: #153022
+
+ -- Jonathan Nieder <jrnieder@gmail.com> Sun, 03 Mar 2013 23:04:24 -0800
+
eglibc (2.13-38) unstable; urgency=low
[ Adam Conrad ]
diff -u eglibc-2.13/debian/patches/series eglibc-2.13/debian/patches/series
--- eglibc-2.13/debian/patches/series
+++ eglibc-2.13/debian/patches/series
@@ -376,0 +377 @@
+any/cvs-exp-rounding-mode.diff
only in patch2:
unchanged:
--- eglibc-2.13.orig/debian/patches/any/cvs-exp-rounding-mode.diff
+++ eglibc-2.13/debian/patches/any/cvs-exp-rounding-mode.diff
@@ -0,0 +1,400 @@
+2012-03-02 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.
+---
+ ChangeLog | 14 +++++
+ math/libm-test.inc | 112 ++++++++++++++++++++++++++++++++++++++
+ sysdeps/i386/fpu/libm-test-ulps | 67 +++++++++++++++++++++++
+ sysdeps/ieee754/dbl-64/e_exp.c | 38 ++++++++-----
+ sysdeps/x86_64/fpu/libm-test-ulps | 51 +++++++++++++++++
+ 5 files changed, 268 insertions(+), 14 deletions(-)
+
+diff --git a/math/libm-test.inc b/math/libm-test.inc
+index c6ed7a39..02f51f2e 100644
+--- a/math/libm-test.inc
++++ b/math/libm-test.inc
+@@ -2527,6 +2527,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;
+@@ -6255,6 +6363,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 4b1a9e73..74dbb600 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
+@@ -1138,6 +1183,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 717469e2..19eb134d 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 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
+@@ -39,6 +39,7 @@
+ #include "MathLib.h"
+ #include "uexp.tbl"
+ #include "math_private.h"
++#include <fenv.h>
+
+ double __slowexp(double);
+
+@@ -53,6 +54,11 @@ double __ieee754_exp(double x) {
+ int4 k;
+ #endif
+ int4 i,j,m,n,ex;
++ fenv_t env;
++ double retval;
++
++ feholdexcept (&env);
++ fesetround (FE_TONEAREST);
+
+ junk1.x = x;
+ m = junk1.i[HIGH_HALF];
+@@ -85,18 +91,19 @@ double __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;
+@@ -121,8 +128,8 @@ double __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;
+@@ -135,15 +142,19 @@ double __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:
++ feupdateenv (&env);
++ return retval;
+ }
+
+ /************************************************************************/
+diff --git a/sysdeps/x86_64/fpu/libm-test-ulps b/sysdeps/x86_64/fpu/libm-test-ulps
+index 0ced4be7..4b93a522 100644
+--- a/sysdeps/x86_64/fpu/libm-test-ulps
++++ b/sysdeps/x86_64/fpu/libm-test-ulps
+@@ -502,6 +502,41 @@ ifloat: 2
+ double: 6
+ idouble: 6
+
++# 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
+@@ -1203,6 +1238,22 @@ ifloat: 2
+ double: 6
+ idouble: 6
+
++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
Reply to: