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  /* __ieee754_cosh(x)
      14   * Method :
      15   * mathematically cosh(x) if defined to be (exp(x)+exp(-x))/2
      16   *	1. Replace x by |x| (cosh(x) = cosh(-x)).
      17   *	2.
      18   *							[ exp(x) - 1 ]^2
      19   *	    0        <= x <= ln2/2  :  cosh(x) := 1 + -------------------
      20   *							   2*exp(x)
      21   *
      22   *						  exp(x) +  1/exp(x)
      23   *	    ln2/2    <= x <= 40     :  cosh(x) := -------------------
      24   *							  2
      25   *	    40       <= x <= lnovft :  cosh(x) := exp(x)/2
      26   *	    lnovft   <= x <= ln2ovft:  cosh(x) := exp(x/2)/2 * exp(x/2)
      27   *	    ln2ovft  <  x	    :  cosh(x) := huge*huge (overflow)
      28   *
      29   * Special cases:
      30   *	cosh(x) is |x| if x is +INF, -INF, or NaN.
      31   *	only cosh(0)=1 is exact for finite x.
      32   */
      33  
      34  #include <math.h>
      35  #include <math_private.h>
      36  #include <libm-alias-finite.h>
      37  
      38  static const long double one = 1.0L, half=0.5L, huge = 1.0e300L;
      39  
      40  long double
      41  __ieee754_coshl (long double x)
      42  {
      43  	long double t,w;
      44  	int64_t ix;
      45  	double xhi;
      46  
      47      /* High word of |x|. */
      48  	xhi = ldbl_high (x);
      49  	EXTRACT_WORDS64 (ix, xhi);
      50  	ix &= 0x7fffffffffffffffLL;
      51  
      52      /* x is INF or NaN */
      53  	if(ix>=0x7ff0000000000000LL) return x*x;
      54  
      55      /* |x| in [0,0.5*ln2], return 1+expm1(|x|)^2/(2*exp(|x|)) */
      56  	if(ix<0x3fd62e42fefa39efLL) {
      57  	    if (ix<0x3c80000000000000LL) return one;	/* cosh(tiny) = 1 */
      58  	    t = __expm1l(fabsl(x));
      59  	    w = one+t;
      60  	    return one+(t*t)/(w+w);
      61  	}
      62  
      63      /* |x| in [0.5*ln2,40], return (exp(|x|)+1/exp(|x|)/2; */
      64  	if (ix < 0x4044000000000000LL) {
      65  		t = __ieee754_expl(fabsl(x));
      66  		return half*t+half/t;
      67  	}
      68  
      69      /* |x| in [40, log(maxdouble)] return half*exp(|x|) */
      70  	if (ix < 0x40862e42fefa39efLL)  return half*__ieee754_expl(fabsl(x));
      71  
      72      /* |x| in [log(maxdouble), overflowthresold] */
      73  	if (ix < 0x408633ce8fb9f87fLL) {
      74  	    w = __ieee754_expl(half*fabsl(x));
      75  	    t = half*w;
      76  	    return t*w;
      77  	}
      78  
      79      /* |x| > overflowthresold, cosh(x) overflow */
      80  	return huge*huge;
      81  }
      82  libm_alias_finite (__ieee754_coshl, __coshl)