1226458Sdas/*- 2226458Sdas * Copyright (c) 2011 David Schultz 3226458Sdas * All rights reserved. 4226458Sdas * 5226458Sdas * Redistribution and use in source and binary forms, with or without 6226458Sdas * modification, are permitted provided that the following conditions 7226458Sdas * are met: 8226458Sdas * 1. Redistributions of source code must retain the above copyright 9226458Sdas * notice unmodified, this list of conditions, and the following 10226458Sdas * disclaimer. 11226458Sdas * 2. Redistributions in binary form must reproduce the above copyright 12226458Sdas * notice, this list of conditions and the following disclaimer in the 13226458Sdas * documentation and/or other materials provided with the distribution. 14226458Sdas * 15226458Sdas * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR 16226458Sdas * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES 17226458Sdas * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. 18226458Sdas * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT, 19226458Sdas * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT 20226458Sdas * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, 21226458Sdas * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY 22226458Sdas * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT 23226458Sdas * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF 24226458Sdas * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 25226458Sdas */ 26226458Sdas 27226458Sdas/* 28226458Sdas * Hyperbolic tangent of a complex argument z = x + i y. 29226458Sdas * 30226458Sdas * The algorithm is from: 31226458Sdas * 32226458Sdas * W. Kahan. Branch Cuts for Complex Elementary Functions or Much 33226458Sdas * Ado About Nothing's Sign Bit. In The State of the Art in 34226458Sdas * Numerical Analysis, pp. 165 ff. Iserles and Powell, eds., 1987. 35226458Sdas * 36226458Sdas * Method: 37226458Sdas * 38226458Sdas * Let t = tan(x) 39226458Sdas * beta = 1/cos^2(y) 40226458Sdas * s = sinh(x) 41226458Sdas * rho = cosh(x) 42226458Sdas * 43226458Sdas * We have: 44226458Sdas * 45226458Sdas * tanh(z) = sinh(z) / cosh(z) 46226458Sdas * 47226458Sdas * sinh(x) cos(y) + i cosh(x) sin(y) 48226458Sdas * = --------------------------------- 49226458Sdas * cosh(x) cos(y) + i sinh(x) sin(y) 50226458Sdas * 51226458Sdas * cosh(x) sinh(x) / cos^2(y) + i tan(y) 52226458Sdas * = ------------------------------------- 53226458Sdas * 1 + sinh^2(x) / cos^2(y) 54226458Sdas * 55226458Sdas * beta rho s + i t 56226458Sdas * = ---------------- 57226458Sdas * 1 + beta s^2 58226458Sdas * 59226458Sdas * Modifications: 60226458Sdas * 61226458Sdas * I omitted the original algorithm's handling of overflow in tan(x) after 62226458Sdas * verifying with nearpi.c that this can't happen in IEEE single or double 63226458Sdas * precision. I also handle large x differently. 64226458Sdas */ 65226458Sdas 66226458Sdas#include <sys/cdefs.h> 67226458Sdas__FBSDID("$FreeBSD$"); 68226458Sdas 69226458Sdas#include <complex.h> 70226458Sdas#include <math.h> 71226458Sdas 72226458Sdas#include "math_private.h" 73226458Sdas 74226458Sdasdouble complex 75226458Sdasctanh(double complex z) 76226458Sdas{ 77226458Sdas double x, y; 78226458Sdas double t, beta, s, rho, denom; 79226458Sdas uint32_t hx, ix, lx; 80226458Sdas 81226458Sdas x = creal(z); 82226458Sdas y = cimag(z); 83226458Sdas 84226458Sdas EXTRACT_WORDS(hx, lx, x); 85226458Sdas ix = hx & 0x7fffffff; 86226458Sdas 87226458Sdas /* 88226458Sdas * ctanh(NaN + i 0) = NaN + i 0 89226458Sdas * 90226458Sdas * ctanh(NaN + i y) = NaN + i NaN for y != 0 91226458Sdas * 92226458Sdas * The imaginary part has the sign of x*sin(2*y), but there's no 93226458Sdas * special effort to get this right. 94226458Sdas * 95226458Sdas * ctanh(+-Inf +- i Inf) = +-1 +- 0 96226458Sdas * 97226458Sdas * ctanh(+-Inf + i y) = +-1 + 0 sin(2y) for y finite 98226458Sdas * 99226458Sdas * The imaginary part of the sign is unspecified. This special 100226458Sdas * case is only needed to avoid a spurious invalid exception when 101226458Sdas * y is infinite. 102226458Sdas */ 103226458Sdas if (ix >= 0x7ff00000) { 104226458Sdas if ((ix & 0xfffff) | lx) /* x is NaN */ 105226458Sdas return (cpack(x, (y == 0 ? y : x * y))); 106226458Sdas SET_HIGH_WORD(x, hx - 0x40000000); /* x = copysign(1, x) */ 107226458Sdas return (cpack(x, copysign(0, isinf(y) ? y : sin(y) * cos(y)))); 108226458Sdas } 109226458Sdas 110226458Sdas /* 111226600Sdas * ctanh(x + i NAN) = NaN + i NaN 112226600Sdas * ctanh(x +- i Inf) = NaN + i NaN 113226600Sdas */ 114226600Sdas if (!isfinite(y)) 115226600Sdas return (cpack(y - y, y - y)); 116226600Sdas 117226600Sdas /* 118226458Sdas * ctanh(+-huge + i +-y) ~= +-1 +- i 2sin(2y)/exp(2x), using the 119226458Sdas * approximation sinh^2(huge) ~= exp(2*huge) / 4. 120226458Sdas * We use a modified formula to avoid spurious overflow. 121226458Sdas */ 122226458Sdas if (ix >= 0x40360000) { /* x >= 22 */ 123226458Sdas double exp_mx = exp(-fabs(x)); 124226458Sdas return (cpack(copysign(1, x), 125226458Sdas 4 * sin(y) * cos(y) * exp_mx * exp_mx)); 126226458Sdas } 127226458Sdas 128226458Sdas /* Kahan's algorithm */ 129226458Sdas t = tan(y); 130226458Sdas beta = 1.0 + t * t; /* = 1 / cos^2(y) */ 131226458Sdas s = sinh(x); 132226458Sdas rho = sqrt(1 + s * s); /* = cosh(x) */ 133226458Sdas denom = 1 + beta * s * s; 134226458Sdas return (cpack((beta * rho * s) / denom, t / denom)); 135226458Sdas} 136226458Sdas 137226458Sdasdouble complex 138226458Sdasctan(double complex z) 139226458Sdas{ 140226458Sdas 141226458Sdas /* ctan(z) = -I * ctanh(I * z) */ 142226458Sdas z = ctanh(cpack(-cimag(z), creal(z))); 143226458Sdas return (cpack(cimag(z), -creal(z))); 144226458Sdas} 145