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) = NaN + i 0 89 * 90 * ctanh(NaN + i y) = NaN + i NaN 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 +- 0 96 * 97 * ctanh(+-Inf + i y) = +-1 + 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 (cpack(x, (y == 0 ? y : x * y))); 106 SET_HIGH_WORD(x, hx - 0x40000000); /* x = copysign(1, x) */ 107 return (cpack(x, copysign(0, isinf(y) ? y : sin(y) * cos(y)))); 108 } 109 110 /* 111 * ctanh(x + i NAN) = NaN + i NaN 112 * ctanh(x +- i Inf) = NaN + i NaN 113 */ 114 if (!isfinite(y)) 115 return (cpack(y - y, y - y)); 116 117 /* 118 * ctanh(+-huge + i +-y) ~= +-1 +- i 2sin(2y)/exp(2x), using the 119 * approximation sinh^2(huge) ~= exp(2*huge) / 4. 120 * We use a modified formula to avoid spurious overflow. 121 */ 122 if (ix >= 0x40360000) { /* x >= 22 */ 123 double exp_mx = exp(-fabs(x)); 124 return (cpack(copysign(1, x), 125 4 * sin(y) * cos(y) * exp_mx * exp_mx)); 126 } 127 128 /* Kahan's algorithm */ 129 t = tan(y); 130 beta = 1.0 + t * t; /* = 1 / cos^2(y) */ 131 s = sinh(x); 132 rho = sqrt(1 + s * s); /* = cosh(x) */ 133 denom = 1 + beta * s * s; 134 return (cpack((beta * rho * s) / denom, t / denom)); 135 } 136 137 double complex 138 ctan(double complex z) 139 { 140 141 /* ctan(z) = -I * ctanh(I * z) */ 142 z = ctanh(cpack(-cimag(z), creal(z))); 143 return (cpack(cimag(z), -creal(z))); 144 } 145