[352] | 1 | /* origin: FreeBSD /usr/src/lib/msun/src/s_ctanh.c */
|
---|
| 2 | /*-
|
---|
| 3 | * Copyright (c) 2011 David Schultz
|
---|
| 4 | * All rights reserved.
|
---|
| 5 | *
|
---|
| 6 | * Redistribution and use in source and binary forms, with or without
|
---|
| 7 | * modification, are permitted provided that the following conditions
|
---|
| 8 | * are met:
|
---|
| 9 | * 1. Redistributions of source code must retain the above copyright
|
---|
| 10 | * notice unmodified, this list of conditions, and the following
|
---|
| 11 | * disclaimer.
|
---|
| 12 | * 2. Redistributions in binary form must reproduce the above copyright
|
---|
| 13 | * notice, this list of conditions and the following disclaimer in the
|
---|
| 14 | * documentation and/or other materials provided with the distribution.
|
---|
| 15 | *
|
---|
| 16 | * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
|
---|
| 17 | * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
|
---|
| 18 | * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
|
---|
| 19 | * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
|
---|
| 20 | * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
|
---|
| 21 | * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
|
---|
| 22 | * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
|
---|
| 23 | * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
|
---|
| 24 | * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
|
---|
| 25 | * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
|
---|
| 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 "libm.h"
|
---|
| 67 |
|
---|
| 68 | double complex ctanh(double complex z)
|
---|
| 69 | {
|
---|
| 70 | double x, y;
|
---|
| 71 | double t, beta, s, rho, denom;
|
---|
| 72 | uint32_t hx, ix, lx;
|
---|
| 73 |
|
---|
| 74 | x = creal(z);
|
---|
| 75 | y = cimag(z);
|
---|
| 76 |
|
---|
| 77 | EXTRACT_WORDS(hx, lx, x);
|
---|
| 78 | ix = hx & 0x7fffffff;
|
---|
| 79 |
|
---|
| 80 | /*
|
---|
| 81 | * ctanh(NaN + i 0) = NaN + i 0
|
---|
| 82 | *
|
---|
| 83 | * ctanh(NaN + i y) = NaN + i NaN for y != 0
|
---|
| 84 | *
|
---|
| 85 | * The imaginary part has the sign of x*sin(2*y), but there's no
|
---|
| 86 | * special effort to get this right.
|
---|
| 87 | *
|
---|
| 88 | * ctanh(+-Inf +- i Inf) = +-1 +- 0
|
---|
| 89 | *
|
---|
| 90 | * ctanh(+-Inf + i y) = +-1 + 0 sin(2y) for y finite
|
---|
| 91 | *
|
---|
| 92 | * The imaginary part of the sign is unspecified. This special
|
---|
| 93 | * case is only needed to avoid a spurious invalid exception when
|
---|
| 94 | * y is infinite.
|
---|
| 95 | */
|
---|
| 96 | if (ix >= 0x7ff00000) {
|
---|
| 97 | if ((ix & 0xfffff) | lx) /* x is NaN */
|
---|
| 98 | return CMPLX(x, (y == 0 ? y : x * y));
|
---|
| 99 | SET_HIGH_WORD(x, hx - 0x40000000); /* x = copysign(1, x) */
|
---|
| 100 | return CMPLX(x, copysign(0, isinf(y) ? y : sin(y) * cos(y)));
|
---|
| 101 | }
|
---|
| 102 |
|
---|
| 103 | /*
|
---|
| 104 | * ctanh(+-0 + i NAN) = +-0 + i NaN
|
---|
| 105 | * ctanh(+-0 +- i Inf) = +-0 + i NaN
|
---|
| 106 | * ctanh(x + i NAN) = NaN + i NaN
|
---|
| 107 | * ctanh(x +- i Inf) = NaN + i NaN
|
---|
| 108 | */
|
---|
| 109 | if (!isfinite(y))
|
---|
| 110 | return CMPLX(x ? y - y : x, y - y);
|
---|
| 111 |
|
---|
| 112 | /*
|
---|
| 113 | * ctanh(+-huge + i +-y) ~= +-1 +- i 2sin(2y)/exp(2x), using the
|
---|
| 114 | * approximation sinh^2(huge) ~= exp(2*huge) / 4.
|
---|
| 115 | * We use a modified formula to avoid spurious overflow.
|
---|
| 116 | */
|
---|
| 117 | if (ix >= 0x40360000) { /* x >= 22 */
|
---|
| 118 | double exp_mx = exp(-fabs(x));
|
---|
| 119 | return CMPLX(copysign(1, x), 4 * sin(y) * cos(y) * exp_mx * exp_mx);
|
---|
| 120 | }
|
---|
| 121 |
|
---|
| 122 | /* Kahan's algorithm */
|
---|
| 123 | t = tan(y);
|
---|
| 124 | beta = 1.0 + t * t; /* = 1 / cos^2(y) */
|
---|
| 125 | s = sinh(x);
|
---|
| 126 | rho = sqrt(1 + s * s); /* = cosh(x) */
|
---|
| 127 | denom = 1 + beta * s * s;
|
---|
| 128 | return CMPLX((beta * rho * s) / denom, t / denom);
|
---|
| 129 | }
|
---|