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