554edb23ff
Similar to various other bugs in this area, some expm1 implementations do not raise the underflow exception for subnormal arguments, when the result is tiny and inexact. This patch forces the exception in a similar way to previous fixes. (The issue does not apply to the ldbl-* implementations or to those for x86 / x86_64 long double. The change to sysdeps/ieee754/dbl-64/wordsize-64/e_cosh.c is one I missed when previously fixing bug 16354; the bug in that implementation was previously latent, but the expm1 fixes stopped it being latent and so required it to be fixed to avoid spurious underflows from cosh.) Tested for x86_64 and x86. [BZ #16353] * sysdeps/i386/fpu/s_expm1.S (dbl_min): New object. (__expm1): Force underflow exception for arguments with small absolute value. * sysdeps/i386/fpu/s_expm1f.S (flt_min): New object. (__expm1f): Force underflow exception for arguments with small absolute value. * sysdeps/ieee754/dbl-64/s_expm1.c: Include <float.h>. (__expm1): Force underflow exception for arguments with small absolute value. * sysdeps/ieee754/flt-32/s_expm1f.c: Include <float.h>. (__expm1f): Force underflow exception for arguments with small absolute value. * sysdeps/ieee754/dbl-64/wordsize-64/e_cosh.c (__ieee754_cosh): Check for small arguments before calling __expm1. * math/auto-libm-test-in: Do not mark underflow exceptions as possibly missing for bug 16353. * math/auto-libm-test-out: Regenerated.
85 lines
2.2 KiB
C
85 lines
2.2 KiB
C
/* Optimized by Ulrich Drepper <drepper@gmail.com>, 2011 */
|
|
/*
|
|
* ====================================================
|
|
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
|
|
*
|
|
* Developed at SunPro, a Sun Microsystems, Inc. business.
|
|
* Permission to use, copy, modify, and distribute this
|
|
* software is freely granted, provided that this notice
|
|
* is preserved.
|
|
* ====================================================
|
|
*/
|
|
|
|
/* __ieee754_cosh(x)
|
|
* Method :
|
|
* mathematically cosh(x) if defined to be (exp(x)+exp(-x))/2
|
|
* 1. Replace x by |x| (cosh(x) = cosh(-x)).
|
|
* 2.
|
|
* [ exp(x) - 1 ]^2
|
|
* 0 <= x <= ln2/2 : cosh(x) := 1 + -------------------
|
|
* 2*exp(x)
|
|
*
|
|
* exp(x) + 1/exp(x)
|
|
* ln2/2 <= x <= 22 : cosh(x) := -------------------
|
|
* 2
|
|
* 22 <= x <= lnovft : cosh(x) := exp(x)/2
|
|
* lnovft <= x <= ln2ovft: cosh(x) := exp(x/2)/2 * exp(x/2)
|
|
* ln2ovft < x : cosh(x) := huge*huge (overflow)
|
|
*
|
|
* Special cases:
|
|
* cosh(x) is |x| if x is +INF, -INF, or NaN.
|
|
* only cosh(0)=1 is exact for finite x.
|
|
*/
|
|
|
|
#include <math.h>
|
|
#include <math_private.h>
|
|
|
|
static const double one = 1.0, half=0.5, huge = 1.0e300;
|
|
|
|
double
|
|
__ieee754_cosh (double x)
|
|
{
|
|
double t,w;
|
|
int32_t ix;
|
|
|
|
/* High word of |x|. */
|
|
GET_HIGH_WORD(ix,x);
|
|
ix &= 0x7fffffff;
|
|
|
|
/* |x| in [0,22] */
|
|
if (ix < 0x40360000) {
|
|
/* |x| in [0,0.5*ln2], return 1+expm1(|x|)^2/(2*exp(|x|)) */
|
|
if(ix<0x3fd62e43) {
|
|
if (ix<0x3c800000) /* cosh(tiny) = 1 */
|
|
return one;
|
|
t = __expm1(fabs(x));
|
|
w = one+t;
|
|
return one+(t*t)/(w+w);
|
|
}
|
|
|
|
/* |x| in [0.5*ln2,22], return (exp(|x|)+1/exp(|x|)/2; */
|
|
t = __ieee754_exp(fabs(x));
|
|
return half*t+half/t;
|
|
}
|
|
|
|
/* |x| in [22, log(maxdouble)] return half*exp(|x|) */
|
|
if (ix < 0x40862e42) return half*__ieee754_exp(fabs(x));
|
|
|
|
/* |x| in [log(maxdouble), overflowthresold] */
|
|
int64_t fix;
|
|
EXTRACT_WORDS64(fix, x);
|
|
fix &= UINT64_C(0x7fffffffffffffff);
|
|
if (fix <= UINT64_C(0x408633ce8fb9f87d)) {
|
|
w = __ieee754_exp(half*fabs(x));
|
|
t = half*w;
|
|
return t*w;
|
|
}
|
|
|
|
/* x is INF or NaN */
|
|
if(ix>=0x7ff00000) return x*x;
|
|
|
|
/* |x| > overflowthresold, cosh(x) overflow */
|
|
return huge*huge;
|
|
}
|
|
strong_alias (__ieee754_cosh, __cosh_finite)
|