1 /*- 2 * SPDX-License-Identifier: BSD-2-Clause 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 <complex.h> 69 #include <math.h> 70 71 #include "math_private.h" 72 73 double complex 74 ctanh(double complex z) 75 { 76 double x, y; 77 double t, beta, s, rho, denom; 78 uint32_t hx, ix, lx; 79 80 x = creal(z); 81 y = cimag(z); 82 83 EXTRACT_WORDS(hx, lx, x); 84 ix = hx & 0x7fffffff; 85 86 /* 87 * ctanh(NaN +- I 0) = d(NaN) +- I 0 88 * 89 * ctanh(NaN + I y) = d(NaN,y) + I d(NaN,y) for y != 0 90 * 91 * The imaginary part has the sign of x*sin(2*y), but there's no 92 * special effort to get this right. 93 * 94 * ctanh(+-Inf +- I Inf) = +-1 +- I 0 95 * 96 * ctanh(+-Inf + I y) = +-1 + I 0 sin(2y) for y finite 97 * 98 * The imaginary part of the sign is unspecified. This special 99 * case is only needed to avoid a spurious invalid exception when 100 * y is infinite. 101 */ 102 if (ix >= 0x7ff00000) { 103 if ((ix & 0xfffff) | lx) /* x is NaN */ 104 return (CMPLX(nan_mix(x, y), 105 y == 0 ? y : nan_mix(x, y))); 106 SET_HIGH_WORD(x, hx - 0x40000000); /* x = copysign(1, x) */ 107 return (CMPLX(x, copysign(0, isinf(y) ? y : sin(y) * cos(y)))); 108 } 109 110 /* 111 * ctanh(+-0 + i NAN) = +-0 + i NaN 112 * ctanh(+-0 +- i Inf) = +-0 + i NaN 113 * ctanh(x + i NAN) = NaN + i NaN 114 * ctanh(x +- i Inf) = NaN + i NaN 115 */ 116 if (!isfinite(y)) 117 return (CMPLX(x ? y - y : x, y - y)); 118 119 /* 120 * ctanh(+-huge +- I y) ~= +-1 +- I 2sin(2y)/exp(2x), using the 121 * approximation sinh^2(huge) ~= exp(2*huge) / 4. 122 * We use a modified formula to avoid spurious overflow. 123 */ 124 if (ix >= 0x40360000) { /* |x| >= 22 */ 125 double exp_mx = exp(-fabs(x)); 126 return (CMPLX(copysign(1, x), 127 4 * sin(y) * cos(y) * exp_mx * exp_mx)); 128 } 129 130 /* Kahan's algorithm */ 131 t = tan(y); 132 beta = 1.0 + t * t; /* = 1 / cos^2(y) */ 133 s = sinh(x); 134 rho = sqrt(1 + s * s); /* = cosh(x) */ 135 denom = 1 + beta * s * s; 136 return (CMPLX((beta * rho * s) / denom, t / denom)); 137 } 138 139 double complex 140 ctan(double complex z) 141 { 142 143 /* ctan(z) = -I * ctanh(I * z) = I * conj(ctanh(I * conj(z))) */ 144 z = ctanh(CMPLX(cimag(z), creal(z))); 145 return (CMPLX(cimag(z), creal(z))); 146 } 147