1*627f7eb2Smrg /*
2*627f7eb2Smrg * ====================================================
3*627f7eb2Smrg * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
4*627f7eb2Smrg *
5*627f7eb2Smrg * Developed at SunPro, a Sun Microsystems, Inc. business.
6*627f7eb2Smrg * Permission to use, copy, modify, and distribute this
7*627f7eb2Smrg * software is freely granted, provided that this notice
8*627f7eb2Smrg * is preserved.
9*627f7eb2Smrg * ====================================================
10*627f7eb2Smrg */
11*627f7eb2Smrg
12*627f7eb2Smrg /* Changes for 128-bit long double are
13*627f7eb2Smrg Copyright (C) 2001 Stephen L. Moshier <moshier@na-net.ornl.gov>
14*627f7eb2Smrg and are incorporated herein by permission of the author. The author
15*627f7eb2Smrg reserves the right to distribute this material elsewhere under different
16*627f7eb2Smrg copying permissions. These modifications are distributed here under
17*627f7eb2Smrg the following terms:
18*627f7eb2Smrg
19*627f7eb2Smrg This library is free software; you can redistribute it and/or
20*627f7eb2Smrg modify it under the terms of the GNU Lesser General Public
21*627f7eb2Smrg License as published by the Free Software Foundation; either
22*627f7eb2Smrg version 2.1 of the License, or (at your option) any later version.
23*627f7eb2Smrg
24*627f7eb2Smrg This library is distributed in the hope that it will be useful,
25*627f7eb2Smrg but WITHOUT ANY WARRANTY; without even the implied warranty of
26*627f7eb2Smrg MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
27*627f7eb2Smrg Lesser General Public License for more details.
28*627f7eb2Smrg
29*627f7eb2Smrg You should have received a copy of the GNU Lesser General Public
30*627f7eb2Smrg License along with this library; if not, see
31*627f7eb2Smrg <http://www.gnu.org/licenses/>. */
32*627f7eb2Smrg
33*627f7eb2Smrg /* coshq(x)
34*627f7eb2Smrg * Method :
35*627f7eb2Smrg * mathematically coshl(x) if defined to be (exp(x)+exp(-x))/2
36*627f7eb2Smrg * 1. Replace x by |x| (coshl(x) = coshl(-x)).
37*627f7eb2Smrg * 2.
38*627f7eb2Smrg * [ exp(x) - 1 ]^2
39*627f7eb2Smrg * 0 <= x <= ln2/2 : coshl(x) := 1 + -------------------
40*627f7eb2Smrg * 2*exp(x)
41*627f7eb2Smrg *
42*627f7eb2Smrg * exp(x) + 1/exp(x)
43*627f7eb2Smrg * ln2/2 <= x <= 22 : coshl(x) := -------------------
44*627f7eb2Smrg * 2
45*627f7eb2Smrg * 22 <= x <= lnovft : coshl(x) := expq(x)/2
46*627f7eb2Smrg * lnovft <= x <= ln2ovft: coshl(x) := expq(x/2)/2 * expq(x/2)
47*627f7eb2Smrg * ln2ovft < x : coshl(x) := huge*huge (overflow)
48*627f7eb2Smrg *
49*627f7eb2Smrg * Special cases:
50*627f7eb2Smrg * coshl(x) is |x| if x is +INF, -INF, or NaN.
51*627f7eb2Smrg * only coshl(0)=1 is exact for finite x.
52*627f7eb2Smrg */
53*627f7eb2Smrg
54*627f7eb2Smrg #include "quadmath-imp.h"
55*627f7eb2Smrg
56*627f7eb2Smrg static const __float128 one = 1.0, half = 0.5, huge = 1.0e4900Q,
57*627f7eb2Smrg ovf_thresh = 1.1357216553474703894801348310092223067821E4Q;
58*627f7eb2Smrg
59*627f7eb2Smrg __float128
coshq(__float128 x)60*627f7eb2Smrg coshq (__float128 x)
61*627f7eb2Smrg {
62*627f7eb2Smrg __float128 t, w;
63*627f7eb2Smrg int32_t ex;
64*627f7eb2Smrg ieee854_float128 u;
65*627f7eb2Smrg
66*627f7eb2Smrg u.value = x;
67*627f7eb2Smrg ex = u.words32.w0 & 0x7fffffff;
68*627f7eb2Smrg
69*627f7eb2Smrg /* Absolute value of x. */
70*627f7eb2Smrg u.words32.w0 = ex;
71*627f7eb2Smrg
72*627f7eb2Smrg /* x is INF or NaN */
73*627f7eb2Smrg if (ex >= 0x7fff0000)
74*627f7eb2Smrg return x * x;
75*627f7eb2Smrg
76*627f7eb2Smrg /* |x| in [0,0.5*ln2], return 1+expm1q(|x|)^2/(2*expq(|x|)) */
77*627f7eb2Smrg if (ex < 0x3ffd62e4) /* 0.3465728759765625 */
78*627f7eb2Smrg {
79*627f7eb2Smrg if (ex < 0x3fb80000) /* |x| < 2^-116 */
80*627f7eb2Smrg return one; /* cosh(tiny) = 1 */
81*627f7eb2Smrg t = expm1q (u.value);
82*627f7eb2Smrg w = one + t;
83*627f7eb2Smrg
84*627f7eb2Smrg return one + (t * t) / (w + w);
85*627f7eb2Smrg }
86*627f7eb2Smrg
87*627f7eb2Smrg /* |x| in [0.5*ln2,40], return (exp(|x|)+1/exp(|x|)/2; */
88*627f7eb2Smrg if (ex < 0x40044000)
89*627f7eb2Smrg {
90*627f7eb2Smrg t = expq (u.value);
91*627f7eb2Smrg return half * t + half / t;
92*627f7eb2Smrg }
93*627f7eb2Smrg
94*627f7eb2Smrg /* |x| in [22, ln(maxdouble)] return half*exp(|x|) */
95*627f7eb2Smrg if (ex <= 0x400c62e3) /* 11356.375 */
96*627f7eb2Smrg return half * expq (u.value);
97*627f7eb2Smrg
98*627f7eb2Smrg /* |x| in [log(maxdouble), overflowthresold] */
99*627f7eb2Smrg if (u.value <= ovf_thresh)
100*627f7eb2Smrg {
101*627f7eb2Smrg w = expq (half * u.value);
102*627f7eb2Smrg t = half * w;
103*627f7eb2Smrg return t * w;
104*627f7eb2Smrg }
105*627f7eb2Smrg
106*627f7eb2Smrg /* |x| > overflowthresold, cosh(x) overflow */
107*627f7eb2Smrg return huge * huge;
108*627f7eb2Smrg }
109