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