Fix log (1) in round-downward mode (bug 16731).

According to ISO C Annex F, log (1) should be +0 in all rounding
modes, but some implementations in glibc wrongly return -0 in
round-downward mode (mapping to log1p (x - 1) is problematic because 1
- 1 is -0 in round-downward mode, and log1p (-0) is -0).  This patch
fixes this.  (It helps with some implementations of other functions
such as acosh, log2 and log10 that call out to log, but not enough to
enable all-rounding-modes testing for those functions without further
fixes to other implementations of them.)

Tested x86_64 and x86 and ulps updated accordingly, and did spot tests
for mips64 for the ldbl-128 fix, and i586 for the sysdeps/i386/fpu
implementations shadowed by those in sysdeps/i386/i686/fpu.

	[BZ #16731]
	* sysdeps/i386/fpu/e_log.S (__ieee754_log): Take absolute value
	when x - 1 is zero.
	* sysdeps/i386/fpu/e_logf.S (__ieee754_logf): Likewise.
	* sysdeps/i386/fpu/e_logl.S (__ieee754_logl): Likewise.
	* sysdeps/i386/i686/fpu/e_logl.S (__ieee754_logl): Likewise.
	* sysdeps/ieee754/dbl-64/e_log.c (__ieee754_log): Return +0 when
	argument is 1.
	* sysdeps/ieee754/ldbl-128/e_logl.c (__ieee754_logl): Likewise.
	* sysdeps/x86_64/fpu/e_logl.S: Take absolute value when x - 1 is
	zero.
	* math/libm-test.inc (log_test): Use ALL_RM_TEST.
	* sysdeps/i386/fpu/libm-test-ulps: Update.
	* sysdeps/x86_64/fpu/libm-test-ulps: Likewise.
This commit is contained in:
Joseph Myers 2014-03-21 18:13:58 +00:00
parent a387428ca7
commit f7be737659
12 changed files with 88 additions and 9 deletions

View File

@ -1,3 +1,20 @@
2014-03-21 Joseph Myers <joseph@codesourcery.com>
[BZ #16731]
* sysdeps/i386/fpu/e_log.S (__ieee754_log): Take absolute value
when x - 1 is zero.
* sysdeps/i386/fpu/e_logf.S (__ieee754_logf): Likewise.
* sysdeps/i386/fpu/e_logl.S (__ieee754_logl): Likewise.
* sysdeps/i386/i686/fpu/e_logl.S (__ieee754_logl): Likewise.
* sysdeps/ieee754/dbl-64/e_log.c (__ieee754_log): Return +0 when
argument is 1.
* sysdeps/ieee754/ldbl-128/e_logl.c (__ieee754_logl): Likewise.
* sysdeps/x86_64/fpu/e_logl.S: Take absolute value when x - 1 is
zero.
* math/libm-test.inc (log_test): Use ALL_RM_TEST.
* sysdeps/i386/fpu/libm-test-ulps: Update.
* sysdeps/x86_64/fpu/libm-test-ulps: Likewise.
2014-03-21 Siddhesh Poyarekar <siddhesh@redhat.com>
* scripts/bench.pl: Remove file.

2
NEWS
View File

@ -11,7 +11,7 @@ Version 2.20
15347, 15804, 15894, 16447, 16532, 16545, 16574, 16600, 16609, 16610,
16611, 16613, 16623, 16632, 16639, 16642, 16649, 16670, 16674, 16677,
16680, 16683, 16689, 16695, 16701, 16706, 16707.
16680, 16683, 16689, 16695, 16701, 16706, 16707, 16731.
* Running the testsuite no longer terminates as soon as a test fails.
Instead, a file tests.sum (xtests.sum from "make xcheck") is generated,

View File

@ -7786,9 +7786,7 @@ static const struct test_f_f_data log_test_data[] =
static void
log_test (void)
{
START (log, 0);
RUN_TEST_LOOP_f_f (log, log_test_data, );
END;
ALL_RM_TEST (log, 0, log_test_data, RUN_TEST_LOOP_f_f, END);
}

View File

@ -46,7 +46,13 @@ ENTRY(__ieee754_log)
fnstsw // x-1 : x : log(2)
andb $0x45, %ah
jz 2f
fstp %st(1) // x-1 : log(2)
fxam
fnstsw
andb $0x45, %ah
cmpb $0x40, %ah
jne 5f
fabs // log(1) is +0 in all rounding modes.
5: fstp %st(1) // x-1 : log(2)
fyl2xp1 // log(x)
ret

View File

@ -47,7 +47,13 @@ ENTRY(__ieee754_logf)
fnstsw // x-1 : x : log(2)
andb $0x45, %ah
jz 2f
fstp %st(1) // x-1 : log(2)
fxam
fnstsw
andb $0x45, %ah
cmpb $0x40, %ah
jne 5f
fabs // log(1) is +0 in all rounding modes.
5: fstp %st(1) // x-1 : log(2)
fyl2xp1 // log(x)
ret

View File

@ -47,7 +47,13 @@ ENTRY(__ieee754_logl)
fnstsw // x-1 : x : log(2)
andb $0x45, %ah
jz 2f
fstp %st(1) // x-1 : log(2)
fxam
fnstsw
andb $0x45, %ah
cmpb $0x40, %ah
jne 5f
fabs // log(1) is +0 in all rounding modes.
5: fstp %st(1) // x-1 : log(2)
fyl2xp1 // log(x)
ret

View File

@ -1096,6 +1096,18 @@ Function: "log1p":
ildouble: 1
ldouble: 1
Function: "log_downward":
ildouble: 1
ldouble: 1
Function: "log_towardzero":
ildouble: 1
ldouble: 1
Function: "log_upward":
ildouble: 1
ldouble: 1
Function: "pow":
ildouble: 1
ldouble: 1

View File

@ -46,7 +46,13 @@ ENTRY(__ieee754_logl)
fcomip %st(1) // |x-1| : x-1 : x : log(2)
fstp %st(0) // x-1 : x : log(2)
jc 2f
fstp %st(1) // x-1 : log(2)
fxam
fnstsw
andb $0x45, %ah
cmpb $0x40, %ah
jne 4f
fabs // log(1) is +0 in all rounding modes.
4: fstp %st(1) // x-1 : log(2)
fyl2xp1 // log(x)
ret

View File

@ -96,6 +96,10 @@ __ieee754_log (double x)
if (__glibc_likely (ABS (w) > U03))
goto case_03;
/* log (1) is +0 in all rounding modes. */
if (w == 0.0)
return 0.0;
/*--- Stage I, the case abs(x-1) < 0.03 */
t8 = MHALF * w;

View File

@ -240,6 +240,8 @@ __ieee754_logl(long double x)
/* On this interval the table is not used due to cancellation error. */
if ((x <= 1.0078125L) && (x >= 0.9921875L))
{
if (x == 1.0L)
return 0.0L;
z = x - 1.0L;
k = 64;
t.value = 1.0L;

View File

@ -45,7 +45,13 @@ ENTRY(__ieee754_logl)
fnstsw // x-1 : x : log(2)
andb $0x45, %ah
jz 2f
fstp %st(1) // x-1 : log(2)
fxam
fnstsw
andb $0x45, %ah
cmpb $0x40, %ah
jne 5f
fabs // log(1) is +0 in all rounding modes.
5: fstp %st(1) // x-1 : log(2)
fyl2xp1 // log(x)
ret

View File

@ -1170,6 +1170,22 @@ ifloat: 1
ildouble: 1
ldouble: 1
Function: "log_downward":
float: 1
ifloat: 1
ildouble: 1
ldouble: 1
Function: "log_towardzero":
ildouble: 1
ldouble: 1
Function: "log_upward":
float: 1
ifloat: 1
ildouble: 1
ldouble: 1
Function: "pow":
float: 1
ifloat: 1