aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorJoseph Myers <joseph@codesourcery.com>2012-03-05 12:22:46 +0000
committerJoseph Myers <joseph@codesourcery.com>2012-03-05 12:22:46 +0000
commitb7cd39e8f8c5cf2844f20eb03f545d19c4c25987 (patch)
tree1b6c55bdacf30fe15197282cf5d556fe708266ed
parentca811b2256d2e48c7288219e9e11dcbab3000f19 (diff)
downloadglibc-b7cd39e8f8c5cf2844f20eb03f545d19c4c25987.zip
glibc-b7cd39e8f8c5cf2844f20eb03f545d19c4c25987.tar.gz
glibc-b7cd39e8f8c5cf2844f20eb03f545d19c4c25987.tar.bz2
Fix pow in non-default rounding modes (bug 3976).
-rw-r--r--ChangeLog14
-rw-r--r--NEWS10
-rw-r--r--math/libm-test.inc109
-rw-r--r--sysdeps/i386/fpu/libm-test-ulps60
-rw-r--r--sysdeps/ieee754/dbl-64/e_pow.c13
-rw-r--r--sysdeps/x86_64/fpu/libm-test-ulps48
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