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 | }
|
---|