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)
6 #ifndef BOOST_MATH_COMPLEX_ATANH_INCLUDED
7 #define BOOST_MATH_COMPLEX_ATANH_INCLUDED
9 #ifndef BOOST_MATH_COMPLEX_DETAILS_INCLUDED
10 # include <boost/math/complex/details.hpp>
12 #ifndef BOOST_MATH_LOG1P_INCLUDED
13 # include <boost/math/special_functions/log1p.hpp>
15 #include <boost/assert.hpp>
17 #ifdef BOOST_NO_STDC_NAMESPACE
18 namespace std{ using ::sqrt; using ::fabs; using ::acos; using ::asin; using ::atan; using ::atan2; }
21 namespace boost{ namespace math{
24 std::complex<T> atanh(const std::complex<T>& z)
29 // Eric W. Weisstein. "Inverse Hyperbolic Tangent."
30 // From MathWorld--A Wolfram Web Resource.
31 // http://mathworld.wolfram.com/InverseHyperbolicTangent.html
33 // Also: The Wolfram Functions Site,
34 // http://functions.wolfram.com/ElementaryFunctions/ArcTanh/
36 // Also "Abramowitz and Stegun. Handbook of Mathematical Functions."
37 // at : http://jove.prohosting.com/~skripty/toc.htm
39 // See also: https://svn.boost.org/trac/boost/ticket/7291
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>();
52 #pragma warning(disable:4127)
55 T x = std::fabs(z.real());
56 T y = std::fabs(z.imag());
58 T real, imag; // our results
60 T safe_upper = detail::safe_max(two);
61 T safe_lower = detail::safe_min(static_cast<T>(2));
64 // Begin by handling the special cases specified in C99:
66 if((boost::math::isnan)(x))
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));
73 return std::complex<T>(x, x);
75 else if((boost::math::isnan)(y))
78 return std::complex<T>(x, y);
79 if((boost::math::isinf)(x))
80 return std::complex<T>(0, y);
82 return std::complex<T>(y, y);
84 else if((x > safe_lower) && (x < safe_upper) && (y > safe_lower) && (y < safe_upper))
90 // The real part is given by:
92 // real(atanh(z)) == log1p(4*x / ((x-1)*(x-1) + y^2))
94 real = boost::math::log1p(four * x / (mxm1*mxm1 + yy));
96 if((boost::math::signbit)(z.real()))
97 real = (boost::math::changesign)(real);
99 imag = std::atan2((y * two), (mxm1*(one+x) - yy));
102 imag = (boost::math::changesign)(imag);
107 // This section handles exception cases that would normally cause
108 // underflow or overflow in the main formulas.
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.
117 // x-1 = x to machine precision:
118 if((boost::math::isinf)(x) || (boost::math::isinf)(y))
122 else if(y >= safe_upper)
124 // Big x and y: divide through by x*y:
125 real = boost::math::log1p((four/y) / (x/y + y/x));
129 // Big x: divide through by x:
130 real = boost::math::log1p(four / (x + y*y/x));
134 // Big x small y, as above but neglect y^2/x:
135 real = boost::math::log1p(four/x);
138 else if(y >= safe_upper)
142 // Big y, medium x, divide through by y:
143 real = boost::math::log1p((four*x/y) / (y + mxm1*mxm1/y));
147 // Small or medium x, large y:
153 // y is small, calculate divisor carefully:
157 real = boost::math::log1p(four*x/div);
160 real = boost::math::changesign(two * (std::log(y) - log_two));
163 if((boost::math::signbit)(z.real()))
164 real = (boost::math::changesign)(real);
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.
172 if((x >= safe_upper) || (y >= safe_upper))
176 else if(x <= safe_lower)
179 // If both x and y are small then atan(2y),
180 // otherwise just x^2 is negligible in the divisor:
183 imag = std::atan2(two*y, one);
186 if((y == zero) && (x == zero))
189 imag = std::atan2(two*y, one - y*y);
195 // y^2 is negligible:
197 if((y == zero) && (x == one))
200 imag = std::atan2(two*y, mxm1*(one+x));
203 if((boost::math::signbit)(z.imag()))
204 imag = (boost::math::changesign)(imag);
206 return std::complex<T>(real, imag);
214 #endif // BOOST_MATH_COMPLEX_ATANH_INCLUDED