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

diff --git a/ChangeLog b/ChangeLog
index 6875313..5c9871b 100644 (file)
--- a/ChangeLog
+++ b/ChangeLog
@@ -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):

diff --git a/sysdeps/ieee754/ldbl-128ibm/e_coshl.c b/sysdeps/ieee754/ldbl-128ibm/e_coshl.c
index ebc9436..569b841 100644 (file)
--- a/sysdeps/ieee754/ldbl-128ibm/e_coshl.c
+++ b/sysdeps/ieee754/ldbl-128ibm/e_coshl.c
@@ -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] */

diff --git a/sysdeps/ieee754/ldbl-128ibm/e_sinhl.c b/sysdeps/ieee754/ldbl-128ibm/e_sinhl.c
index b8e86b1..4b53d90 100644 (file)
--- a/sysdeps/ieee754/ldbl-128ibm/e_sinhl.c
+++ b/sysdeps/ieee754/ldbl-128ibm/e_sinhl.c
@@ -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] */