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  #if defined(LIBM_SCCS) && !defined(lint)
      13  static char rcsid[] = "$NetBSD: e_cosh.c,v 1.7 1995/05/10 20:44:58 jtc Exp $";
      14  #endif
      15  
      16  /* __ieee754_coshl(x)
      17   * Method :
      18   * mathematically coshl(x) if defined to be (exp(x)+exp(-x))/2
      19   *	1. Replace x by |x| (coshl(x) = coshl(-x)).
      20   *	2.
      21   *							[ exp(x) - 1 ]^2
      22   *	    0        <= x <= ln2/2  :  coshl(x) := 1 + -------------------
      23   *							   2*exp(x)
      24   *
      25   *						   exp(x) +  1/exp(x)
      26   *	    ln2/2    <= x <= 22     :  coshl(x) := -------------------
      27   *							   2
      28   *	    22       <= x <= lnovft :  coshl(x) := expl(x)/2
      29   *	    lnovft   <= x <= ln2ovft:  coshl(x) := expl(x/2)/2 * expl(x/2)
      30   *	    ln2ovft  <  x	    :  coshl(x) := huge*huge (overflow)
      31   *
      32   * Special cases:
      33   *	coshl(x) is |x| if x is +INF, -INF, or NaN.
      34   *	only coshl(0)=1 is exact for finite x.
      35   */
      36  
      37  #include <math.h>
      38  #include <math_private.h>
      39  #include <libm-alias-finite.h>
      40  
      41  static const long double one = 1.0, half=0.5, huge = 1.0e4900L;
      42  
      43  long double
      44  __ieee754_coshl (long double x)
      45  {
      46  	long double t,w;
      47  	int32_t ex;
      48  	uint32_t mx,lx;
      49  
      50      /* High word of |x|. */
      51  	GET_LDOUBLE_WORDS(ex,mx,lx,x);
      52  	ex &= 0x7fff;
      53  
      54      /* |x| in [0,22] */
      55  	if (ex < 0x4003 || (ex == 0x4003 && mx < 0xb0000000u)) {
      56  	    /* |x| in [0,0.5*ln2], return 1+expm1l(|x|)^2/(2*expl(|x|)) */
      57  		if(ex < 0x3ffd || (ex == 0x3ffd && mx < 0xb17217f7u)) {
      58  		    if (ex<0x3fbc) return one;	/* cosh(tiny) = 1 */
      59  		    t = __expm1l(fabsl(x));
      60  		    w = one+t;
      61  		    return one+(t*t)/(w+w);
      62  		}
      63  
      64  	    /* |x| in [0.5*ln2,22], return (exp(|x|)+1/exp(|x|)/2; */
      65  		t = __ieee754_expl(fabsl(x));
      66  		return half*t+half/t;
      67  	}
      68  
      69      /* |x| in [22, ln(maxdouble)] return half*exp(|x|) */
      70  	if (ex < 0x400c || (ex == 0x400c && mx < 0xb1700000u))
      71  		return half*__ieee754_expl(fabsl(x));
      72  
      73      /* |x| in [log(maxdouble), log(2*maxdouble)) */
      74  	if (ex == 0x400c && (mx < 0xb174ddc0u
      75  			     || (mx == 0xb174ddc0u && lx < 0x31aec0ebu)))
      76  	{
      77  	    w = __ieee754_expl(half*fabsl(x));
      78  	    t = half*w;
      79  	    return t*w;
      80  	}
      81  
      82      /* x is INF or NaN */
      83  	if(ex==0x7fff) return x*x;
      84  
      85      /* |x| >= log(2*maxdouble), cosh(x) overflow */
      86  	return huge*huge;
      87  }
      88  libm_alias_finite (__ieee754_coshl, __coshl)