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

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.

	[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.

diff --git a/ChangeLog b/ChangeLog
index 62192f3b9350fd3f03701fe9905d6f6fa8104e83..a4376290adfa902b8229a43acb2471e297e248f3 100644 (file)
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,5 +1,17 @@
 2014-05-14  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.
+
        [BZ #16516]
        * sysdeps/ieee754/dbl-64/s_erf.c (efx8): Remove variable.
        (__erf): Scale by 16 instead of 8 in potentially underflowing

diff --git a/NEWS b/NEWS
index 9293f802f8fda0af6cb459a16ed15fd6adff23ff..974d2c80a9325ad9b8ea74b3026990eb81d6549d 100644 (file)
--- a/NEWS
+++ b/NEWS
@@ -17,7 +17,7 @@ Version 2.20
   16713, 16714, 16731, 16739, 16740, 16743, 16754, 16758, 16759, 16760,
   16770, 16786, 16789, 16791, 16799, 16800, 16815, 16823, 16824, 16831,
   16838, 16854, 16876, 16877, 16885, 16888, 16890, 16912, 16916, 16922,
-  16932.
+  16927, 16932.
 
 * The minimum Linux kernel version that this version of the GNU C Library
   can be used with is 2.6.32.

diff --git a/math/libm-test.inc b/math/libm-test.inc
index a4bf0b88c3cc6f18c305c31702cdc560a22a8c9c..b4177e8f8e0288f3e86e7790abaabb97784e928f 100644 (file)
--- 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 98a329179bf8bc74d65d2de3816e6034e387d463..c66781c98592472aa05ced9dfadc97e73a6e2ae8 100644 (file)
--- 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 db9cf337b26ba1e271375bb7b23399e1a43c4cea..fa35d50fd07c4b4e6bc2aa34e93d2fdff19aa193 100644 (file)
--- 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 a832155d17dc9288681fb632e16913c8f27e653f..38d8110550b26a2dc2f0f335cf9f308cf77740d3 100644 (file)
--- 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 aea6c51b2335c49cac4b035c734b894349b3ce95..ccef44ac638786739e8dbca6eb88a974170bfcaa 100644 (file)
--- 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 bbaef680bd520fa415073be987ff3cb84746fcb0..cf9a6db0ef3a3114fbdc035db26fdbfd5c71bdfe 100644 (file)
--- 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 4ba83a46fd4a9a29b4cd011d2bbe1e4d4d1908d9..ad8ae9ce5c6d9100664154ac976ebddfb754ceb0 100644 (file)
--- 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