// (C) Copyright John Maddock 2005. | |
// Use, modification and distribution are subject to the | |
// Boost Software License, Version 1.0. (See accompanying file | |
// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) | |
#ifndef BOOST_MATH_COMPLEX_ATANH_INCLUDED | |
#define BOOST_MATH_COMPLEX_ATANH_INCLUDED | |
#ifndef BOOST_MATH_COMPLEX_DETAILS_INCLUDED | |
# include <boost/math/complex/details.hpp> | |
#endif | |
#ifndef BOOST_MATH_LOG1P_INCLUDED | |
# include <boost/math/special_functions/log1p.hpp> | |
#endif | |
#include <boost/assert.hpp> | |
#ifdef BOOST_NO_STDC_NAMESPACE | |
namespace std{ using ::sqrt; using ::fabs; using ::acos; using ::asin; using ::atan; using ::atan2; } | |
#endif | |
namespace boost{ namespace math{ | |
template<class T> | |
std::complex<T> atanh(const std::complex<T>& z) | |
{ | |
// | |
// References: | |
// | |
// Eric W. Weisstein. "Inverse Hyperbolic Tangent." | |
// From MathWorld--A Wolfram Web Resource. | |
// http://mathworld.wolfram.com/InverseHyperbolicTangent.html | |
// | |
// Also: The Wolfram Functions Site, | |
// http://functions.wolfram.com/ElementaryFunctions/ArcTanh/ | |
// | |
// Also "Abramowitz and Stegun. Handbook of Mathematical Functions." | |
// at : http://jove.prohosting.com/~skripty/toc.htm | |
// | |
static const T half_pi = static_cast<T>(1.57079632679489661923132169163975144L); | |
static const T pi = static_cast<T>(3.141592653589793238462643383279502884197L); | |
static const T one = static_cast<T>(1.0L); | |
static const T two = static_cast<T>(2.0L); | |
static const T four = static_cast<T>(4.0L); | |
static const T zero = static_cast<T>(0); | |
static const T a_crossover = static_cast<T>(0.3L); | |
T x = std::fabs(z.real()); | |
T y = std::fabs(z.imag()); | |
T real, imag; // our results | |
T safe_upper = detail::safe_max(two); | |
T safe_lower = detail::safe_min(static_cast<T>(2)); | |
// | |
// Begin by handling the special cases specified in C99: | |
// | |
if(detail::test_is_nan(x)) | |
{ | |
if(detail::test_is_nan(y)) | |
return std::complex<T>(x, x); | |
else if(std::numeric_limits<T>::has_infinity && (y == std::numeric_limits<T>::infinity())) | |
return std::complex<T>(0, ((z.imag() < 0) ? -half_pi : half_pi)); | |
else | |
return std::complex<T>(x, x); | |
} | |
else if(detail::test_is_nan(y)) | |
{ | |
if(x == 0) | |
return std::complex<T>(x, y); | |
if(std::numeric_limits<T>::has_infinity && (x == std::numeric_limits<T>::infinity())) | |
return std::complex<T>(0, y); | |
else | |
return std::complex<T>(y, y); | |
} | |
else if((x > safe_lower) && (x < safe_upper) && (y > safe_lower) && (y < safe_upper)) | |
{ | |
T xx = x*x; | |
T yy = y*y; | |
T x2 = x * two; | |
/// | |
// The real part is given by: | |
// | |
// real(atanh(z)) == log((1 + x^2 + y^2 + 2x) / (1 + x^2 + y^2 - 2x)) | |
// | |
// However, when x is either large (x > 1/E) or very small | |
// (x < E) then this effectively simplifies | |
// to log(1), leading to wildly inaccurate results. | |
// By dividing the above (top and bottom) by (1 + x^2 + y^2) we get: | |
// | |
// real(atanh(z)) == log((1 + (2x / (1 + x^2 + y^2))) / (1 - (-2x / (1 + x^2 + y^2)))) | |
// | |
// which is much more sensitive to the value of x, when x is not near 1 | |
// (remember we can compute log(1+x) for small x very accurately). | |
// | |
// The cross-over from one method to the other has to be determined | |
// experimentally, the value used below appears correct to within a | |
// factor of 2 (and there are larger errors from other parts | |
// of the input domain anyway). | |
// | |
T alpha = two*x / (one + xx + yy); | |
if(alpha < a_crossover) | |
{ | |
real = boost::math::log1p(alpha) - boost::math::log1p(-alpha); | |
} | |
else | |
{ | |
T xm1 = x - one; | |
real = boost::math::log1p(x2 + xx + yy) - std::log(xm1*xm1 + yy); | |
} | |
real /= four; | |
if(z.real() < 0) | |
real = -real; | |
imag = std::atan2((y * two), (one - xx - yy)); | |
imag /= two; | |
if(z.imag() < 0) | |
imag = -imag; | |
} | |
else | |
{ | |
// | |
// This section handles exception cases that would normally cause | |
// underflow or overflow in the main formulas. | |
// | |
// Begin by working out the real part, we need to approximate | |
// alpha = 2x / (1 + x^2 + y^2) | |
// without either overflow or underflow in the squared terms. | |
// | |
T alpha = 0; | |
if(x >= safe_upper) | |
{ | |
// this is really a test for infinity, | |
// but we may not have the necessary numeric_limits support: | |
if((x > (std::numeric_limits<T>::max)()) || (y > (std::numeric_limits<T>::max)())) | |
{ | |
alpha = 0; | |
} | |
else if(y >= safe_upper) | |
{ | |
// Big x and y: divide alpha through by x*y: | |
alpha = (two/y) / (x/y + y/x); | |
} | |
else if(y > one) | |
{ | |
// Big x: divide through by x: | |
alpha = two / (x + y*y/x); | |
} | |
else | |
{ | |
// Big x small y, as above but neglect y^2/x: | |
alpha = two/x; | |
} | |
} | |
else if(y >= safe_upper) | |
{ | |
if(x > one) | |
{ | |
// Big y, medium x, divide through by y: | |
alpha = (two*x/y) / (y + x*x/y); | |
} | |
else | |
{ | |
// Small x and y, whatever alpha is, it's too small to calculate: | |
alpha = 0; | |
} | |
} | |
else | |
{ | |
// one or both of x and y are small, calculate divisor carefully: | |
T div = one; | |
if(x > safe_lower) | |
div += x*x; | |
if(y > safe_lower) | |
div += y*y; | |
alpha = two*x/div; | |
} | |
if(alpha < a_crossover) | |
{ | |
real = boost::math::log1p(alpha) - boost::math::log1p(-alpha); | |
} | |
else | |
{ | |
// We can only get here as a result of small y and medium sized x, | |
// we can simply neglect the y^2 terms: | |
BOOST_ASSERT(x >= safe_lower); | |
BOOST_ASSERT(x <= safe_upper); | |
//BOOST_ASSERT(y <= safe_lower); | |
T xm1 = x - one; | |
real = std::log(1 + two*x + x*x) - std::log(xm1*xm1); | |
} | |
real /= four; | |
if(z.real() < 0) | |
real = -real; | |
// | |
// Now handle imaginary part, this is much easier, | |
// if x or y are large, then the formula: | |
// atan2(2y, 1 - x^2 - y^2) | |
// evaluates to +-(PI - theta) where theta is negligible compared to PI. | |
// | |
if((x >= safe_upper) || (y >= safe_upper)) | |
{ | |
imag = pi; | |
} | |
else if(x <= safe_lower) | |
{ | |
// | |
// If both x and y are small then atan(2y), | |
// otherwise just x^2 is negligible in the divisor: | |
// | |
if(y <= safe_lower) | |
imag = std::atan2(two*y, one); | |
else | |
{ | |
if((y == zero) && (x == zero)) | |
imag = 0; | |
else | |
imag = std::atan2(two*y, one - y*y); | |
} | |
} | |
else | |
{ | |
// | |
// y^2 is negligible: | |
// | |
if((y == zero) && (x == one)) | |
imag = 0; | |
else | |
imag = std::atan2(two*y, 1 - x*x); | |
} | |
imag /= two; | |
if(z.imag() < 0) | |
imag = -imag; | |
} | |
return std::complex<T>(real, imag); | |
} | |
} } // namespaces | |
#endif // BOOST_MATH_COMPLEX_ATANH_INCLUDED |