math: finished cosh.c cleanup

changed the algorithm: large input is not special cased
(when exp(-x) is small compared to exp(x))
and the threshold values are reevaluated
(fdlibm code had a log(2)/2 cutoff for which i could not find
justification, log(2) seems to be a better threshold and this
was verified empirically)

the new code is simpler, makes smaller binaries and should be
faster for common cases

the old comments were removed as they are no longer true for the
new algorithm and the fdlibm copyright was dropped as well
because there is no common code or idea with the original anymore
except for trivial ones.

diff --git a/src/math/cosh.c b/src/math/cosh.c
index 2fcdea8f8e24cec0ddd34802837482d3bb619b6d..100f8231d8829604b9e293491a9abb1774b4b29a 100644 (file)
--- a/src/math/cosh.c
+++ b/src/math/cosh.c
@@ -1,71 +1,40 @@
-/* origin: FreeBSD /usr/src/lib/msun/src/e_cosh.c */
-/*
- * ====================================================
- * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
- *
- * Developed at SunSoft, a Sun Microsystems, Inc. business.
- * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice
- * is preserved.
- * ====================================================
- */
-/* 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) := inf (overflow)
- *
- * Special cases:
- *      cosh(x) is |x| if x is +INF, -INF, or NaN.
- *      only cosh(0)=1 is exact for finite x.
- */
-
 #include "libm.h"
 
+/* cosh(x) = (exp(x) + 1/exp(x))/2
+ *         = 1 + 0.5*(exp(x)-1)*(exp(x)-1)/exp(x)
+ *         = 1 + x*x/2 + o(x^4)
+ */
 double cosh(double x)
 {
        union {double f; uint64_t i;} u = {.f = x};
-       uint32_t ix;
+       uint32_t w;
        double t;
 
        /* |x| */
        u.i &= (uint64_t)-1/2;
        x = u.f;
-       ix = u.i >> 32;
+       w = u.i >> 32;
 
-       /* |x| in [0,0.5*ln2], return 1+expm1(|x|)^2/(2*exp(|x|)) */
-       if (ix < 0x3fd62e43) {
-               t = expm1(x);
-               if (ix < 0x3c800000)
+       /* |x| < log(2) */
+       if (w < 0x3fe62e42) {
+               if (w < 0x3ff00000 - (26<<20)) {
+                       /* raise inexact if x!=0 */
+                       FORCE_EVAL(x + 0x1p120f);
                        return 1;
+               }
+               t = expm1(x);
                return 1 + t*t/(2*(1+t));
        }
 
-       /* |x| in [0.5*ln2,22], return (exp(|x|)+1/exp(|x|))/2; */
-       if (ix < 0x40360000) {
+       /* |x| < log(DBL_MAX) */
+       if (w < 0x40862e42) {
                t = exp(x);
-               return 0.5*t + 0.5/t;
+               /* note: if x>log(0x1p26) then the 1/t is not needed */
+               return 0.5*(t + 1/t);
        }
 
-       /* |x| in [22, log(maxdouble)] return 0.5*exp(|x|) */
-       if (ix < 0x40862e42)
-               return 0.5*exp(x);
-
-       /* |x| in [log(maxdouble), overflowthresold] */
-       if (ix <= 0x408633ce)
-               return __expo2(x);
-
-       /* overflow (or nan) */
-       x *= 0x1p1023;
-       return x;
+       /* |x| > log(DBL_MAX) or nan */
+       /* note: the result is stored to handle overflow */
+       t = __expo2(x);
+       return t;
 }

diff --git a/src/math/coshf.c b/src/math/coshf.c
index 2bed1258a2ca9a87c88d4973c75dbb65c92499af..b09f2ee5751f4341cfb6cc759055279deda1d7cb 100644 (file)
--- a/src/math/coshf.c
+++ b/src/math/coshf.c
@@ -1,54 +1,33 @@
-/* origin: FreeBSD /usr/src/lib/msun/src/e_coshf.c */
-/*
- * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
- */
-/*
- * ====================================================
- * 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.
- * ====================================================
- */
-
 #include "libm.h"
 
 float coshf(float x)
 {
        union {float f; uint32_t i;} u = {.f = x};
-       uint32_t ix;
+       uint32_t w;
        float t;
 
        /* |x| */
        u.i &= 0x7fffffff;
        x = u.f;
-       ix = u.i;
+       w = u.i;
 
-       /* |x| in [0,0.5*ln2], return 1+expm1(|x|)^2/(2*exp(|x|)) */
-       if (ix < 0x3eb17218) {
-               t = expm1f(x);
-               if (ix < 0x39800000)
+       /* |x| < log(2) */
+       if (w < 0x3f317217) {
+               if (w < 0x3f800000 - (12<<23)) {
+                       FORCE_EVAL(x + 0x1p120f);
                        return 1;
+               }
+               t = expm1f(x);
                return 1 + t*t/(2*(1+t));
        }
 
-       /* |x| in [0.5*ln2,9], return (exp(|x|)+1/exp(|x|))/2; */
-       if (ix < 0x41100000) {
+       /* |x| < log(FLT_MAX) */
+       if (w < 0x42b17217) {
                t = expf(x);
-               return 0.5f*t + 0.5f/t;
+               return 0.5f*(t + 1/t);
        }
 
-       /* |x| in [9, log(maxfloat)] return 0.5f*exp(|x|) */
-       if (ix < 0x42b17217)
-               return 0.5f*expf(x);
-
-       /* |x| in [log(maxfloat), overflowthresold] */
-       if (ix <= 0x42b2d4fc)
-               return __expo2f(x);
-
-       /* |x| > overflowthresold or nan */
-       x *= 0x1p127f;
-       return x;
+       /* |x| > log(FLT_MAX) or nan */
+       t = __expo2f(x);
+       return t;
 }

diff --git a/src/math/coshl.c b/src/math/coshl.c
index 3ea56e00bfd5f8be949787119fcef32b90a2f20b..d09070bbf7a4962d7f1f66826ef58cbaafce6158 100644 (file)
--- a/src/math/coshl.c
+++ b/src/math/coshl.c
@@ -1,35 +1,3 @@
-/* origin: OpenBSD /usr/src/lib/libm/src/ld80/e_coshl.c */
-/*
- * ====================================================
- * 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.
- * ====================================================
- */
-/* coshl(x)
- * Method :
- * mathematically coshl(x) if defined to be (exp(x)+exp(-x))/2
- *      1. Replace x by |x| (coshl(x) = coshl(-x)).
- *      2.
- *                                                      [ exp(x) - 1 ]^2
- *          0        <= x <= ln2/2  :  coshl(x) := 1 + -------------------
- *                                                         2*exp(x)
- *
- *                                                 exp(x) +  1/exp(x)
- *          ln2/2    <= x <= 22     :  coshl(x) := -------------------
- *                                                         2
- *          22       <= x <= lnovft :  coshl(x) := expl(x)/2
- *          lnovft   <= x <= ln2ovft:  coshl(x) := expl(x/2)/2 * expl(x/2)
- *          ln2ovft  <  x           :  coshl(x) := inf (overflow)
- *
- * Special cases:
- *      coshl(x) is |x| if x is +INF, -INF, or NaN.
- *      only coshl(0)=1 is exact for finite x.
- */
-
 #include "libm.h"
 
 #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
@@ -45,41 +13,32 @@ long double coshl(long double x)
                struct{uint64_t m; uint16_t se; uint16_t pad;} i;
        } u = {.f = x};
        unsigned ex = u.i.se & 0x7fff;
+       uint32_t w;
        long double t;
-       uint32_t mx,lx;
 
        /* |x| */
        u.i.se = ex;
        x = u.f;
-       mx = u.i.m >> 32;
-       lx = u.i.m;
+       w = u.i.m >> 32;
 
-       /* |x| in [0,0.5*ln2], return 1+expm1l(|x|)^2/(2*expl(|x|)) */
-       if (ex < 0x3fff-2 || (ex == 0x3fff-2 && mx < 0xb17217f7)) {
-               t = expm1l(x);
-               if (ex < 0x3fff-64)
+       /* |x| < log(2) */
+       if (ex < 0x3fff-1 || (ex == 0x3fff-1 && w < 0xb17217f7)) {
+               if (ex < 0x3fff-32) {
+                       FORCE_EVAL(x + 0x1p120f);
                        return 1;
+               }
+               t = expm1l(x);
                return 1 + t*t/(2*(1+t));
        }
 
-       /* |x| in [0.5*ln2,22], return (exp(|x|)+1/exp(|x|)/2; */
-       if (ex < 0x3fff+4 || (ex == 0x3fff+4 && mx < 0xb0000000)) {
+       /* |x| < log(LDBL_MAX) */
+       if (ex < 0x3fff+13 || (ex == 0x3fff+13 && w < 0xb17217f7)) {
                t = expl(x);
-               return 0.5*t + 0.5/t;
-       }
-
-       /* |x| in [22, ln(maxdouble)] return 0.5*exp(|x|) */
-       if (ex < 0x3fff+13 || (ex == 0x3fff+13 && mx < 0xb1700000))
-               return 0.5*expl(x);
-
-       /* |x| in [log(maxdouble), log(2*maxdouble)) */
-       if (ex == 0x3fff+13 && (mx < 0xb174ddc0 ||
-            (mx == 0xb174ddc0 && lx < 0x31aec0eb))) {
-               t = expl(0.5*x);
-               return 0.5*t*t;
+               return 0.5*(t + 1/t);
        }
 
-       /* |x| >= log(2*maxdouble) or nan */
-       return x*0x1p16383L;
+       /* |x| > log(LDBL_MAX) or nan */
+       t = expl(0.5*x);
+       return 0.5*t*t;
 }
 #endif