Chris@16: //  (C) Copyright John Maddock 2005.
Chris@16: //  Use, modification and distribution are subject to the
Chris@16: //  Boost Software License, Version 1.0. (See accompanying file
Chris@16: //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
Chris@16: 
Chris@16: #ifndef BOOST_MATH_COMPLEX_ATANH_INCLUDED
Chris@16: #define BOOST_MATH_COMPLEX_ATANH_INCLUDED
Chris@16: 
Chris@16: #ifndef BOOST_MATH_COMPLEX_DETAILS_INCLUDED
Chris@16: #  include <boost/math/complex/details.hpp>
Chris@16: #endif
Chris@16: #ifndef BOOST_MATH_LOG1P_INCLUDED
Chris@16: #  include <boost/math/special_functions/log1p.hpp>
Chris@16: #endif
Chris@16: #include <boost/assert.hpp>
Chris@16: 
Chris@16: #ifdef BOOST_NO_STDC_NAMESPACE
Chris@16: namespace std{ using ::sqrt; using ::fabs; using ::acos; using ::asin; using ::atan; using ::atan2; }
Chris@16: #endif
Chris@16: 
Chris@16: namespace boost{ namespace math{
Chris@16: 
Chris@16: template<class T> 
Chris@16: std::complex<T> atanh(const std::complex<T>& z)
Chris@16: {
Chris@16:    //
Chris@16:    // References:
Chris@16:    //
Chris@16:    // Eric W. Weisstein. "Inverse Hyperbolic Tangent." 
Chris@16:    // From MathWorld--A Wolfram Web Resource. 
Chris@16:    // http://mathworld.wolfram.com/InverseHyperbolicTangent.html
Chris@16:    //
Chris@16:    // Also: The Wolfram Functions Site,
Chris@16:    // http://functions.wolfram.com/ElementaryFunctions/ArcTanh/
Chris@16:    //
Chris@16:    // Also "Abramowitz and Stegun. Handbook of Mathematical Functions."
Chris@16:    // at : http://jove.prohosting.com/~skripty/toc.htm
Chris@16:    //
Chris@16:    // See also: https://svn.boost.org/trac/boost/ticket/7291
Chris@16:    //
Chris@16:    
Chris@16:    static const T pi = boost::math::constants::pi<T>();
Chris@16:    static const T half_pi = pi / 2;
Chris@16:    static const T one = static_cast<T>(1.0L);
Chris@16:    static const T two = static_cast<T>(2.0L);
Chris@16:    static const T four = static_cast<T>(4.0L);
Chris@16:    static const T zero = static_cast<T>(0);
Chris@16:    static const T log_two = boost::math::constants::ln_two<T>();
Chris@16: 
Chris@16: #ifdef BOOST_MSVC
Chris@16: #pragma warning(push)
Chris@16: #pragma warning(disable:4127)
Chris@16: #endif
Chris@16: 
Chris@16:    T x = std::fabs(z.real());
Chris@16:    T y = std::fabs(z.imag());
Chris@16: 
Chris@16:    T real, imag;  // our results
Chris@16: 
Chris@16:    T safe_upper = detail::safe_max(two);
Chris@16:    T safe_lower = detail::safe_min(static_cast<T>(2));
Chris@16: 
Chris@16:    //
Chris@16:    // Begin by handling the special cases specified in C99:
Chris@16:    //
Chris@16:    if((boost::math::isnan)(x))
Chris@16:    {
Chris@16:       if((boost::math::isnan)(y))
Chris@16:          return std::complex<T>(x, x);
Chris@16:       else if((boost::math::isinf)(y))
Chris@16:          return std::complex<T>(0, ((boost::math::signbit)(z.imag()) ? -half_pi : half_pi));
Chris@16:       else
Chris@16:          return std::complex<T>(x, x);
Chris@16:    }
Chris@16:    else if((boost::math::isnan)(y))
Chris@16:    {
Chris@16:       if(x == 0)
Chris@16:          return std::complex<T>(x, y);
Chris@16:       if((boost::math::isinf)(x))
Chris@16:          return std::complex<T>(0, y);
Chris@16:       else
Chris@16:          return std::complex<T>(y, y);
Chris@16:    }
Chris@16:    else if((x > safe_lower) && (x < safe_upper) && (y > safe_lower) && (y < safe_upper))
Chris@16:    {
Chris@16: 
Chris@16:       T yy = y*y;
Chris@16:       T mxm1 = one - x;
Chris@16:       ///
Chris@16:       // The real part is given by:
Chris@16:       // 
Chris@16:       // real(atanh(z)) == log1p(4*x / ((x-1)*(x-1) + y^2))
Chris@16:       // 
Chris@16:       real = boost::math::log1p(four * x / (mxm1*mxm1 + yy));
Chris@16:       real /= four;
Chris@16:       if((boost::math::signbit)(z.real()))
Chris@16:          real = (boost::math::changesign)(real);
Chris@16: 
Chris@16:       imag = std::atan2((y * two), (mxm1*(one+x) - yy));
Chris@16:       imag /= two;
Chris@16:       if(z.imag() < 0)
Chris@16:          imag = (boost::math::changesign)(imag);
Chris@16:    }
Chris@16:    else
Chris@16:    {
Chris@16:       //
Chris@16:       // This section handles exception cases that would normally cause
Chris@16:       // underflow or overflow in the main formulas.
Chris@16:       //
Chris@16:       // Begin by working out the real part, we need to approximate
Chris@16:       //    real = boost::math::log1p(4x / ((x-1)^2 + y^2))
Chris@16:       // without either overflow or underflow in the squared terms.
Chris@16:       //
Chris@16:       T mxm1 = one - x;
Chris@16:       if(x >= safe_upper)
Chris@16:       {
Chris@16:          // x-1 = x to machine precision:
Chris@16:          if((boost::math::isinf)(x) || (boost::math::isinf)(y))
Chris@16:          {
Chris@16:             real = 0;
Chris@16:          }
Chris@16:          else if(y >= safe_upper)
Chris@16:          {
Chris@16:             // Big x and y: divide through by x*y:
Chris@16:             real = boost::math::log1p((four/y) / (x/y + y/x));
Chris@16:          }
Chris@16:          else if(y > one)
Chris@16:          {
Chris@16:             // Big x: divide through by x:
Chris@16:             real = boost::math::log1p(four / (x + y*y/x));
Chris@16:          }
Chris@16:          else
Chris@16:          {
Chris@16:             // Big x small y, as above but neglect y^2/x:
Chris@16:             real = boost::math::log1p(four/x);
Chris@16:          }
Chris@16:       }
Chris@16:       else if(y >= safe_upper)
Chris@16:       {
Chris@16:          if(x > one)
Chris@16:          {
Chris@16:             // Big y, medium x, divide through by y:
Chris@16:             real = boost::math::log1p((four*x/y) / (y + mxm1*mxm1/y));
Chris@16:          }
Chris@16:          else
Chris@16:          {
Chris@16:             // Small or medium x, large y:
Chris@16:             real = four*x/y/y;
Chris@16:          }
Chris@16:       }
Chris@16:       else if (x != one)
Chris@16:       {
Chris@16:          // y is small, calculate divisor carefully:
Chris@16:          T div = mxm1*mxm1;
Chris@16:          if(y > safe_lower)
Chris@16:             div += y*y;
Chris@16:          real = boost::math::log1p(four*x/div);
Chris@16:       }
Chris@16:       else
Chris@16:          real = boost::math::changesign(two * (std::log(y) - log_two));
Chris@16: 
Chris@16:       real /= four;
Chris@16:       if((boost::math::signbit)(z.real()))
Chris@16:          real = (boost::math::changesign)(real);
Chris@16: 
Chris@16:       //
Chris@16:       // Now handle imaginary part, this is much easier,
Chris@16:       // if x or y are large, then the formula:
Chris@16:       //    atan2(2y, (1-x)*(1+x) - y^2)
Chris@16:       // evaluates to +-(PI - theta) where theta is negligible compared to PI.
Chris@16:       //
Chris@16:       if((x >= safe_upper) || (y >= safe_upper))
Chris@16:       {
Chris@16:          imag = pi;
Chris@16:       }
Chris@16:       else if(x <= safe_lower)
Chris@16:       {
Chris@16:          //
Chris@16:          // If both x and y are small then atan(2y),
Chris@16:          // otherwise just x^2 is negligible in the divisor:
Chris@16:          //
Chris@16:          if(y <= safe_lower)
Chris@16:             imag = std::atan2(two*y, one);
Chris@16:          else
Chris@16:          {
Chris@16:             if((y == zero) && (x == zero))
Chris@16:                imag = 0;
Chris@16:             else
Chris@16:                imag = std::atan2(two*y, one - y*y);
Chris@16:          }
Chris@16:       }
Chris@16:       else
Chris@16:       {
Chris@16:          //
Chris@16:          // y^2 is negligible:
Chris@16:          //
Chris@16:          if((y == zero) && (x == one))
Chris@16:             imag = 0;
Chris@16:          else
Chris@16:             imag = std::atan2(two*y, mxm1*(one+x));
Chris@16:       }
Chris@16:       imag /= two;
Chris@16:       if((boost::math::signbit)(z.imag()))
Chris@16:          imag = (boost::math::changesign)(imag);
Chris@16:    }
Chris@16:    return std::complex<T>(real, imag);
Chris@16: #ifdef BOOST_MSVC
Chris@16: #pragma warning(pop)
Chris@16: #endif
Chris@16: }
Chris@16: 
Chris@16: } } // namespaces
Chris@16: 
Chris@16: #endif // BOOST_MATH_COMPLEX_ATANH_INCLUDED