This is the mail archive of the
libc-alpha@sourceware.org
mailing list for the glibc project.
Fix pow in non-default rounding modes (bug 3976)
- From: "Joseph S. Myers" <joseph at codesourcery dot com>
- To: libc-alpha at sourceware dot org
- Cc: Andreas Jaeger <aj at suse dot com>
- Date: Fri, 2 Mar 2012 22:30:39 +0000 (UTC)
- Subject: Fix pow in non-default rounding modes (bug 3976)
Similar to the patches for other functions, this patch fixes problems
with pow producing wild results (or crashing) in non-default rounding
modes by making it save and restore the rounding mode and do
computations in round-to-nearest.
Unlike the other functions, pow is structured so that only a small
part of the function, with a single return does any significant
computations; the rest does special cases where the computations are
at most a single multiplication or division and so it's fine (indeed
best) for the rest to operate in the rounding mode with which the
function was called. Thus, this patch only saves and restores the
rounding mode within that one block, making the patch much simpler
than the others.
Tested on x86_64 and x86 and ULPs updated based on those tests.
This completes fixing the cases reported in bug 3976 of wild results
from functions in non-default rounding modes.
2012-03-02 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.
diff --git a/math/libm-test.inc b/math/libm-test.inc
index 61c62dd..7c49d62 100644
--- a/math/libm-test.inc
+++ b/math/libm-test.inc
@@ -5210,6 +5210,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)
{
@@ -6993,6 +7098,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 cb22d45..da962cd 100644
--- a/sysdeps/i386/fpu/libm-test-ulps
+++ b/sysdeps/i386/fpu/libm-test-ulps
@@ -897,6 +897,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
@@ -1733,6 +1769,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 0185f28..194c3aa 100644
--- a/sysdeps/x86_64/fpu/libm-test-ulps
+++ b/sysdeps/x86_64/fpu/libm-test-ulps
@@ -932,6 +932,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
@@ -1718,6 +1748,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
--
Joseph S. Myers
joseph@codesourcery.com