Fix loss of precision in cosh and sinh for IBM long double

This commit is contained in:
Andreas Schwab 2012-03-05 20:19:14 +01:00
parent 1c9f702a07
commit bbb78d030e
3 changed files with 16 additions and 10 deletions

View File

@ -1,3 +1,9 @@
2012-03-05 Andreas Schwab <schwab@linux-m68k.org>
* sysdeps/ieee754/ldbl-128ibm/e_coshl.c: Drop exp(-x) term
only for |x| >= 40.
* sysdeps/ieee754/ldbl-128ibm/e_sinhl.c: Likewise.
2012-03-05 H.J. Lu <hongjiu.lu@intel.com>
* sysdeps/unix/sysv/linux/x86_64/gettimeofday.c (gettimeofday_ifunc):

View File

@ -20,9 +20,9 @@
* 2*exp(x)
*
* exp(x) + 1/exp(x)
* ln2/2 <= x <= 22 : cosh(x) := -------------------
* ln2/2 <= x <= 40 : cosh(x) := -------------------
* 2
* 22 <= x <= lnovft : cosh(x) := exp(x)/2
* 40 <= 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)
*
@ -57,13 +57,13 @@ __ieee754_coshl (long double x)
return one+(t*t)/(w+w);
}
/* |x| in [0.5*ln2,22], return (exp(|x|)+1/exp(|x|)/2; */
if (ix < 0x4036000000000000LL) {
/* |x| in [0.5*ln2,40], return (exp(|x|)+1/exp(|x|)/2; */
if (ix < 0x4044000000000000LL) {
t = __ieee754_expl(fabsl(x));
return half*t+half/t;
}
/* |x| in [22, log(maxdouble)] return half*exp(|x|) */
/* |x| in [40, log(maxdouble)] return half*exp(|x|) */
if (ix < 0x40862e42fefa39efLL) return half*__ieee754_expl(fabsl(x));
/* |x| in [log(maxdouble), overflowthresold] */

View File

@ -16,10 +16,10 @@
* 1. Replace x by |x| (sinh(-x) = -sinh(x)).
* 2.
* E + E/(E+1)
* 0 <= x <= 22 : sinh(x) := --------------, E=expm1(x)
* 0 <= x <= 40 : sinh(x) := --------------, E=expm1(x)
* 2
*
* 22 <= x <= lnovft : sinh(x) := exp(x)/2
* 40 <= x <= lnovft : sinh(x) := exp(x)/2
* lnovft <= x <= ln2ovft: sinh(x) := exp(x/2)/2 * exp(x/2)
* ln2ovft < x : sinh(x) := x*shuge (overflow)
*
@ -48,8 +48,8 @@ __ieee754_sinhl(long double x)
h = 0.5;
if (jx<0) h = -h;
/* |x| in [0,22], return sign(x)*0.5*(E+E/(E+1))) */
if (ix < 0x4036000000000000LL) { /* |x|<22 */
/* |x| in [0,40], return sign(x)*0.5*(E+E/(E+1))) */
if (ix < 0x4044000000000000LL) { /* |x|<40 */
if (ix<0x3e20000000000000LL) /* |x|<2**-29 */
if(shuge+x>one) return x;/* sinhl(tiny) = tiny with inexact */
t = __expm1l(fabsl(x));
@ -57,7 +57,7 @@ __ieee754_sinhl(long double x)
return h*(t+t/(t+one));
}
/* |x| in [22, log(maxdouble)] return 0.5*exp(|x|) */
/* |x| in [40, log(maxdouble)] return 0.5*exp(|x|) */
if (ix < 0x40862e42fefa39efLL) return h*__ieee754_expl(fabsl(x));
/* |x| in [log(maxdouble), overflowthresold] */