POK
e_cosh.c
1 /*
2  * POK header
3  *
4  * The following file is a part of the POK project. Any modification should
5  * made according to the POK licence. You CANNOT use this file or a part of
6  * this file is this part of a file for your own project
7  *
8  * For more information on the POK licence, please see our LICENCE FILE
9  *
10  * Please follow the coding guidelines described in doc/CODING_GUIDELINES
11  *
12  * Copyright (c) 2007-2009 POK team
13  *
14  * Created by julien on Fri Jan 30 14:41:34 2009
15  */
16 
17 /* @(#)e_cosh.c 5.1 93/09/24 */
18 /*
19  * ====================================================
20  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
21  *
22  * Developed at SunPro, a Sun Microsystems, Inc. business.
23  * Permission to use, copy, modify, and distribute this
24  * software is freely granted, provided that this notice
25  * is preserved.
26  * ====================================================
27  */
28 
29 /* __ieee754_cosh(x)
30  * Method :
31  * mathematically cosh(x) if defined to be (exp(x)+exp(-x))/2
32  * 1. Replace x by |x| (cosh(x) = cosh(-x)).
33  * 2.
34  * [ exp(x) - 1 ]^2
35  * 0 <= x <= ln2/2 : cosh(x) := 1 + -------------------
36  * 2*exp(x)
37  *
38  * exp(x) + 1/exp(x)
39  * ln2/2 <= x <= 22 : cosh(x) := -------------------
40  * 2
41  * 22 <= x <= lnovft : cosh(x) := exp(x)/2
42  * lnovft <= x <= ln2ovft: cosh(x) := exp(x/2)/2 * exp(x/2)
43  * ln2ovft < x : cosh(x) := huge*huge (overflow)
44  *
45  * Special cases:
46  * cosh(x) is |x| if x is +INF, -INF, or NaN.
47  * only cosh(0)=1 is exact for finite x.
48  */
49 
50 #ifdef POK_NEEDS_LIBMATH
51 #include "math_private.h"
52 #include <libm.h>
53 
54 static const double one = 1.0, half=0.5, huge = 1.0e300;
55 
56 double
57 __ieee754_cosh(double x)
58 {
59  double t,w;
60  int32_t ix;
61  uint32_t lx;
62 
63  /* High word of |x|. */
64  GET_HIGH_WORD(ix,x);
65  ix &= 0x7fffffff;
66 
67  /* x is INF or NaN */
68  if(ix>=0x7ff00000) return x*x;
69 
70  /* |x| in [0,0.5*ln2], return 1+expm1(|x|)^2/(2*exp(|x|)) */
71  if(ix<0x3fd62e43) {
72  t = expm1(fabs(x));
73  w = one+t;
74  if (ix<0x3c800000) return w; /* cosh(tiny) = 1 */
75  return one+(t*t)/(w+w);
76  }
77 
78  /* |x| in [0.5*ln2,22], return (exp(|x|)+1/exp(|x|)/2; */
79  if (ix < 0x40360000) {
80  t = __ieee754_exp(fabs(x));
81  return half*t+half/t;
82  }
83 
84  /* |x| in [22, log(maxdouble)] return half*exp(|x|) */
85  if (ix < 0x40862E42) return half*__ieee754_exp(fabs(x));
86 
87  /* |x| in [log(maxdouble), overflowthresold] */
88  GET_LOW_WORD(lx,x);
89  if (ix<0x408633CE ||
90  ((ix==0x408633ce)&&(lx<=(uint32_t)0x8fb9f87d))) {
91  w = __ieee754_exp(half*fabs(x));
92  t = half*w;
93  return t*w;
94  }
95 
96  /* |x| > overflowthresold, cosh(x) overflow */
97  return huge*huge;
98 }
99 #endif