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

6 files changed, 246 insertions, 8 deletions

diff --git a/ChangeLog b/ChangeLog
index 8326f7c..bcb2d50 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,4 +1,16 @@
-2012-03-02  Joseph Myers  <joseph@codesourcery.com>
+2012-03-05  Joseph Myers  <joseph@codesourcery.com>
+
+	[BZ #3976]
+	* sysdeps/ieee754/dbl-64/e_pow.c: Include <fenv.h>.
+	(__ieee754_pow): Save and restore rounding mode and use
+	round-to-nearest for main computations.
+	* math/libm-test.inc (pow_test_tonearest): New function.
+	(pow_test_towardzero): Likewise.
+	(pow_test_downward): Likewise.
+	(pow_test_upward): Likewise.
+	(main): Call the new functions.
+	* sysdeps/i386/fpu/libm-test-ulps: Update.
+	* sysdeps/x86_64/fpu/libm-test-ulps: Likewise.
 
 	[BZ #3976]
 	* math/libm-test.inc (cosh_test_tonearest): New function.
diff --git a/NEWS b/NEWS
index 5d7d093..867e6be 100644
--- a/NEWS
+++ b/NEWS
@@ -9,11 +9,11 @@ Version 2.16
 
 * The following bugs are resolved with this release:
 
-  174, 350, 411, 2541, 2547, 2548, 3335, 3992, 4026, 4108, 4596, 4822, 5077,
-  5461, 5805, 5993, 6884, 6907, 9739, 9902, 10110, 10135, 10140, 10210,
-  11174, 11322, 11365, 11494, 12047, 13058, 13525, 13526, 13527, 13528,
-  13529, 13530, 13531, 13532, 13533, 13547, 13551, 13552, 13553, 13555,
-  13559, 13583, 13618, 13637, 13695, 13704, 13706, 13738, 13786
+  174, 350, 411, 2541, 2547, 2548, 3335, 3976, 3992, 4026, 4108, 4596, 4822,
+  5077, 5461, 5805, 5993, 6884, 6907, 9739, 9902, 10110, 10135, 10140,
+  10210, 11174, 11322, 11365, 11494, 12047, 13058, 13525, 13526, 13527,
+  13528, 13529, 13530, 13531, 13532, 13533, 13547, 13551, 13552, 13553,
+  13555, 13559, 13583, 13618, 13637, 13695, 13704, 13706, 13738, 13786
 
 * ISO C11 support:
 
diff --git a/math/libm-test.inc b/math/libm-test.inc
index 684955e..9bdbc4c 100644
--- a/math/libm-test.inc
+++ b/math/libm-test.inc
@@ -5318,6 +5318,111 @@ pow_test (void)
   END (pow);
 }
 
+
+static void
+pow_test_tonearest (void)
+{
+  int save_round_mode;
+  errno = 0;
+  FUNC(pow) (0, 0);
+  if (errno == ENOSYS)
+    /* Function not implemented.  */
+    return;
+
+  START (pow_tonearest);
+
+  save_round_mode = fegetround ();
+
+  if (!fesetround (FE_TONEAREST))
+    {
+      TEST_ff_f (pow, 1.0625L, 1.125L, 1.070582293028761362162622578677070098674L);
+      TEST_ff_f (pow, 1.5L, 1.03125L, 1.519127098714743184071644334163037684948L);
+    }
+
+  fesetround (save_round_mode);
+
+  END (pow_tonearest);
+}
+
+
+static void
+pow_test_towardzero (void)
+{
+  int save_round_mode;
+  errno = 0;
+  FUNC(pow) (0, 0);
+  if (errno == ENOSYS)
+    /* Function not implemented.  */
+    return;
+
+  START (pow_towardzero);
+
+  save_round_mode = fegetround ();
+
+  if (!fesetround (FE_TOWARDZERO))
+    {
+      TEST_ff_f (pow, 1.0625L, 1.125L, 1.070582293028761362162622578677070098674L);
+      TEST_ff_f (pow, 1.5L, 1.03125L, 1.519127098714743184071644334163037684948L);
+    }
+
+  fesetround (save_round_mode);
+
+  END (pow_towardzero);
+}
+
+
+static void
+pow_test_downward (void)
+{
+  int save_round_mode;
+  errno = 0;
+  FUNC(pow) (0, 0);
+  if (errno == ENOSYS)
+    /* Function not implemented.  */
+    return;
+
+  START (pow_downward);
+
+  save_round_mode = fegetround ();
+
+  if (!fesetround (FE_DOWNWARD))
+    {
+      TEST_ff_f (pow, 1.0625L, 1.125L, 1.070582293028761362162622578677070098674L);
+      TEST_ff_f (pow, 1.5L, 1.03125L, 1.519127098714743184071644334163037684948L);
+    }
+
+  fesetround (save_round_mode);
+
+  END (pow_downward);
+}
+
+
+static void
+pow_test_upward (void)
+{
+  int save_round_mode;
+  errno = 0;
+  FUNC(pow) (0, 0);
+  if (errno == ENOSYS)
+    /* Function not implemented.  */
+    return;
+
+  START (pow_upward);
+
+  save_round_mode = fegetround ();
+
+  if (!fesetround (FE_UPWARD))
+    {
+      TEST_ff_f (pow, 1.0625L, 1.125L, 1.070582293028761362162622578677070098674L);
+      TEST_ff_f (pow, 1.5L, 1.03125L, 1.519127098714743184071644334163037684948L);
+    }
+
+  fesetround (save_round_mode);
+
+  END (pow_upward);
+}
+
+
 static void
 remainder_test (void)
 {
@@ -7218,6 +7323,10 @@ main (int argc, char **argv)
   fabs_test ();
   hypot_test ();
   pow_test ();
+  pow_test_tonearest ();
+  pow_test_towardzero ();
+  pow_test_downward ();
+  pow_test_upward ();
   sqrt_test ();
 
   /* Error and gamma functions:  */
diff --git a/sysdeps/i386/fpu/libm-test-ulps b/sysdeps/i386/fpu/libm-test-ulps
index e17bc53..049d6d1 100644
--- a/sysdeps/i386/fpu/libm-test-ulps
+++ b/sysdeps/i386/fpu/libm-test-ulps
@@ -933,6 +933,42 @@ ifloat: 1
 ildouble: 1
 ldouble: 1
 
+# pow_downward
+Test "pow_downward (1.0625, 1.125) == 1.070582293028761362162622578677070098674":
+double: 1
+idouble: 1
+ildouble: 1
+ldouble: 1
+Test "pow_downward (1.5, 1.03125) == 1.519127098714743184071644334163037684948":
+float: 1
+ifloat: 1
+ildouble: 1
+ldouble: 1
+
+# pow_towardzero
+Test "pow_towardzero (1.0625, 1.125) == 1.070582293028761362162622578677070098674":
+double: 1
+idouble: 1
+ildouble: 1
+ldouble: 1
+Test "pow_towardzero (1.5, 1.03125) == 1.519127098714743184071644334163037684948":
+float: 1
+ifloat: 1
+ildouble: 1
+ldouble: 1
+
+# pow_upward
+Test "pow_upward (1.0625, 1.125) == 1.070582293028761362162622578677070098674":
+float: 1
+ifloat: 1
+ildouble: 1
+ldouble: 1
+Test "pow_upward (1.5, 1.03125) == 1.519127098714743184071644334163037684948":
+double: 1
+idouble: 1
+ildouble: 1
+ldouble: 1
+
 # sin_downward
 Test "sin_downward (1) == 0.8414709848078965066525023216302989996226":
 ildouble: 1
@@ -1838,6 +1874,30 @@ ifloat: 1
 ildouble: 1
 ldouble: 1
 
+Function: "pow_downward":
+double: 1
+float: 1
+idouble: 1
+ifloat: 1
+ildouble: 1
+ldouble: 1
+
+Function: "pow_towardzero":
+double: 1
+float: 1
+idouble: 1
+ifloat: 1
+ildouble: 1
+ldouble: 1
+
+Function: "pow_upward":
+double: 1
+float: 1
+idouble: 1
+ifloat: 1
+ildouble: 1
+ldouble: 1
+
 Function: "sin_downward":
 double: 1
 float: 1
diff --git a/sysdeps/ieee754/dbl-64/e_pow.c b/sysdeps/ieee754/dbl-64/e_pow.c
index 28435fd..f668b4b 100644
--- a/sysdeps/ieee754/dbl-64/e_pow.c
+++ b/sysdeps/ieee754/dbl-64/e_pow.c
@@ -1,7 +1,7 @@
 /*
  * IBM Accurate Mathematical Library
  * written by International Business Machines Corp.
- * Copyright (C) 2001, 2002, 2004, 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
@@ -41,6 +41,7 @@
 #include "MathLib.h"
 #include "upow.tbl"
 #include "math_private.h"
+#include <fenv.h>
 
 #ifndef SECTION
 # define SECTION
@@ -84,6 +85,11 @@ __ieee754_pow(double x, double y) {
        (u.i[HIGH_HALF]==0 && u.i[LOW_HALF]!=0))  &&
 				      /*   2^-1023< x<= 2^-1023 * 0x1.0000ffffffff */
       (v.i[HIGH_HALF]&0x7fffffff) < 0x4ff00000) {              /* if y<-1 or y>1   */
+    fenv_t env;
+    double retval;
+
+    libc_feholdexcept_setround (&env, FE_TONEAREST);
+
     z = log1(x,&aa,&error);                                 /* x^y  =e^(y log (X)) */
     t = y*134217729.0;
     y1 = t - (t-y);
@@ -97,7 +103,10 @@ __ieee754_pow(double x, double y) {
     a2 = (a-a1)+aa;
     error = error*ABS(y);
     t = __exp1(a1,a2,1.9e16*error);     /* return -10 or 0 if wasn't computed exactly */
-    return (t>0)?t:power1(x,y);
+    retval = (t>0)?t:power1(x,y);
+
+    libc_feupdateenv (&env);
+    return retval;
   }
 
   if (x == 0) {
diff --git a/sysdeps/x86_64/fpu/libm-test-ulps b/sysdeps/x86_64/fpu/libm-test-ulps
index dd9e130..269dca6 100644
--- a/sysdeps/x86_64/fpu/libm-test-ulps
+++ b/sysdeps/x86_64/fpu/libm-test-ulps
@@ -976,6 +976,36 @@ Test "log1p (-0.25) == -0.287682072451780927439219005993827432":
 float: 1
 ifloat: 1
 
+# pow_downward
+Test "pow_downward (1.0625, 1.125) == 1.070582293028761362162622578677070098674":
+ildouble: 1
+ldouble: 1
+Test "pow_downward (1.5, 1.03125) == 1.519127098714743184071644334163037684948":
+float: 1
+ifloat: 1
+ildouble: 1
+ldouble: 1
+
+# pow_towardzero
+Test "pow_towardzero (1.0625, 1.125) == 1.070582293028761362162622578677070098674":
+ildouble: 1
+ldouble: 1
+Test "pow_towardzero (1.5, 1.03125) == 1.519127098714743184071644334163037684948":
+float: 1
+ifloat: 1
+ildouble: 1
+ldouble: 1
+
+# pow_upward
+Test "pow_upward (1.0625, 1.125) == 1.070582293028761362162622578677070098674":
+float: 1
+ifloat: 1
+ildouble: 1
+ldouble: 1
+Test "pow_upward (1.5, 1.03125) == 1.519127098714743184071644334163037684948":
+ildouble: 1
+ldouble: 1
+
 # sin_downward
 Test "sin_downward (1) == 0.8414709848078965066525023216302989996226":
 ildouble: 1
@@ -1834,6 +1864,24 @@ Function: "log1p":
 float: 1
 ifloat: 1
 
+Function: "pow_downward":
+float: 1
+ifloat: 1
+ildouble: 1
+ldouble: 1
+
+Function: "pow_towardzero":
+float: 1
+ifloat: 1
+ildouble: 1
+ldouble: 1
+
+Function: "pow_upward":
+float: 1
+ifloat: 1
+ildouble: 1
+ldouble: 1
+
 Function: "sin_downward":
 float: 1
 ifloat: 1