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