gonzui

    1: /* @(#)e_cosh.c 5.1 93/09/24 */
    2: /*
    3:  * ====================================================
    4:  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
    5:  *
    6:  * Developed at SunPro, a Sun Microsystems, Inc. business.
    7:  * Permission to use, copy, modify, and distribute this
    8:  * software is freely granted, provided that this notice 
    9:  * is preserved.
   10:  * ====================================================
   11:  */
   12: 
   13: /* LINTLIBRARY */
   14: 
   15: /* cosh(x)
   16:  * Method : 
   17:  * mathematically cosh(x) if defined to be (exp(x)+exp(-x))/2
   18:  *      1. Replace x by |x| (cosh(x) = cosh(-x)). 
   19:  *      2. 
   20:  *                                                      [ exp(x) - 1 ]^2 
   21:  *          0        <= x <= ln2/2  :  cosh(x) := 1 + -------------------
   22:  *                                                            2*exp(x)
   23:  *
   24:  *                                                exp(x) +  1/exp(x)
   25:  *          ln2/2    <= x <= 22     :  cosh(x) := -------------------
   26:  *                                                           2
   27:  *          22       <= x <= lnovft :  cosh(x) := exp(x)/2 
   28:  *          lnovft   <= x <= ln2ovft:  cosh(x) := exp(x/2)/2 * exp(x/2)
   29:  *          ln2ovft  <  x        :  cosh(x) := huge*huge (overflow)
   30:  *
   31:  * Special cases:
   32:  *      cosh(x) is |x| if x is +INF, -INF, or NaN.
   33:  *      only cosh(0)=1 is exact for finite x.
   34:  */
   35: 
   36: #include <sys/cdefs.h>
   37: #include <float.h>
   38: #include <math.h>
   39: 
   40: #include "math_private.h"
   41: 
   42: static const double one = 1.0, half=0.5, huge = 1.0e300;
   43: 
   44: double
   45: cosh(double x)
   46: {       
   47:         double t,w;
   48:         int32_t ix;
   49:         u_int32_t lx;
   50: 
   51:     /* High word of |x|. */
   52:         GET_HIGH_WORD(ix,x);
   53:         ix &= 0x7fffffff;
   54: 
   55:     /* x is INF or NaN */
   56:         if(ix>=0x7ff00000) return x*x; 
   57: 
   58:     /* |x| in [0,0.5*ln2], return 1+expm1(|x|)^2/(2*exp(|x|)) */
   59:         if(ix<0x3fd62e43) {
   60:             t = expm1(fabs(x));
   61:             w = one+t;
   62:             if (ix<0x3c800000) return w;       /* cosh(tiny) = 1 */
   63:             return one+(t*t)/(w+w);
   64:         }
   65: 
   66:     /* |x| in [0.5*ln2,22], return (exp(|x|)+1/exp(|x|)/2; */
   67:         if (ix < 0x40360000) {
   68:                 t = exp(fabs(x));
   69:                 return half*t+half/t;
   70:         }
   71: 
   72:     /* |x| in [22, log(maxdouble)] return half*exp(|x|) */
   73:         if (ix < 0x40862E42)  return half*exp(fabs(x));
   74: 
   75:     /* |x| in [log(maxdouble), overflowthresold] */
   76:         GET_LOW_WORD(lx,x);
   77:         if (ix<0x408633CE || 
   78:               ((ix==0x408633ce)&&(lx<=(u_int32_t)0x8fb9f87d))) {
   79:             w = exp(half*fabs(x));
   80:             t = half*w;
   81:             return t*w;
   82:         }
   83: 
   84:     /* |x| > overflowthresold, cosh(x) overflow */
   85:         return huge*huge;
   86: }
   87: 
   88: #if     LDBL_MANT_DIG == 53
   89: #ifdef  lint
   90: /* PROTOLIB1 */
   91: long double coshl(long double);
   92: #else   /* lint */
   93: __weak_alias(coshl, cosh);
   94: #endif  /* lint */
   95: #endif  /* LDBL_MANT_DIG == 53 */