1 /*- 2 * Copyright (c) 2011 David Schultz 3 * All rights reserved. 4 * 5 * Redistribution and use in source and binary forms, with or without 6 * modification, are permitted provided that the following conditions 7 * are met: 8 * 1. Redistributions of source code must retain the above copyright 9 * notice unmodified, this list of conditions, and the following 10 * disclaimer. 11 * 2. Redistributions in binary form must reproduce the above copyright 12 * notice, this list of conditions and the following disclaimer in the 13 * documentation and/or other materials provided with the distribution. 14 * 15 * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR 16 * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES 17 * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. 18 * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT, 19 * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT 20 * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, 21 * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY 22 * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT 23 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF 24 * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 25 */ 26 27 /* 28 * Hyperbolic tangent of a complex argument z = x + I y. 29 * 30 * The algorithm is from: 31 * 32 * W. Kahan. Branch Cuts for Complex Elementary Functions or Much 33 * Ado About Nothing's Sign Bit. In The State of the Art in 34 * Numerical Analysis, pp. 165 ff. Iserles and Powell, eds., 1987. 35 * 36 * Method: 37 * 38 * Let t = tan(x) 39 * beta = 1/cos^2(y) 40 * s = sinh(x) 41 * rho = cosh(x) 42 * 43 * We have: 44 * 45 * tanh(z) = sinh(z) / cosh(z) 46 * 47 * sinh(x) cos(y) + I cosh(x) sin(y) 48 * = --------------------------------- 49 * cosh(x) cos(y) + I sinh(x) sin(y) 50 * 51 * cosh(x) sinh(x) / cos^2(y) + I tan(y) 52 * = ------------------------------------- 53 * 1 + sinh^2(x) / cos^2(y) 54 * 55 * beta rho s + I t 56 * = ---------------- 57 * 1 + beta s^2 58 * 59 * Modifications: 60 * 61 * I omitted the original algorithm's handling of overflow in tan(x) after 62 * verifying with nearpi.c that this can't happen in IEEE single or double 63 * precision. I also handle large x differently. 64 */ 65 66 #include <sys/cdefs.h> 67 __FBSDID("$FreeBSD$"); 68 69 #include <complex.h> 70 #include <math.h> 71 72 #include "math_private.h" 73 74 double complex 75 ctanh(double complex z) 76 { 77 double x, y; 78 double t, beta, s, rho, denom; 79 uint32_t hx, ix, lx; 80 81 x = creal(z); 82 y = cimag(z); 83 84 EXTRACT_WORDS(hx, lx, x); 85 ix = hx & 0x7fffffff; 86 87 /* 88 * ctanh(NaN +- I 0) = d(NaN) +- I 0 89 * 90 * ctanh(NaN + I y) = d(NaN,y) + I d(NaN,y) for y != 0 91 * 92 * The imaginary part has the sign of x*sin(2*y), but there's no 93 * special effort to get this right. 94 * 95 * ctanh(+-Inf +- I Inf) = +-1 +- I 0 96 * 97 * ctanh(+-Inf + I y) = +-1 + I 0 sin(2y) for y finite 98 * 99 * The imaginary part of the sign is unspecified. This special 100 * case is only needed to avoid a spurious invalid exception when 101 * y is infinite. 102 */ 103 if (ix >= 0x7ff00000) { 104 if ((ix & 0xfffff) | lx) /* x is NaN */ 105 return (CMPLX((x + 0) * (y + 0), 106 y == 0 ? y : (x + 0) * (y + 0))); 107 SET_HIGH_WORD(x, hx - 0x40000000); /* x = copysign(1, x) */ 108 return (CMPLX(x, copysign(0, isinf(y) ? y : sin(y) * cos(y)))); 109 } 110 111 /* 112 * ctanh(x + I NaN) = d(NaN) + I d(NaN) 113 * ctanh(x +- I Inf) = dNaN + I dNaN 114 */ 115 if (!isfinite(y)) 116 return (CMPLX(y - y, y - y)); 117 118 /* 119 * ctanh(+-huge +- I y) ~= +-1 +- I 2sin(2y)/exp(2x), using the 120 * approximation sinh^2(huge) ~= exp(2*huge) / 4. 121 * We use a modified formula to avoid spurious overflow. 122 */ 123 if (ix >= 0x40360000) { /* |x| >= 22 */ 124 double exp_mx = exp(-fabs(x)); 125 return (CMPLX(copysign(1, x), 126 4 * sin(y) * cos(y) * exp_mx * exp_mx)); 127 } 128 129 /* Kahan's algorithm */ 130 t = tan(y); 131 beta = 1.0 + t * t; /* = 1 / cos^2(y) */ 132 s = sinh(x); 133 rho = sqrt(1 + s * s); /* = cosh(x) */ 134 denom = 1 + beta * s * s; 135 return (CMPLX((beta * rho * s) / denom, t / denom)); 136 } 137 138 double complex 139 ctan(double complex z) 140 { 141 142 /* ctan(z) = -I * ctanh(I * z) = I * conj(ctanh(I * conj(z))) */ 143 z = ctanh(CMPLX(cimag(z), creal(z))); 144 return (CMPLX(cimag(z), creal(z))); 145 } 146