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