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