xref: /illumos-gate/usr/src/lib/libm/common/C/tanh.c (revision 32991bed)
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 __tanh = tanh
31 
32 /* INDENT OFF */
33 /*
34  * TANH(X)
35  * RETURN THE HYPERBOLIC TANGENT OF X
36  * code based on 4.3bsd
37  * Modified by K.C. Ng for sun 4.0, Jan 31, 1987
38  *
39  * Method :
40  *	1. reduce x to non-negative by tanh(-x) = - tanh(x).
41  *	2.
42  *	    0      <  x <=  1.e-10 :  tanh(x) := x
43  *					          -expm1(-2x)
44  *	    1.e-10 <  x <=  1      :  tanh(x) := --------------
45  *					         expm1(-2x) + 2
46  *							  2
47  *	    1      <= x <=  22.0   :  tanh(x) := 1 -  ---------------
48  *						      expm1(2x) + 2
49  *	    22.0   <  x <= INF     :  tanh(x) := 1.
50  *
51  *	Note: 22 was chosen so that fl(1.0+2/(expm1(2*22)+2)) == 1.
52  *
53  * Special cases:
54  *	tanh(NaN) is NaN;
55  *	only tanh(0)=0 is exact for finite argument.
56  */
57 
58 #include "libm.h"
59 #include "libm_protos.h"
60 #include <math.h>
61 
62 static const double
63 	one = 1.0,
64 	two = 2.0,
65 	small = 1.0e-10,
66 	big = 1.0e10;
67 /* INDENT ON */
68 
69 double
70 tanh(double x) {
71 	double t, y, z;
72 	int signx;
73 	volatile double dummy;
74 
75 	if (isnan(x))
76 		return (x * x);	/* + -> * for Cheetah */
77 	signx = signbit(x);
78 	t = fabs(x);
79 	z = one;
80 	if (t <= 22.0) {
81 		if (t > one)
82 			z = one - two / (expm1(t + t) + two);
83 		else if (t > small) {
84 			y = expm1(-t - t);
85 			z = -y / (y + two);
86 		} else {
87 			/* raise the INEXACT flag for non-zero t */
88 			dummy = t + big;
89 #ifdef lint
90 			dummy = dummy;
91 #endif
92 			return (x);
93 		}
94 	} else if (!finite(t))
95 		return (copysign(1.0, x));
96 	else
97 		return (signx == 1 ? -z + small * small : z - small * small);
98 
99 	return (signx == 1 ? -z : z);
100 }
101