| Message ID | Pine.LNX.4.64.1406162136430.14606@digraph.polyomino.org.uk |
|---|---|
| State | New |
| Headers | show |
Ping. This patch <https://sourceware.org/ml/libc-alpha/2014-06/msg00396.html> is pending review.
On 06/16/2014 11:38 PM, Joseph S. Myers wrote: > This patch fixes bug 16354, spurious underflows from cosh when a tiny > argument is passed to expm1 and expm1 correctly underflows although > the final result of cosh should be 1. As noted in that bug, some > cases are latent because of expm1 implementations not raising > underflow (bug 16353), but all the implementations are fixed > similarly. They already contained checks for tiny arguments, but the > checks were too late to avoid underflow from expm1 (although they > would avoid underflow from subsequent squaring of the result of > expm1); they are moved before the expm1 calls. > > The thresholds used for considering arguments tiny are not > particularly consistent in how they relate to the precision of the > floating-point format in question. They are, however, all sufficient > to ensure that the round-to-nearest result of cosh is indeed 1 below > the threshold (although sometimes they are smaller than necessary). > But the previous logic did not return 1, but the previously computed 1 > + expm1(abs(x)) value. And the thresholds in the ldbl-128 and > ldbl-128ibm code (0x1p-71L - I suspect 0x3f8b was intended in the code > instead of 0x3fb8 - and (roughly) 0x1p-55L) are not sufficient for > that value to be 1. So by moving the test for tiny arguments, and > consequently returning 1 directly now the expm1 value hasn't been > computed by that point, this patch also fixes bug 17061, the (large > number of ulps) inaccuracy for small arguments in those > implementations. Tests for that bug are duly added. Thanks, Andreas
diff --git a/math/auto-libm-test-in b/math/auto-libm-test-in index 6edad5a..c29abbf 100644 --- a/math/auto-libm-test-in +++ b/math/auto-libm-test-in @@ -624,11 +624,14 @@ cosh 50 # GCC bug 59666: results on directed rounding may be incorrect. cosh max no-test-inline xfail-rounding:ldbl-128ibm cosh -max no-test-inline xfail-rounding:ldbl-128ibm -# Bug 16354: spurious underflow may occur. -cosh min spurious-underflow -cosh -min spurious-underflow -cosh min_subnorm spurious-underflow -cosh -min_subnorm spurious-underflow +cosh min +cosh -min +cosh min_subnorm +cosh -min_subnorm +cosh 0x1p-56 +cosh -0x1p-56 +cosh 0x1p-72 +cosh -0x1p-72 # Test values either side of overflow for each floating-point format. cosh 0x5.96a7ep+4 cosh 0x5.96a7e8p+4 diff --git a/sysdeps/i386/fpu/libm-test-ulps b/sysdeps/i386/fpu/libm-test-ulps index d7424a6..148584e 100644 --- a/sysdeps/i386/fpu/libm-test-ulps +++ b/sysdeps/i386/fpu/libm-test-ulps @@ -895,6 +895,7 @@ ldouble: 2 Function: "cosh_upward": double: 1 +float: 1 idouble: 1 ifloat: 1 ildouble: 2 diff --git a/sysdeps/ieee754/dbl-64/e_cosh.c b/sysdeps/ieee754/dbl-64/e_cosh.c index 6caf943..af3910d 100644 --- a/sysdeps/ieee754/dbl-64/e_cosh.c +++ b/sysdeps/ieee754/dbl-64/e_cosh.c @@ -53,10 +53,10 @@ __ieee754_cosh (double x) /* |x| in [0,0.5*ln2], return 1+expm1(|x|)^2/(2*exp(|x|)) */ if (ix < 0x3fd62e43) { + if (ix < 0x3c800000) + return one; /* cosh(tiny) = 1 */ t = __expm1 (fabs (x)); w = one + t; - if (ix < 0x3c800000) - return w; /* cosh(tiny) = 1 */ return one + (t * t) / (w + w); } diff --git a/sysdeps/ieee754/flt-32/e_coshf.c b/sysdeps/ieee754/flt-32/e_coshf.c index 7eeea4d..dedda47 100644 --- a/sysdeps/ieee754/flt-32/e_coshf.c +++ b/sysdeps/ieee754/flt-32/e_coshf.c @@ -33,9 +33,9 @@ __ieee754_coshf (float x) if (ix < 0x41b00000) { /* |x| in [0,0.5*ln2], return 1+expm1(|x|)^2/(2*exp(|x|)) */ if(ix<0x3eb17218) { + if (ix<0x24000000) return one; /* cosh(tiny) = 1 */ t = __expm1f(fabsf(x)); w = one+t; - if (ix<0x24000000) return w; /* cosh(tiny) = 1 */ return one+(t*t)/(w+w); } diff --git a/sysdeps/ieee754/ldbl-128/e_coshl.c b/sysdeps/ieee754/ldbl-128/e_coshl.c index 59f2030..488c318 100644 --- a/sysdeps/ieee754/ldbl-128/e_coshl.c +++ b/sysdeps/ieee754/ldbl-128/e_coshl.c @@ -77,10 +77,10 @@ __ieee754_coshl (long double x) /* |x| in [0,0.5*ln2], return 1+expm1l(|x|)^2/(2*expl(|x|)) */ if (ex < 0x3ffd62e4) /* 0.3465728759765625 */ { + if (ex < 0x3fb80000) /* |x| < 2^-116 */ + return one; /* cosh(tiny) = 1 */ t = __expm1l (u.value); w = one + t; - if (ex < 0x3fb80000) /* |x| < 2^-116 */ - return w; /* cosh(tiny) = 1 */ return one + (t * t) / (w + w); } diff --git a/sysdeps/ieee754/ldbl-128ibm/e_coshl.c b/sysdeps/ieee754/ldbl-128ibm/e_coshl.c index 92313e2..327b2ab 100644 --- a/sysdeps/ieee754/ldbl-128ibm/e_coshl.c +++ b/sysdeps/ieee754/ldbl-128ibm/e_coshl.c @@ -53,9 +53,9 @@ __ieee754_coshl (long double x) /* |x| in [0,0.5*ln2], return 1+expm1(|x|)^2/(2*exp(|x|)) */ if(ix<0x3fd62e42fefa39efLL) { + if (ix<0x3c80000000000000LL) return one; /* cosh(tiny) = 1 */ t = __expm1l(fabsl(x)); w = one+t; - if (ix<0x3c80000000000000LL) return w; /* cosh(tiny) = 1 */ return one+(t*t)/(w+w); } diff --git a/sysdeps/ieee754/ldbl-96/e_coshl.c b/sysdeps/ieee754/ldbl-96/e_coshl.c index 995c992..dd22cae 100644 --- a/sysdeps/ieee754/ldbl-96/e_coshl.c +++ b/sysdeps/ieee754/ldbl-96/e_coshl.c @@ -54,9 +54,9 @@ __ieee754_coshl (long double x) if (ex < 0x4003 || (ex == 0x4003 && mx < 0xb0000000u)) { /* |x| in [0,0.5*ln2], return 1+expm1l(|x|)^2/(2*expl(|x|)) */ if(ex < 0x3ffd || (ex == 0x3ffd && mx < 0xb17217f7u)) { + if (ex<0x3fbc) return one; /* cosh(tiny) = 1 */ t = __expm1l(fabsl(x)); w = one+t; - if (ex<0x3fbc) return w; /* cosh(tiny) = 1 */ return one+(t*t)/(w+w); }