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 __ctanhl = ctanhl 31 32 #include "libm.h" /* expl/expm1l/fabsl/isinfl/isnanl/sincosl/sinl/tanhl */ 33 #include "complex_wrapper.h" 34 #include "longdouble.h" 35 36 /* INDENT OFF */ 37 static const long double four = 4.0L, two = 2.0L, one = 1.0L, zero = 0.0L; 38 /* INDENT ON */ 39 40 ldcomplex 41 ctanhl(ldcomplex z) { 42 long double r, u, v, t, x, y, S, C; 43 int hx, ix, hy, iy; 44 ldcomplex ans; 45 46 x = LD_RE(z); 47 y = LD_IM(z); 48 hx = HI_XWORD(x); 49 ix = hx & 0x7fffffff; 50 hy = HI_XWORD(y); 51 iy = hy & 0x7fffffff; 52 x = fabsl(x); 53 y = fabsl(y); 54 55 if (y == zero) { /* ctanh(x,0) = (x,0) for x = 0 or NaN */ 56 LD_RE(ans) = tanhl(x); 57 LD_IM(ans) = zero; 58 } else if (iy >= 0x7fff0000) { /* y is inf or NaN */ 59 if (ix < 0x7fff0000) /* catanh(finite x,inf/nan) is nan */ 60 LD_RE(ans) = LD_IM(ans) = y - y; 61 else if (isinfl(x)) { /* x is inf */ 62 LD_RE(ans) = one; 63 LD_IM(ans) = zero; 64 } else { 65 LD_RE(ans) = x + y; 66 LD_IM(ans) = y - y; 67 } 68 } else if (ix >= 0x4004e000) { 69 /* INDENT OFF */ 70 /* 71 * |x| > 60 = prec/2 (14,28,34,60) 72 * ctanh z ~ 1 + i (sin2y)/(exp(2x)) 73 */ 74 /* INDENT ON */ 75 LD_RE(ans) = one; 76 if (iy < 0x7ffe0000) /* t = sin(2y) */ 77 S = sinl(y + y); 78 else { 79 (void) sincosl(y, &S, &C); 80 S = (S + S) * C; 81 } 82 if (ix >= 0x7ffe0000) { /* |x| > max/2 */ 83 if (ix >= 0x7fff0000) { /* |x| is inf or NaN */ 84 if (isnanl(x)) /* x is NaN */ 85 LD_RE(ans) = LD_IM(ans) = x + y; 86 else 87 LD_IM(ans) = zero * S; /* x is inf */ 88 } else 89 LD_IM(ans) = S * expl(-x); /* underflow */ 90 } else 91 LD_IM(ans) = (S + S) * expl(-(x + x)); 92 /* 2 sin 2y / exp(2x) */ 93 } else { 94 /* INDENT OFF */ 95 /* 96 * t*t+2t 97 * ctanh z = --------------------------- 98 * t*t+[4(t+1)(cos y)](cos y) 99 * 100 * [4(t+1)(cos y)]*(sin y) 101 * i -------------------------- 102 * t*t+[4(t+1)(cos y)](cos y) 103 */ 104 /* INDENT ON */ 105 sincosl(y, &S, &C); 106 t = expm1l(x + x); 107 r = (four * C) * (t + one); 108 u = t * t; 109 v = one / (u + r * C); 110 LD_RE(ans) = (u + two * t) * v; 111 LD_IM(ans) = (r * S) * v; 112 } 113 if (hx < 0) 114 LD_RE(ans) = -LD_RE(ans); 115 if (hy < 0) 116 LD_IM(ans) = -LD_IM(ans); 117 return (ans); 118 } 119