1  /*
       2   * ====================================================
       3   * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
       4   *
       5   * Developed at SunPro, a Sun Microsystems, Inc. business.
       6   * Permission to use, copy, modify, and distribute this
       7   * software is freely granted, provided that this notice
       8   * is preserved.
       9   * ====================================================
      10   */
      11  
      12  /* __ieee754_cosh(x)
      13   * Method :
      14   * mathematically cosh(x) if defined to be (exp(x)+exp(-x))/2
      15   *	1. Replace x by |x| (cosh(x) = cosh(-x)).
      16   *	2.
      17   *							[ exp(x) - 1 ]^2
      18   *	    0        <= x <= ln2/2  :  cosh(x) := 1 + -------------------
      19   *							   2*exp(x)
      20   *
      21   *						  exp(x) +  1/exp(x)
      22   *	    ln2/2    <= x <= 22     :  cosh(x) := -------------------
      23   *							  2
      24   *	    22       <= x <= lnovft :  cosh(x) := exp(x)/2
      25   *	    lnovft   <= x <= ln2ovft:  cosh(x) := exp(x/2)/2 * exp(x/2)
      26   *	    ln2ovft  <  x	    :  cosh(x) := huge*huge (overflow)
      27   *
      28   * Special cases:
      29   *	cosh(x) is |x| if x is +INF, -INF, or NaN.
      30   *	only cosh(0)=1 is exact for finite x.
      31   */
      32  
      33  #include <math.h>
      34  #include <math_private.h>
      35  #include <math-narrow-eval.h>
      36  #include <libm-alias-finite.h>
      37  
      38  static const double one = 1.0, half=0.5, huge = 1.0e300;
      39  
      40  double
      41  __ieee754_cosh (double x)
      42  {
      43  	double t,w;
      44  	int32_t ix;
      45  
      46      /* High word of |x|. */
      47  	GET_HIGH_WORD(ix,x);
      48  	ix &= 0x7fffffff;
      49  
      50      /* |x| in [0,22] */
      51  	if (ix < 0x40360000) {
      52  	    /* |x| in [0,0.5*ln2], return 1+expm1(|x|)^2/(2*exp(|x|)) */
      53  		if(ix<0x3fd62e43) {
      54  		    if (ix<0x3c800000)			/* cosh(tiny) = 1 */
      55  		      return one;
      56  		    t = __expm1(fabs(x));
      57  		    w = one+t;
      58  		    return one+(t*t)/(w+w);
      59  		}
      60  
      61  	    /* |x| in [0.5*ln2,22], return (exp(|x|)+1/exp(|x|)/2; */
      62  		t = __ieee754_exp(fabs(x));
      63  		return half*t+half/t;
      64  	}
      65  
      66      /* |x| in [22, log(maxdouble)] return half*exp(|x|) */
      67  	if (ix < 0x40862e42)  return half*__ieee754_exp(fabs(x));
      68  
      69      /* |x| in [log(maxdouble), overflowthresold] */
      70  	int64_t fix;
      71  	EXTRACT_WORDS64(fix, x);
      72  	fix &= UINT64_C(0x7fffffffffffffff);
      73  	if (fix <= UINT64_C(0x408633ce8fb9f87d)) {
      74  	    w = __ieee754_exp(half*fabs(x));
      75  	    t = half*w;
      76  	    return t*w;
      77  	}
      78  
      79      /* x is INF or NaN */
      80  	if(ix>=0x7ff00000) return x*x;
      81  
      82      /* |x| > overflowthresold, cosh(x) overflow */
      83  	return math_narrow_eval (huge * huge);
      84  }
      85  libm_alias_finite (__ieee754_cosh, __cosh)