Fix exp in non-default rounding modes (bug 3976).

ChangeLog

@@ -1,3 +1,17 @@
+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.
+
 2012-03-01  Chris Demetriou  <cgd@google.com>
 
 	* sysdeps/gnu/errlist-compat.awk: Don't depend on AWK internals to

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 ();

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

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;
+    if  (res == (res+cor*err_0)) { retval = res*binexp.x; goto ret; }
-    else return __slowexp(x); /*if error is over bound */
+    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) { retval = x+x; goto ret; }               /* x is NaN */
-    if (n < infint) return ( (x>0) ? (hhuge*hhuge) : (tiny*tiny) );
+    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  */
+    if (junk1.i[LOW_HALF] != 0) { retval = x+x; goto ret; } /*  x is NaN  */
-    return ((x>0)?inf.x:zero );             /* |x| = inf;  return either inf or 0 */
+    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;
+      if  (res == (res+cor*err_0)) { retval = res*binexp.x; goto ret; }
-      else return __slowexp(x); /*if error is over bound */
+      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;
+    if (res == (res+cor*err_0)) { retval = res*binexp.x*t256.x; goto ret; }
-    else return __slowexp(x);
+    else { retval = __slowexp(x); goto ret; }
   }
+ ret:
+  libc_feupdateenv (&env);
+  return retval;
 }
 #ifndef __ieee754_exp
 strong_alias (__ieee754_exp, __exp_finite)

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