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 /*
27 * Copyright 2006 Sun Microsystems, Inc. All rights reserved.
28 * Use is subject to license terms.
29 */
30
31 #pragma weak __ctanhf = ctanhf
32
33 #include "libm.h" /* expf/expm1f/fabsf/sincosf/sinf/tanhf */
34 #include "complex_wrapper.h"
35
36 static const float four = 4.0F, two = 2.0F, one = 1.0F, zero = 0.0F;
37
38
39 fcomplex
40 ctanhf(fcomplex z)
41 {
42 float r, u, v, t, x, y, S, C;
43 int hx, ix, hy, iy;
44 fcomplex ans;
45
46 x = F_RE(z);
47 y = F_IM(z);
48 hx = THE_WORD(x);
49 ix = hx & 0x7fffffff;
50 hy = THE_WORD(y);
51 iy = hy & 0x7fffffff;
52 x = fabsf(x);
53 y = fabsf(y);
54
55 if (iy == 0) { /* ctanh(x,0) = (x,0) for x = 0 or NaN */
56 F_RE(ans) = tanhf(x);
57 F_IM(ans) = zero;
58 } else if (iy >= 0x7f800000) { /* y is inf or NaN */
59 if (ix < 0x7f800000) { /* catanh(finite x,inf/nan) is nan */
60 F_RE(ans) = F_IM(ans) = y - y;
61 } else if (ix == 0x7f800000) { /* x is inf */
62 F_RE(ans) = one;
63 F_IM(ans) = zero;
64 } else {
65 F_RE(ans) = x + y;
66 F_IM(ans) = y - y;
67 }
68 } else if (ix >= 0x41600000) {
69 /*
70 * |x| > 14 = prec/2 (14,28,34,60)
71 * ctanh z ~ 1 + i (sin2y)/(exp(2x))
72 */
73 F_RE(ans) = one;
74
75 if (iy < 0x7f000000) { /* t = sin(2y) */
76 S = sinf(y + y);
77 } else {
78 (void) sincosf(y, &S, &C);
79 S = (S + S) * C;
80 }
81
82 if (ix >= 0x7f000000) { /* |x| > max/2 */
83 if (ix >= 0x7f800000) { /* |x| is inf or NaN */
84 if (ix > 0x7f800000) /* x is NaN */
85 F_RE(ans) = F_IM(ans) = x + y;
86 else
87 F_IM(ans) = zero * S; /* x is inf */
88 } else {
89 F_IM(ans) = S * expf(-x); /* underflow */
90 }
91 } else {
92 F_IM(ans) = (S + S) * expf(-(x + x));
93 }
94
95 /* 2 sin 2y / exp(2x) */
96 } else {
97 /* BEGIN CSTYLED */
98 /*
99 * t*t+2t
100 * ctanh z = ---------------------------
101 * t*t+[4(t+1)(cos y)](cos y)
102 *
103 * [4(t+1)(cos y)]*(sin y)
104 * i --------------------------
105 * t*t+[4(t+1)(cos y)](cos y)
106 */
107 /* END CSTYLED */
108 (void) sincosf(y, &S, &C);
109 t = expm1f(x + x);
110 r = (four * C) * (t + one);
111 u = t * t;
112 v = one / (u + r * C);
113 F_RE(ans) = (u + two * t) * v;
114 F_IM(ans) = (r * S) * v;
115 }
116
117 if (hx < 0)
118 F_RE(ans) = -F_RE(ans);
119
120 if (hy < 0)
121 F_IM(ans) = -F_IM(ans);
122
123 return (ans);
124 }