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/*
13 * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net>
14 *
15 * Permission to use, copy, modify, and distribute this software for any
16 * purpose with or without fee is hereby granted, provided that the above
17 * copyright notice and this permission notice appear in all copies.
18 *
19 * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
20 * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
21 * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
22 * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
23 * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
24 * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
25 * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
26 */
27
28/* coshl(x)
29 * Method :
30 * mathematically coshl(x) if defined to be (exp(x)+exp(-x))/2
31 *      1. Replace x by |x| (coshl(x) = coshl(-x)).
32 *      2.
33 *                                                      [ exp(x) - 1 ]^2
34 *          0        <= x <= ln2/2  :  coshl(x) := 1 + -------------------
35 *                                                         2*exp(x)
36 *
37 *                                                 exp(x) +  1/exp(x)
38 *          ln2/2    <= x <= 22     :  coshl(x) := -------------------
39 *                                                         2
40 *          22       <= x <= lnovft :  coshl(x) := expl(x)/2
41 *          lnovft   <= x <= ln2ovft:  coshl(x) := expl(x/2)/2 * expl(x/2)
42 *          ln2ovft  <  x           :  coshl(x) := huge*huge (overflow)
43 *
44 * Special cases:
45 *      coshl(x) is |x| if x is +INF, -INF, or NaN.
46 *      only coshl(0)=1 is exact for finite x.
47 */
48
49#include <math.h>
50
51#include "math_private.h"
52
53static const long double one = 1.0, half = 0.5, huge = 1.0e4900L,
54ovf_thresh = 1.1357216553474703894801348310092223067821E4L;
55
56long double
57coshl(long double x)
58{
59  long double t, w;
60  int32_t ex;
61  ieee_quad_shape_type u;
62
63  u.value = x;
64  ex = u.parts32.mswhi & 0x7fffffff;
65
66  /* Absolute value of x.  */
67  u.parts32.mswhi = ex;
68
69  /* x is INF or NaN */
70  if (ex >= 0x7fff0000)
71    return x * x;
72
73  /* |x| in [0,0.5*ln2], return 1+expm1l(|x|)^2/(2*expl(|x|)) */
74  if (ex < 0x3ffd62e4) /* 0.3465728759765625 */
75    {
76      t = expm1l (u.value);
77      w = one + t;
78      if (ex < 0x3fb80000) /* |x| < 2^-116 */
79	return w;		/* cosh(tiny) = 1 */
80
81      return one + (t * t) / (w + w);
82    }
83
84  /* |x| in [0.5*ln2,40], return (exp(|x|)+1/exp(|x|)/2; */
85  if (ex < 0x40044000)
86    {
87      t = expl (u.value);
88      return half * t + half / t;
89    }
90
91  /* |x| in [22, ln(maxdouble)] return half*exp(|x|) */
92  if (ex <= 0x400c62e3) /* 11356.375 */
93    return half * expl (u.value);
94
95  /* |x| in [log(maxdouble), overflowthresold] */
96  if (u.value <= ovf_thresh)
97    {
98      w = expl (half * u.value);
99      t = half * w;
100      return t * w;
101    }
102
103  /* |x| > overflowthresold, cosh(x) overflow */
104  return huge * huge;
105}
106DEF_STD(coshl);
107