]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
1 | // (C) Copyright John Maddock 2005. |
2 | // Use, modification and distribution are subject to the | |
3 | // Boost Software License, Version 1.0. (See accompanying file | |
4 | // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) | |
5 | ||
6 | #ifndef BOOST_MATH_COMPLEX_ATANH_INCLUDED | |
7 | #define BOOST_MATH_COMPLEX_ATANH_INCLUDED | |
8 | ||
9 | #ifndef BOOST_MATH_COMPLEX_DETAILS_INCLUDED | |
10 | # include <boost/math/complex/details.hpp> | |
11 | #endif | |
12 | #ifndef BOOST_MATH_LOG1P_INCLUDED | |
13 | # include <boost/math/special_functions/log1p.hpp> | |
14 | #endif | |
15 | #include <boost/assert.hpp> | |
16 | ||
17 | #ifdef BOOST_NO_STDC_NAMESPACE | |
18 | namespace std{ using ::sqrt; using ::fabs; using ::acos; using ::asin; using ::atan; using ::atan2; } | |
19 | #endif | |
20 | ||
21 | namespace boost{ namespace math{ | |
22 | ||
23 | template<class T> | |
24 | std::complex<T> atanh(const std::complex<T>& z) | |
25 | { | |
26 | // | |
27 | // References: | |
28 | // | |
29 | // Eric W. Weisstein. "Inverse Hyperbolic Tangent." | |
30 | // From MathWorld--A Wolfram Web Resource. | |
31 | // http://mathworld.wolfram.com/InverseHyperbolicTangent.html | |
32 | // | |
33 | // Also: The Wolfram Functions Site, | |
34 | // http://functions.wolfram.com/ElementaryFunctions/ArcTanh/ | |
35 | // | |
36 | // Also "Abramowitz and Stegun. Handbook of Mathematical Functions." | |
37 | // at : http://jove.prohosting.com/~skripty/toc.htm | |
38 | // | |
39 | // See also: https://svn.boost.org/trac/boost/ticket/7291 | |
40 | // | |
41 | ||
42 | static const T pi = boost::math::constants::pi<T>(); | |
43 | static const T half_pi = pi / 2; | |
44 | static const T one = static_cast<T>(1.0L); | |
45 | static const T two = static_cast<T>(2.0L); | |
46 | static const T four = static_cast<T>(4.0L); | |
47 | static const T zero = static_cast<T>(0); | |
48 | static const T log_two = boost::math::constants::ln_two<T>(); | |
49 | ||
50 | #ifdef BOOST_MSVC | |
51 | #pragma warning(push) | |
52 | #pragma warning(disable:4127) | |
53 | #endif | |
54 | ||
55 | T x = std::fabs(z.real()); | |
56 | T y = std::fabs(z.imag()); | |
57 | ||
58 | T real, imag; // our results | |
59 | ||
60 | T safe_upper = detail::safe_max(two); | |
61 | T safe_lower = detail::safe_min(static_cast<T>(2)); | |
62 | ||
63 | // | |
64 | // Begin by handling the special cases specified in C99: | |
65 | // | |
66 | if((boost::math::isnan)(x)) | |
67 | { | |
68 | if((boost::math::isnan)(y)) | |
69 | return std::complex<T>(x, x); | |
70 | else if((boost::math::isinf)(y)) | |
71 | return std::complex<T>(0, ((boost::math::signbit)(z.imag()) ? -half_pi : half_pi)); | |
72 | else | |
73 | return std::complex<T>(x, x); | |
74 | } | |
75 | else if((boost::math::isnan)(y)) | |
76 | { | |
77 | if(x == 0) | |
78 | return std::complex<T>(x, y); | |
79 | if((boost::math::isinf)(x)) | |
80 | return std::complex<T>(0, y); | |
81 | else | |
82 | return std::complex<T>(y, y); | |
83 | } | |
84 | else if((x > safe_lower) && (x < safe_upper) && (y > safe_lower) && (y < safe_upper)) | |
85 | { | |
86 | ||
87 | T yy = y*y; | |
88 | T mxm1 = one - x; | |
89 | /// | |
90 | // The real part is given by: | |
91 | // | |
92 | // real(atanh(z)) == log1p(4*x / ((x-1)*(x-1) + y^2)) | |
93 | // | |
94 | real = boost::math::log1p(four * x / (mxm1*mxm1 + yy)); | |
95 | real /= four; | |
96 | if((boost::math::signbit)(z.real())) | |
97 | real = (boost::math::changesign)(real); | |
98 | ||
99 | imag = std::atan2((y * two), (mxm1*(one+x) - yy)); | |
100 | imag /= two; | |
101 | if(z.imag() < 0) | |
102 | imag = (boost::math::changesign)(imag); | |
103 | } | |
104 | else | |
105 | { | |
106 | // | |
107 | // This section handles exception cases that would normally cause | |
108 | // underflow or overflow in the main formulas. | |
109 | // | |
110 | // Begin by working out the real part, we need to approximate | |
111 | // real = boost::math::log1p(4x / ((x-1)^2 + y^2)) | |
112 | // without either overflow or underflow in the squared terms. | |
113 | // | |
114 | T mxm1 = one - x; | |
115 | if(x >= safe_upper) | |
116 | { | |
117 | // x-1 = x to machine precision: | |
118 | if((boost::math::isinf)(x) || (boost::math::isinf)(y)) | |
119 | { | |
120 | real = 0; | |
121 | } | |
122 | else if(y >= safe_upper) | |
123 | { | |
124 | // Big x and y: divide through by x*y: | |
125 | real = boost::math::log1p((four/y) / (x/y + y/x)); | |
126 | } | |
127 | else if(y > one) | |
128 | { | |
129 | // Big x: divide through by x: | |
130 | real = boost::math::log1p(four / (x + y*y/x)); | |
131 | } | |
132 | else | |
133 | { | |
134 | // Big x small y, as above but neglect y^2/x: | |
135 | real = boost::math::log1p(four/x); | |
136 | } | |
137 | } | |
138 | else if(y >= safe_upper) | |
139 | { | |
140 | if(x > one) | |
141 | { | |
142 | // Big y, medium x, divide through by y: | |
143 | real = boost::math::log1p((four*x/y) / (y + mxm1*mxm1/y)); | |
144 | } | |
145 | else | |
146 | { | |
147 | // Small or medium x, large y: | |
148 | real = four*x/y/y; | |
149 | } | |
150 | } | |
151 | else if (x != one) | |
152 | { | |
153 | // y is small, calculate divisor carefully: | |
154 | T div = mxm1*mxm1; | |
155 | if(y > safe_lower) | |
156 | div += y*y; | |
157 | real = boost::math::log1p(four*x/div); | |
158 | } | |
159 | else | |
160 | real = boost::math::changesign(two * (std::log(y) - log_two)); | |
161 | ||
162 | real /= four; | |
163 | if((boost::math::signbit)(z.real())) | |
164 | real = (boost::math::changesign)(real); | |
165 | ||
166 | // | |
167 | // Now handle imaginary part, this is much easier, | |
168 | // if x or y are large, then the formula: | |
169 | // atan2(2y, (1-x)*(1+x) - y^2) | |
170 | // evaluates to +-(PI - theta) where theta is negligible compared to PI. | |
171 | // | |
172 | if((x >= safe_upper) || (y >= safe_upper)) | |
173 | { | |
174 | imag = pi; | |
175 | } | |
176 | else if(x <= safe_lower) | |
177 | { | |
178 | // | |
179 | // If both x and y are small then atan(2y), | |
180 | // otherwise just x^2 is negligible in the divisor: | |
181 | // | |
182 | if(y <= safe_lower) | |
183 | imag = std::atan2(two*y, one); | |
184 | else | |
185 | { | |
186 | if((y == zero) && (x == zero)) | |
187 | imag = 0; | |
188 | else | |
189 | imag = std::atan2(two*y, one - y*y); | |
190 | } | |
191 | } | |
192 | else | |
193 | { | |
194 | // | |
195 | // y^2 is negligible: | |
196 | // | |
197 | if((y == zero) && (x == one)) | |
198 | imag = 0; | |
199 | else | |
200 | imag = std::atan2(two*y, mxm1*(one+x)); | |
201 | } | |
202 | imag /= two; | |
203 | if((boost::math::signbit)(z.imag())) | |
204 | imag = (boost::math::changesign)(imag); | |
205 | } | |
206 | return std::complex<T>(real, imag); | |
207 | #ifdef BOOST_MSVC | |
208 | #pragma warning(pop) | |
209 | #endif | |
210 | } | |
211 | ||
212 | } } // namespaces | |
213 | ||
214 | #endif // BOOST_MATH_COMPLEX_ATANH_INCLUDED |