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

This commit is contained in:
Joseph Myers 2012-03-05 12:22:46 +00:00
parent ca811b2256
commit b7cd39e8f8
6 changed files with 246 additions and 8 deletions

View File

@ -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] [BZ #3976]
* math/libm-test.inc (cosh_test_tonearest): New function. * math/libm-test.inc (cosh_test_tonearest): New function.

10
NEWS
View File

@ -9,11 +9,11 @@ Version 2.16
* The following bugs are resolved with this release: * The following bugs are resolved with this release:
174, 350, 411, 2541, 2547, 2548, 3335, 3992, 4026, 4108, 4596, 4822, 5077, 174, 350, 411, 2541, 2547, 2548, 3335, 3976, 3992, 4026, 4108, 4596, 4822,
5461, 5805, 5993, 6884, 6907, 9739, 9902, 10110, 10135, 10140, 10210, 5077, 5461, 5805, 5993, 6884, 6907, 9739, 9902, 10110, 10135, 10140,
11174, 11322, 11365, 11494, 12047, 13058, 13525, 13526, 13527, 13528, 10210, 11174, 11322, 11365, 11494, 12047, 13058, 13525, 13526, 13527,
13529, 13530, 13531, 13532, 13533, 13547, 13551, 13552, 13553, 13555, 13528, 13529, 13530, 13531, 13532, 13533, 13547, 13551, 13552, 13553,
13559, 13583, 13618, 13637, 13695, 13704, 13706, 13738, 13786 13555, 13559, 13583, 13618, 13637, 13695, 13704, 13706, 13738, 13786
* ISO C11 support: * ISO C11 support:

View File

@ -5318,6 +5318,111 @@ pow_test (void)
END (pow); 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 static void
remainder_test (void) remainder_test (void)
{ {
@ -7218,6 +7323,10 @@ main (int argc, char **argv)
fabs_test (); fabs_test ();
hypot_test (); hypot_test ();
pow_test (); pow_test ();
pow_test_tonearest ();
pow_test_towardzero ();
pow_test_downward ();
pow_test_upward ();
sqrt_test (); sqrt_test ();
/* Error and gamma functions: */ /* Error and gamma functions: */

View File

@ -933,6 +933,42 @@ ifloat: 1
ildouble: 1 ildouble: 1
ldouble: 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 # sin_downward
Test "sin_downward (1) == 0.8414709848078965066525023216302989996226": Test "sin_downward (1) == 0.8414709848078965066525023216302989996226":
ildouble: 1 ildouble: 1
@ -1838,6 +1874,30 @@ ifloat: 1
ildouble: 1 ildouble: 1
ldouble: 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": Function: "sin_downward":
double: 1 double: 1
float: 1 float: 1

View File

@ -1,7 +1,7 @@
/* /*
* IBM Accurate Mathematical Library * IBM Accurate Mathematical Library
* written by International Business Machines Corp. * 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 * 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 * it under the terms of the GNU Lesser General Public License as published by
@ -41,6 +41,7 @@
#include "MathLib.h" #include "MathLib.h"
#include "upow.tbl" #include "upow.tbl"
#include "math_private.h" #include "math_private.h"
#include <fenv.h>
#ifndef SECTION #ifndef SECTION
# define SECTION # define SECTION
@ -84,6 +85,11 @@ __ieee754_pow(double x, double y) {
(u.i[HIGH_HALF]==0 && u.i[LOW_HALF]!=0)) && (u.i[HIGH_HALF]==0 && u.i[LOW_HALF]!=0)) &&
/* 2^-1023< x<= 2^-1023 * 0x1.0000ffffffff */ /* 2^-1023< x<= 2^-1023 * 0x1.0000ffffffff */
(v.i[HIGH_HALF]&0x7fffffff) < 0x4ff00000) { /* if y<-1 or y>1 */ (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)) */ z = log1(x,&aa,&error); /* x^y =e^(y log (X)) */
t = y*134217729.0; t = y*134217729.0;
y1 = t - (t-y); y1 = t - (t-y);
@ -97,7 +103,10 @@ __ieee754_pow(double x, double y) {
a2 = (a-a1)+aa; a2 = (a-a1)+aa;
error = error*ABS(y); error = error*ABS(y);
t = __exp1(a1,a2,1.9e16*error); /* return -10 or 0 if wasn't computed exactly */ 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) { if (x == 0) {

View File

@ -976,6 +976,36 @@ Test "log1p (-0.25) == -0.287682072451780927439219005993827432":
float: 1 float: 1
ifloat: 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 # sin_downward
Test "sin_downward (1) == 0.8414709848078965066525023216302989996226": Test "sin_downward (1) == 0.8414709848078965066525023216302989996226":
ildouble: 1 ildouble: 1
@ -1834,6 +1864,24 @@ Function: "log1p":
float: 1 float: 1
ifloat: 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": Function: "sin_downward":
float: 1 float: 1
ifloat: 1 ifloat: 1