//  (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