xref: /illumos-gate/usr/src/lib/libm/common/LD/coshl.c (revision ddc0e0b5)
1 /*
2  * CDDL HEADER START
3  *
4  * The contents of this file are subject to the terms of the
5  * Common Development and Distribution License (the "License").
6  * You may not use this file except in compliance with the License.
7  *
8  * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
9  * or http://www.opensolaris.org/os/licensing.
10  * See the License for the specific language governing permissions
11  * and limitations under the License.
12  *
13  * When distributing Covered Code, include this CDDL HEADER in each
14  * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
15  * If applicable, add the following below this CDDL HEADER, with the
16  * fields enclosed by brackets "[]" replaced with your own identifying
17  * information: Portions Copyright [yyyy] [name of copyright owner]
18  *
19  * CDDL HEADER END
20  */
21 
22 /*
23  * Copyright 2011 Nexenta Systems, Inc.  All rights reserved.
24  */
25 /*
26  * Copyright 2006 Sun Microsystems, Inc.  All rights reserved.
27  * Use is subject to license terms.
28  */
29 
30 #pragma weak __coshl = coshl
31 
32 #include "libm.h"
33 #include "longdouble.h"
34 
35 /*
36  * COSH(X)
37  * RETURN THE HYPERBOLIC COSINE OF X
38  *
39  * Method :
40  *	1. Replace x by |x| (COSH(x) = COSH(-x)).
41  *	2.
42  *		                                        [ EXP(x) - 1 ]^2
43  *	    0        <= x <= 0.3465  :  COSH(x) := 1 + -------------------
44  *							   2*EXP(x)
45  *
46  *		                                   EXP(x) +  1/EXP(x)
47  *	    0.3465   <= x <= thresh  :  COSH(x) := -------------------
48  *							   2
49  *	    thresh   <= x <= lnovft  :  COSH(x) := EXP(x)/2
50  *	    lnovft   <= x <  INF     :  COSH(x) := SCALBN(EXP(x-MEP1*ln2),ME)
51  *
52  *
53  * here
54  *	0.3465		a number that is near one half of ln2.
55  *	thresh		a number such that
56  *				EXP(thresh)+EXP(-thresh)=EXP(thresh)
57  *	lnovft		logarithm of the overflow threshold
58  *			= MEP1*ln2 chopped to machine precision.
59  *	ME		maximum exponent
60  *	MEP1		maximum exponent plus 1
61  *
62  * Special cases:
63  *	COSH(x) is |x| if x is +INF, -INF, or NaN.
64  *	only COSH(0)=1 is exact for finite x.
65  */
66 
67 static const long double C[] = {
68 	0.5L,
69 	1.0L,
70 	0.3465L,
71 	45.0L,
72 	1.135652340629414394879149e+04L,
73 	7.004447686242549087858985e-16L,
74 	2.710505431213761085018632e-20L,		/* 2^-65 */
75 };
76 
77 #define	half	C[0]
78 #define	one	C[1]
79 #define	thr1	C[2]
80 #define	thr2	C[3]
81 #define	lnovft	C[4]
82 #define	lnovlo	C[5]
83 #define	tinyl	C[6]
84 
85 long double
coshl(long double x)86 coshl(long double x) {
87 	long double w, t;
88 
89 	w = fabsl(x);
90 	if (!finitel(w))
91 		return (w + w);	/* x is INF or NaN */
92 	if (w < thr1) {
93 		if (w < tinyl)
94 			return (one + w);	/* inexact+directed rounding */
95 		t = expm1l(w);
96 		w = one + t;
97 		w = one + (t * t) / (w + w);
98 		return (w);
99 	}
100 	if (w < thr2) {
101 		t = expl(w);
102 		return (half * (t + one / t));
103 	}
104 	if (w <= lnovft)
105 		return (half * expl(w));
106 	return (scalbnl(expl((w - lnovft) - lnovlo), 16383));
107 }
108