xref: /illumos-gate/usr/src/lib/libm/common/C/sinh.c (revision 66582b606a8194f7f3ba5b3a3a6dca5b0d346361)
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  * Copyright 2011 Nexenta Systems, Inc.  All rights reserved.
23  */
24 /*
25  * Copyright 2005 Sun Microsystems, Inc.  All rights reserved.
26  * Use is subject to license terms.
27  */
28 
29 #pragma weak __sinh = sinh
30 
31 /* INDENT OFF */
32 /*
33  * sinh(x)
34  * Code originated from 4.3bsd.
35  * Modified by K.C. Ng for SUN 4.0 libm.
36  * Method :
37  *	1. reduce x to non-negative by sinh(-x) = - sinh(x).
38  *	2.
39  *
40  *                                   expm1(x) + expm1(x)/(expm1(x)+1)
41  *   0 <= x <= lnovft   : sinh(x) := --------------------------------
42  *								       2
43  *  lnovft <= x <  INF  : sinh(x) := exp(x-1024*ln2)*2**1023
44  *
45  *
46  * Special cases:
47  *	sinh(x) is x if x is +INF, -INF, or NaN.
48  *	only sinh(0)=0 is exact for finite argument.
49  *
50  */
51 /* INDENT ON */
52 
53 #include "libm.h"
54 
55 static const double
56 	ln2hi = 6.93147180369123816490e-01,
57 	ln2lo = 1.90821492927058770002e-10,
58 	lnovft = 7.09782712893383973096e+02;
59 
60 double
61 sinh(double x) {
62 	double	ox, r, t;
63 
64 	ox = x;
65 	r = fabs(x);
66 	if (!finite(x))
67 		return (x * r);
68 	if (r < lnovft) {
69 		t = expm1(r);
70 		r = copysign((t + t / (1.0 + t)) * 0.5, x);
71 	} else {
72 		if (r < 1000.0)
73 			x = copysign(exp((r - 1024 * ln2hi) - 1024 * ln2lo), x);
74 		r = scalbn(x, 1023);
75 	}
76 	if (!finite(r))
77 		r = _SVID_libm_err(ox, ox, 25);
78 	return (r);
79 }
80