15b2ba9d3SPiotr Jasiukajtis /*
25b2ba9d3SPiotr Jasiukajtis * CDDL HEADER START
35b2ba9d3SPiotr Jasiukajtis *
45b2ba9d3SPiotr Jasiukajtis * The contents of this file are subject to the terms of the
55b2ba9d3SPiotr Jasiukajtis * Common Development and Distribution License (the "License").
65b2ba9d3SPiotr Jasiukajtis * You may not use this file except in compliance with the License.
75b2ba9d3SPiotr Jasiukajtis *
85b2ba9d3SPiotr Jasiukajtis * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
95b2ba9d3SPiotr Jasiukajtis * or http://www.opensolaris.org/os/licensing.
105b2ba9d3SPiotr Jasiukajtis * See the License for the specific language governing permissions
115b2ba9d3SPiotr Jasiukajtis * and limitations under the License.
125b2ba9d3SPiotr Jasiukajtis *
135b2ba9d3SPiotr Jasiukajtis * When distributing Covered Code, include this CDDL HEADER in each
145b2ba9d3SPiotr Jasiukajtis * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
155b2ba9d3SPiotr Jasiukajtis * If applicable, add the following below this CDDL HEADER, with the
165b2ba9d3SPiotr Jasiukajtis * fields enclosed by brackets "[]" replaced with your own identifying
175b2ba9d3SPiotr Jasiukajtis * information: Portions Copyright [yyyy] [name of copyright owner]
185b2ba9d3SPiotr Jasiukajtis *
195b2ba9d3SPiotr Jasiukajtis * CDDL HEADER END
205b2ba9d3SPiotr Jasiukajtis */
215b2ba9d3SPiotr Jasiukajtis
225b2ba9d3SPiotr Jasiukajtis /*
235b2ba9d3SPiotr Jasiukajtis * Copyright 2011 Nexenta Systems, Inc. All rights reserved.
245b2ba9d3SPiotr Jasiukajtis */
255b2ba9d3SPiotr Jasiukajtis /*
265b2ba9d3SPiotr Jasiukajtis * Copyright 2006 Sun Microsystems, Inc. All rights reserved.
275b2ba9d3SPiotr Jasiukajtis * Use is subject to license terms.
285b2ba9d3SPiotr Jasiukajtis */
295b2ba9d3SPiotr Jasiukajtis
30*a9d3dcd5SRichard Lowe #pragma weak __ctanh = ctanh
315b2ba9d3SPiotr Jasiukajtis
325b2ba9d3SPiotr Jasiukajtis /* INDENT OFF */
335b2ba9d3SPiotr Jasiukajtis /*
345b2ba9d3SPiotr Jasiukajtis * dcomplex ctanh(dcomplex z);
355b2ba9d3SPiotr Jasiukajtis *
365b2ba9d3SPiotr Jasiukajtis * tanh x + i tan y sinh 2x + i sin 2y
375b2ba9d3SPiotr Jasiukajtis * ctanh z = --------------------- = --------------------
385b2ba9d3SPiotr Jasiukajtis * 1 + i tanh(x)tan(y) cosh 2x + cos 2y
395b2ba9d3SPiotr Jasiukajtis *
405b2ba9d3SPiotr Jasiukajtis * For |x| >= prec/2 (14,28,34,60 for single, double, double extended, quad),
415b2ba9d3SPiotr Jasiukajtis * we use
425b2ba9d3SPiotr Jasiukajtis *
435b2ba9d3SPiotr Jasiukajtis * 1 2x 2 sin 2y
445b2ba9d3SPiotr Jasiukajtis * cosh 2x = sinh 2x = --- e and hence ctanh z = 1 + i -----------;
455b2ba9d3SPiotr Jasiukajtis * 2 2x
465b2ba9d3SPiotr Jasiukajtis * e
475b2ba9d3SPiotr Jasiukajtis *
485b2ba9d3SPiotr Jasiukajtis * otherwise, to avoid cancellation, for |x| < prec/2,
495b2ba9d3SPiotr Jasiukajtis * 2x 2
505b2ba9d3SPiotr Jasiukajtis * (e - 1) 2 2
515b2ba9d3SPiotr Jasiukajtis * cosh 2x + cos 2y = 1 + ------------ + cos y - sin y
525b2ba9d3SPiotr Jasiukajtis * 2x
535b2ba9d3SPiotr Jasiukajtis * 2 e
545b2ba9d3SPiotr Jasiukajtis *
555b2ba9d3SPiotr Jasiukajtis * 1 2x 2 -2x 2
565b2ba9d3SPiotr Jasiukajtis * = --- (e - 1) e + 2 cos y
575b2ba9d3SPiotr Jasiukajtis * 2
585b2ba9d3SPiotr Jasiukajtis * and
595b2ba9d3SPiotr Jasiukajtis *
605b2ba9d3SPiotr Jasiukajtis * [ 2x ]
615b2ba9d3SPiotr Jasiukajtis * 1 [ 2x e - 1 ]
625b2ba9d3SPiotr Jasiukajtis * sinh 2x = --- [ e - 1 + --------- ]
635b2ba9d3SPiotr Jasiukajtis * 2 [ 2x ]
645b2ba9d3SPiotr Jasiukajtis * [ e ]
655b2ba9d3SPiotr Jasiukajtis * 2x
665b2ba9d3SPiotr Jasiukajtis * Implementation notes: let t = expm1(2x) = e - 1, then
675b2ba9d3SPiotr Jasiukajtis *
685b2ba9d3SPiotr Jasiukajtis * 1 [ t*t 2 ] 1 [ t ]
695b2ba9d3SPiotr Jasiukajtis * cosh 2x + cos 2y = --- * [ ----- + 4 cos y ]; sinh 2x = --- * [ t + --- ]
705b2ba9d3SPiotr Jasiukajtis * 2 [ t+1 ] 2 [ t+1 ]
715b2ba9d3SPiotr Jasiukajtis *
725b2ba9d3SPiotr Jasiukajtis * Hence,
735b2ba9d3SPiotr Jasiukajtis *
745b2ba9d3SPiotr Jasiukajtis *
755b2ba9d3SPiotr Jasiukajtis * t*t+2t [4(t+1)(cos y)]*(sin y)
765b2ba9d3SPiotr Jasiukajtis * ctanh z = --------------------------- + i --------------------------
775b2ba9d3SPiotr Jasiukajtis * t*t+[4(t+1)(cos y)](cos y) t*t+[4(t+1)(cos y)](cos y)
785b2ba9d3SPiotr Jasiukajtis *
795b2ba9d3SPiotr Jasiukajtis * EXCEPTION (conform to ISO/IEC 9899:1999(E)):
805b2ba9d3SPiotr Jasiukajtis * ctanh(0,0)=(0,0)
815b2ba9d3SPiotr Jasiukajtis * ctanh(x,inf) = (NaN,NaN) for finite x
825b2ba9d3SPiotr Jasiukajtis * ctanh(x,NaN) = (NaN,NaN) for finite x
835b2ba9d3SPiotr Jasiukajtis * ctanh(inf,y) = 1+ i*0*sin(2y) for positive-signed finite y
845b2ba9d3SPiotr Jasiukajtis * ctanh(inf,inf) = (1, +-0)
855b2ba9d3SPiotr Jasiukajtis * ctanh(inf,NaN) = (1, +-0)
865b2ba9d3SPiotr Jasiukajtis * ctanh(NaN,0) = (NaN,0)
875b2ba9d3SPiotr Jasiukajtis * ctanh(NaN,y) = (NaN,NaN) for non-zero y
885b2ba9d3SPiotr Jasiukajtis * ctanh(NaN,NaN) = (NaN,NaN)
895b2ba9d3SPiotr Jasiukajtis */
905b2ba9d3SPiotr Jasiukajtis /* INDENT ON */
915b2ba9d3SPiotr Jasiukajtis
925b2ba9d3SPiotr Jasiukajtis #include "libm.h" /* exp/expm1/fabs/sin/tanh/sincos */
935b2ba9d3SPiotr Jasiukajtis #include "complex_wrapper.h"
945b2ba9d3SPiotr Jasiukajtis
955b2ba9d3SPiotr Jasiukajtis static const double four = 4.0, two = 2.0, one = 1.0, zero = 0.0;
965b2ba9d3SPiotr Jasiukajtis
975b2ba9d3SPiotr Jasiukajtis dcomplex
ctanh(dcomplex z)985b2ba9d3SPiotr Jasiukajtis ctanh(dcomplex z) {
995b2ba9d3SPiotr Jasiukajtis double t, r, v, u, x, y, S, C;
1005b2ba9d3SPiotr Jasiukajtis int hx, ix, lx, hy, iy, ly;
1015b2ba9d3SPiotr Jasiukajtis dcomplex ans;
1025b2ba9d3SPiotr Jasiukajtis
1035b2ba9d3SPiotr Jasiukajtis x = D_RE(z);
1045b2ba9d3SPiotr Jasiukajtis y = D_IM(z);
1055b2ba9d3SPiotr Jasiukajtis hx = HI_WORD(x);
1065b2ba9d3SPiotr Jasiukajtis lx = LO_WORD(x);
1075b2ba9d3SPiotr Jasiukajtis ix = hx & 0x7fffffff;
1085b2ba9d3SPiotr Jasiukajtis hy = HI_WORD(y);
1095b2ba9d3SPiotr Jasiukajtis ly = LO_WORD(y);
1105b2ba9d3SPiotr Jasiukajtis iy = hy & 0x7fffffff;
1115b2ba9d3SPiotr Jasiukajtis x = fabs(x);
1125b2ba9d3SPiotr Jasiukajtis y = fabs(y);
1135b2ba9d3SPiotr Jasiukajtis
1145b2ba9d3SPiotr Jasiukajtis if ((iy | ly) == 0) { /* ctanh(x,0) = (x,0) for x = 0 or NaN */
1155b2ba9d3SPiotr Jasiukajtis D_RE(ans) = tanh(x);
1165b2ba9d3SPiotr Jasiukajtis D_IM(ans) = zero;
1175b2ba9d3SPiotr Jasiukajtis } else if (iy >= 0x7ff00000) { /* y is inf or NaN */
1185b2ba9d3SPiotr Jasiukajtis if (ix < 0x7ff00000) /* catanh(finite x,inf/nan) is nan */
1195b2ba9d3SPiotr Jasiukajtis D_RE(ans) = D_IM(ans) = y - y;
1205b2ba9d3SPiotr Jasiukajtis else if (((ix - 0x7ff00000) | lx) == 0) { /* x is inf */
1215b2ba9d3SPiotr Jasiukajtis D_RE(ans) = one;
1225b2ba9d3SPiotr Jasiukajtis D_IM(ans) = zero;
1235b2ba9d3SPiotr Jasiukajtis } else {
1245b2ba9d3SPiotr Jasiukajtis D_RE(ans) = x + y;
1255b2ba9d3SPiotr Jasiukajtis D_IM(ans) = y - y;
1265b2ba9d3SPiotr Jasiukajtis }
1275b2ba9d3SPiotr Jasiukajtis } else if (ix >= 0x403c0000) {
1285b2ba9d3SPiotr Jasiukajtis /*
1295b2ba9d3SPiotr Jasiukajtis * |x| > 28 = prec/2 (14,28,34,60)
1305b2ba9d3SPiotr Jasiukajtis * ctanh z ~ 1 + i (sin2y)/(exp(2x))
1315b2ba9d3SPiotr Jasiukajtis */
1325b2ba9d3SPiotr Jasiukajtis D_RE(ans) = one;
1335b2ba9d3SPiotr Jasiukajtis if (iy < 0x7fe00000) /* t = sin(2y) */
1345b2ba9d3SPiotr Jasiukajtis S = sin(y + y);
1355b2ba9d3SPiotr Jasiukajtis else {
1365b2ba9d3SPiotr Jasiukajtis (void) sincos(y, &S, &C);
1375b2ba9d3SPiotr Jasiukajtis S = (S + S) * C;
1385b2ba9d3SPiotr Jasiukajtis }
1395b2ba9d3SPiotr Jasiukajtis if (ix >= 0x7fe00000) { /* |x| > max/2 */
1405b2ba9d3SPiotr Jasiukajtis if (ix >= 0x7ff00000) { /* |x| is inf or NaN */
1415b2ba9d3SPiotr Jasiukajtis if (((ix - 0x7ff00000) | lx) != 0)
1425b2ba9d3SPiotr Jasiukajtis D_RE(ans) = D_IM(ans) = x + y;
1435b2ba9d3SPiotr Jasiukajtis /* x is NaN */
1445b2ba9d3SPiotr Jasiukajtis else
1455b2ba9d3SPiotr Jasiukajtis D_IM(ans) = zero * S; /* x is inf */
1465b2ba9d3SPiotr Jasiukajtis } else
1475b2ba9d3SPiotr Jasiukajtis D_IM(ans) = S * exp(-x); /* underflow */
1485b2ba9d3SPiotr Jasiukajtis } else
1495b2ba9d3SPiotr Jasiukajtis D_IM(ans) = (S + S) * exp(-(x + x));
1505b2ba9d3SPiotr Jasiukajtis /* 2 sin 2y / exp(2x) */
1515b2ba9d3SPiotr Jasiukajtis } else {
1525b2ba9d3SPiotr Jasiukajtis /* INDENT OFF */
1535b2ba9d3SPiotr Jasiukajtis /*
1545b2ba9d3SPiotr Jasiukajtis * t*t+2t
1555b2ba9d3SPiotr Jasiukajtis * ctanh z = --------------------------- +
1565b2ba9d3SPiotr Jasiukajtis * t*t+[4(t+1)(cos y)](cos y)
1575b2ba9d3SPiotr Jasiukajtis *
1585b2ba9d3SPiotr Jasiukajtis * [4(t+1)(cos y)]*(sin y)
1595b2ba9d3SPiotr Jasiukajtis * i --------------------------
1605b2ba9d3SPiotr Jasiukajtis * t*t+[4(t+1)(cos y)](cos y)
1615b2ba9d3SPiotr Jasiukajtis */
1625b2ba9d3SPiotr Jasiukajtis /* INDENT ON */
1635b2ba9d3SPiotr Jasiukajtis (void) sincos(y, &S, &C);
1645b2ba9d3SPiotr Jasiukajtis t = expm1(x + x);
1655b2ba9d3SPiotr Jasiukajtis r = (four * C) * (t + one);
1665b2ba9d3SPiotr Jasiukajtis u = t * t;
1675b2ba9d3SPiotr Jasiukajtis v = one / (u + r * C);
1685b2ba9d3SPiotr Jasiukajtis D_RE(ans) = (u + two * t) * v;
1695b2ba9d3SPiotr Jasiukajtis D_IM(ans) = (r * S) * v;
1705b2ba9d3SPiotr Jasiukajtis }
1715b2ba9d3SPiotr Jasiukajtis if (hx < 0)
1725b2ba9d3SPiotr Jasiukajtis D_RE(ans) = -D_RE(ans);
1735b2ba9d3SPiotr Jasiukajtis if (hy < 0)
1745b2ba9d3SPiotr Jasiukajtis D_IM(ans) = -D_IM(ans);
1755b2ba9d3SPiotr Jasiukajtis return (ans);
1765b2ba9d3SPiotr Jasiukajtis }
177