//  (C) Copyright John Maddock 2005.

 //  Distributed under 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_ACOS_INCLUDED

 #define BOOST_MATH_COMPLEX_ACOS_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> acos(const std::complex<T>& z)

 {

    //

    // This implementation is a transcription of the pseudo-code in:

    //

    // "Implementing the Complex Arcsine and Arccosine Functions using Exception Handling."

    // T E Hull, Thomas F Fairgrieve and Ping Tak Peter Tang.

    // ACM Transactions on Mathematical Software, Vol 23, No 3, Sept 1997.

    //

    //

    // These static constants should really be in a maths constants library:

    //

    static const T one = static_cast<T>(1);

    //static const T two = static_cast<T>(2);

    static const T half = static_cast<T>(0.5L);

    static const T a_crossover = static_cast<T>(1.5L);

    static const T b_crossover = static_cast<T>(0.6417L);

    static const T s_pi = static_cast<T>(3.141592653589793238462643383279502884197L);

    static const T half_pi = static_cast<T>(1.57079632679489661923132169163975144L);

    static const T log_two = static_cast<T>(0.69314718055994530941723212145817657L);

    static const T quarter_pi = static_cast<T>(0.78539816339744830961566084581987572L);

    //

    // Get real and imaginary parts, discard the signs as we can 

    // figure out the sign of the result later:

    //

    T x = std::fabs(z.real());

    T y = std::fabs(z.imag());

    T real, imag; // these hold our result

    // 

    // Handle special cases specified by the C99 standard,

    // many of these special cases aren't really needed here,

    // but doing it this way prevents overflow/underflow arithmetic

    // in the main body of the logic, which may trip up some machines:

    //

    if(std::numeric_limits<T>::has_infinity && (x == std::numeric_limits<T>::infinity()))

    {

       if(y == std::numeric_limits<T>::infinity())

       {

          real = quarter_pi;

          imag = std::numeric_limits<T>::infinity();

       }

       else if(detail::test_is_nan(y))

       {

          return std::complex<T>(y, -std::numeric_limits<T>::infinity());

       }

       else

       {

          // y is not infinity or nan:

          real = 0;

          imag = std::numeric_limits<T>::infinity();

       }

    }

    else if(detail::test_is_nan(x))

    {

       if(y == std::numeric_limits<T>::infinity())

          return std::complex<T>(x, (z.imag() < 0) ? std::numeric_limits<T>::infinity() :  -std::numeric_limits<T>::infinity());

       return std::complex<T>(x, x);

    }

    else if(std::numeric_limits<T>::has_infinity && (y == std::numeric_limits<T>::infinity()))

    {

       real = half_pi;

       imag = std::numeric_limits<T>::infinity();

    }

    else if(detail::test_is_nan(y))

    {

       return std::complex<T>((x == 0) ? half_pi : y, y);

    }

    else

    {

       //

       // What follows is the regular Hull et al code,

       // begin with the special case for real numbers:

       //

       if((y == 0) && (x <= one))

          return std::complex<T>((x == 0) ? half_pi : std::acos(z.real()));

       //

       // Figure out if our input is within the "safe area" identified by Hull et al.

       // This would be more efficient with portable floating point exception handling;

       // fortunately the quantities M and u identified by Hull et al (figure 3), 

       // match with the max and min methods of numeric_limits<T>.

       //

       T safe_max = detail::safe_max(static_cast<T>(8));

       T safe_min = detail::safe_min(static_cast<T>(4));

       T xp1 = one + x;

       T xm1 = x - one;

       if((x < safe_max) && (x > safe_min) && (y < safe_max) && (y > safe_min))

       {

          T yy = y * y;

          T r = std::sqrt(xp1*xp1 + yy);

          T s = std::sqrt(xm1*xm1 + yy);

          T a = half * (r + s);

          T b = x / a;

          if(b <= b_crossover)

          {

             real = std::acos(b);

          }

          else

          {

             T apx = a + x;

             if(x <= one)

             {

                real = std::atan(std::sqrt(half * apx * (yy /(r + xp1) + (s-xm1)))/x);

             }

             else

             {

                real = std::atan((y * std::sqrt(half * (apx/(r + xp1) + apx/(s+xm1))))/x);

             }

          }

          if(a <= a_crossover)

          {

             T am1;

             if(x < one)

             {

                am1 = half * (yy/(r + xp1) + yy/(s - xm1));

             }

             else

             {

                am1 = half * (yy/(r + xp1) + (s + xm1));

             }

             imag = boost::math::log1p(am1 + std::sqrt(am1 * (a + one)));

          }

          else

          {

             imag = std::log(a + std::sqrt(a*a - one));

          }

       }

       else

       {

          //

          // This is the Hull et al exception handling code from Fig 6 of their paper:

          //

          if(y <= (std::numeric_limits<T>::epsilon() * std::fabs(xm1)))

          {

             if(x < one)

             {

                real = std::acos(x);

                imag = y / std::sqrt(xp1*(one-x));

             }

             else

             {

                real = 0;

                if(((std::numeric_limits<T>::max)() / xp1) > xm1)

                {

                   // xp1 * xm1 won't overflow:

                   imag = boost::math::log1p(xm1 + std::sqrt(xp1*xm1));

                }

                else

                {

                   imag = log_two + std::log(x);

                }

             }

          }

          else if(y <= safe_min)

          {

             // There is an assumption in Hull et al's analysis that

             // if we get here then x == 1.  This is true for all "good"

             // machines where :

             // 

             // E^2 > 8*sqrt(u); with:

             //

             // E =  std::numeric_limits<T>::epsilon()

             // u = (std::numeric_limits<T>::min)()

             //

             // Hull et al provide alternative code for "bad" machines

             // but we have no way to test that here, so for now just assert

             // on the assumption:

             //

             BOOST_ASSERT(x == 1);

             real = std::sqrt(y);

             imag = std::sqrt(y);

          }

          else if(std::numeric_limits<T>::epsilon() * y - one >= x)

          {

             real = half_pi;

             imag = log_two + std::log(y);

          }

          else if(x > one)

          {

             real = std::atan(y/x);

             T xoy = x/y;

             imag = log_two + std::log(y) + half * boost::math::log1p(xoy*xoy);

          }

          else

          {

             real = half_pi;

             T a = std::sqrt(one + y*y);

             imag = half * boost::math::log1p(static_cast<T>(2)*y*(y+a));

          }

       }

    }

    //

    // Finish off by working out the sign of the result:

    //

    if(z.real() < 0)

       real = s_pi - real;

    if(z.imag() > 0)

       imag = -imag;

    return std::complex<T>(real, imag);

 }

 } } // namespaces

 #endif // BOOST_MATH_COMPLEX_ACOS_INCLUDED