|
1 /* E_COSH.C |
|
2 * |
|
3 * Portions Copyright (c) 1993-1999 Nokia Corporation and/or its subsidiary(-ies). |
|
4 * All rights reserved. |
|
5 */ |
|
6 |
|
7 |
|
8 /* @(#)e_cosh.c 5.1 93/09/24 */ |
|
9 /* |
|
10 * ==================================================== |
|
11 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. |
|
12 * |
|
13 * Developed at SunPro, a Sun Microsystems, Inc. business. |
|
14 * Permission to use, copy, modify, and distribute this |
|
15 * software is freely granted, provided that this notice |
|
16 * is preserved. |
|
17 * ==================================================== |
|
18 */ |
|
19 |
|
20 /* __ieee754_cosh(x) |
|
21 * Method : |
|
22 * mathematically cosh(x) if defined to be (exp(x)+exp(-x))/2 |
|
23 * 1. Replace x by |x| (cosh(x) = cosh(-x)). |
|
24 * 2. |
|
25 * [ exp(x) - 1 ]^2 |
|
26 * 0 <= x <= ln2/2 : cosh(x) := 1 + ------------------- |
|
27 * 2*exp(x) |
|
28 * |
|
29 * exp(x) + 1/exp(x) |
|
30 * ln2/2 <= x <= 22 : cosh(x) := ------------------- |
|
31 * 2 |
|
32 * 22 <= x <= lnovft : cosh(x) := exp(x)/2 |
|
33 * lnovft <= x <= ln2ovft: cosh(x) := exp(x/2)/2 * exp(x/2) |
|
34 * ln2ovft < x : cosh(x) := huge*huge (overflow) |
|
35 * |
|
36 * Special cases: |
|
37 * cosh(x) is |x| if x is +INF, -INF, or NaN. |
|
38 * only cosh(0)=1 is exact for finite x. |
|
39 */ |
|
40 |
|
41 #include "FDLIBM.H" |
|
42 |
|
43 static const double one = 1.0, half=0.5, huge = 1.0e300; |
|
44 |
|
45 EXPORT_C double __ieee754_cosh(double x) __SOFTFP |
|
46 { |
|
47 double t,w; |
|
48 __int32_t ix; |
|
49 __uint32_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 = __ieee754_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*__ieee754_exp(fabs(x)); |
|
74 |
|
75 /* |x| in [log(maxdouble), overflowthresold] */ |
|
76 GET_LOW_WORD(lx,x); |
|
77 if (ix<0x408633CE || |
|
78 ((ix==0x408633ce)&&(lx<=(__uint32_t)0x8fb9f87d))) { |
|
79 w = __ieee754_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 } |