diff mbox

Fix acosh (1) in round-downward mode (bug 16927)

Message ID Pine.LNX.4.64.1405091637460.14073@digraph.polyomino.org.uk
State New
Headers show

Commit Message

Joseph Myers May 9, 2014, 4:38 p.m. UTC 
According to C99 and C11 Annex F, acosh (1) should be +0 in all
rounding modes.  However, some implementations in glibc wrongly return
-0 in round-downward mode (which is what you get if you end up
computing log1p (-0), via 1 - 1 being -0 in round-downward mode).
This patch fixes the problem implementations, by correcting the test
for an exact 1 value in the ldbl-96 implementation to allow for the
explicit high bit of the mantissa, and by inserting fabs instructions
in the i386 implementations; tests of acosh are duly converted to
ALL_RM_TEST.  I believe all the other sysdeps/ieee754 implementations
are already OK (I haven't checked the ia64 versions, but if buggy then
that will be obvious from the results of test runs after this patch is
in).

Tested x86_64 and x86 and ulps updated accordingly.

2014-05-09  Joseph Myers  <joseph@codesourcery.com>

	[BZ #16927]
	* sysdeps/i386/fpu/e_acosh.S (__ieee754_acosh): Use fabs on x-1
	value.
	* sysdeps/i386/fpu/e_acoshf.S (__ieee754_acoshf): Likewise.
	* sysdeps/i386/fpu/e_acoshl.S (__ieee754_acoshl): Likewise.
	* sysdeps/ieee754/ldbl-96/e_acoshl.c (__ieee754_acoshl): Correct
	for explicit high bit of mantissa when testing for argument equal
	to 1.
	* math/libm-test.inc (acosh_test): Use ALL_RM_TEST.
	* sysdeps/i386/fpu/libm-test-ulps: Update.
	* sysdeps/x86_64/fpu/libm-test-ulps: Likewise.

Comments

Andreas Jaeger May 14, 2014, 10:26 a.m. UTC | #1 
On 05/09/2014 12:38 PM, Joseph S. Myers wrote:
> According to C99 and C11 Annex F, acosh (1) should be +0 in all
> rounding modes.  However, some implementations in glibc wrongly return
> -0 in round-downward mode (which is what you get if you end up
> computing log1p (-0), via 1 - 1 being -0 in round-downward mode).
> This patch fixes the problem implementations, by correcting the test
> for an exact 1 value in the ldbl-96 implementation to allow for the
> explicit high bit of the mantissa, and by inserting fabs instructions
> in the i386 implementations; tests of acosh are duly converted to
> ALL_RM_TEST.  I believe all the other sysdeps/ieee754 implementations
> are already OK (I haven't checked the ia64 versions, but if buggy then
> that will be obvious from the results of test runs after this patch is
> in).

ok, thanks,
Andreas
diff mbox

Patch

diff --git a/math/libm-test.inc b/math/libm-test.inc
index a4bf0b8..b4177e8 100644
--- a/math/libm-test.inc
+++ b/math/libm-test.inc
@@ -1792,9 +1792,7 @@  static const struct test_f_f_data acosh_test_data[] =
 static void
 acosh_test (void)
 {
-  START (acosh, 0);
-  RUN_TEST_LOOP_f_f (acosh, acosh_test_data, );
-  END;
+  ALL_RM_TEST (acosh, 0, acosh_test_data, RUN_TEST_LOOP_f_f, END);
 }
 
 static const struct test_f_f_data asin_test_data[] =
diff --git a/sysdeps/i386/fpu/e_acosh.S b/sysdeps/i386/fpu/e_acosh.S
index 98a3291..c66781c 100644
--- a/sysdeps/i386/fpu/e_acosh.S
+++ b/sysdeps/i386/fpu/e_acosh.S
@@ -52,6 +52,7 @@  ENTRY(__ieee754_acosh)
 
 	// 1 <= x <= 2 => y = log1p(x-1+sqrt(2*(x-1)+(x-1)^2))
 	fsubl	MO(one)			// x-1 : log(2)
+	fabs				// acosh(1) is +0 in all rounding modes
 	fld	%st			// x-1 : x-1 : log(2)
 	fmul	%st(1)			// (x-1)^2 : x-1 : log(2)
 	fadd	%st(1)			// x-1+(x-1)^2 : x-1 : log(2)
diff --git a/sysdeps/i386/fpu/e_acoshf.S b/sysdeps/i386/fpu/e_acoshf.S
index db9cf33..fa35d50 100644
--- a/sysdeps/i386/fpu/e_acoshf.S
+++ b/sysdeps/i386/fpu/e_acoshf.S
@@ -52,6 +52,7 @@  ENTRY(__ieee754_acoshf)
 
 	// 1 <= x <= 2 => y = log1p(x-1+sqrt(2*(x-1)+(x-1)^2))
 	fsubl	MO(one)			// x-1 : log(2)
+	fabs				// acosh(1) is +0 in all rounding modes
 	fld	%st			// x-1 : x-1 : log(2)
 	fmul	%st(1)			// (x-1)^2 : x-1 : log(2)
 	fadd	%st(1)			// x-1+(x-1)^2 : x-1 : log(2)
diff --git a/sysdeps/i386/fpu/e_acoshl.S b/sysdeps/i386/fpu/e_acoshl.S
index a832155..38d8110 100644
--- a/sysdeps/i386/fpu/e_acoshl.S
+++ b/sysdeps/i386/fpu/e_acoshl.S
@@ -59,6 +59,7 @@  ENTRY(__ieee754_acoshl)
 
 	// 1 <= x <= 2 => y = log1p(x-1+sqrt(2*(x-1)+(x-1)^2))
 	fsubl	MO(one)			// x-1 : log(2)
+	fabs				// acosh(1) is +0 in all rounding modes
 	fld	%st			// x-1 : x-1 : log(2)
 	fmul	%st(1)			// (x-1)^2 : x-1 : log(2)
 	fadd	%st(1)			// x-1+(x-1)^2 : x-1 : log(2)
diff --git a/sysdeps/i386/fpu/libm-test-ulps b/sysdeps/i386/fpu/libm-test-ulps
index aea6c51..ccef44a 100644
--- a/sysdeps/i386/fpu/libm-test-ulps
+++ b/sysdeps/i386/fpu/libm-test-ulps
@@ -21,6 +21,18 @@  Function: "acos_upward":
 ildouble: 1
 ldouble: 1
 
+Function: "acosh_downward":
+ildouble: 2
+ldouble: 2
+
+Function: "acosh_towardzero":
+ildouble: 2
+ldouble: 2
+
+Function: "acosh_upward":
+ildouble: 1
+ldouble: 1
+
 Function: "asin_downward":
 double: 1
 float: 1
diff --git a/sysdeps/ieee754/ldbl-96/e_acoshl.c b/sysdeps/ieee754/ldbl-96/e_acoshl.c
index bbaef68..cf9a6db 100644
--- a/sysdeps/ieee754/ldbl-96/e_acoshl.c
+++ b/sysdeps/ieee754/ldbl-96/e_acoshl.c
@@ -48,7 +48,7 @@  __ieee754_acoshl(long double x)
 		return x+x;
 	    } else
 		return __ieee754_logl(x)+ln2;	/* acoshl(huge)=logl(2x) */
-	} else if(((se-0x3fff)|i0|i1)==0) {
+	} else if(((se-0x3fff)|(i0^0x80000000)|i1)==0) {
 	    return 0.0;			/* acosh(1) = 0 */
 	} else if (se > 0x4000) {	/* 2**28 > x > 2 */
 	    t=x*x;
diff --git a/sysdeps/x86_64/fpu/libm-test-ulps b/sysdeps/x86_64/fpu/libm-test-ulps
index 4ba83a4..ad8ae9c 100644
--- a/sysdeps/x86_64/fpu/libm-test-ulps
+++ b/sysdeps/x86_64/fpu/libm-test-ulps
@@ -31,6 +31,21 @@  Function: "acosh":
 double: 1
 idouble: 1
 
+Function: "acosh_downward":
+float: 1
+ildouble: 1
+ldouble: 2
+
+Function: "acosh_towardzero":
+float: 1
+ildouble: 1
+ldouble: 2
+
+Function: "acosh_upward":
+double: 1
+ildouble: 1
+ldouble: 1
+
 Function: "asin_downward":
 double: 1
 float: 1