From 913d03c864ea2547e97f8d2d30fc71a008d4d103 Mon Sep 17 00:00:00 2001 From: Joseph Myers Date: Wed, 14 May 2014 12:35:40 +0000 Subject: [PATCH] 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. --- ChangeLog | 12 ++++++++++++ NEWS | 2 +- math/libm-test.inc | 4 +--- sysdeps/i386/fpu/e_acosh.S | 1 + sysdeps/i386/fpu/e_acoshf.S | 1 + sysdeps/i386/fpu/e_acoshl.S | 1 + sysdeps/i386/fpu/libm-test-ulps | 12 ++++++++++++ sysdeps/ieee754/ldbl-96/e_acoshl.c | 2 +- sysdeps/x86_64/fpu/libm-test-ulps | 15 +++++++++++++++ 9 files changed, 45 insertions(+), 5 deletions(-) diff --git a/ChangeLog b/ChangeLog index 62192f3b93..a4376290ad 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,5 +1,17 @@ 2014-05-14 Joseph Myers + [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 9293f802f8..974d2c80a9 100644 --- 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 a4bf0b88c3..b4177e8f8e 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 98a329179b..c66781c985 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 db9cf337b2..fa35d50fd0 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 a832155d17..38d8110550 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 aea6c51b23..ccef44ac63 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 bbaef680bd..cf9a6db0ef 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 4ba83a46fd..ad8ae9ce5c 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