1 //  (C) Copyright John Maddock 2005.

     2 //  Distributed under the Boost Software License, Version 1.0. (See accompanying

     3 //  file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)

     4 

     5 #ifndef BOOST_MATH_COMPLEX_ACOS_INCLUDED

     6 #define BOOST_MATH_COMPLEX_ACOS_INCLUDED

     7 

     8 #ifndef BOOST_MATH_COMPLEX_DETAILS_INCLUDED

     9 #  include <boost/math/complex/details.hpp>

    10 #endif

    11 #ifndef BOOST_MATH_LOG1P_INCLUDED

    12 #  include <boost/math/special_functions/log1p.hpp>

    13 #endif

    14 #include <boost/assert.hpp>

    15 

    16 #ifdef BOOST_NO_STDC_NAMESPACE

    17 namespace std{ using ::sqrt; using ::fabs; using ::acos; using ::asin; using ::atan; using ::atan2; }

    18 #endif

    19 

    20 namespace boost{ namespace math{

    21 

    22 template<class T> 

    23 std::complex<T> acos(const std::complex<T>& z)

    24 {

    25    //

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

    27    //

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

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

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

    31    //

    32 

    33    //

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

    35    //

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

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

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

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

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

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

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

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

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

    45    

    46    //

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

    48    // figure out the sign of the result later:

    49    //

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

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

    52 

    53    T real, imag; // these hold our result

    54 

    55    // 

    56    // Handle special cases specified by the C99 standard,

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

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

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

    60    //

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

    62    {

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

    64       {

    65          real = quarter_pi;

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

    67       }

    68       else if(detail::test_is_nan(y))

    69       {

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

    71       }

    72       else

    73       {

    74          // y is not infinity or nan:

    75          real = 0;

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

    77       }

    78    }

    79    else if(detail::test_is_nan(x))

    80    {

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

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

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

    84    }

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

    86    {

    87       real = half_pi;

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

    89    }

    90    else if(detail::test_is_nan(y))

    91    {

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

    93    }

    94    else

    95    {

    96       //

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

    98       // begin with the special case for real numbers:

    99       //

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

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

   102       //

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

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

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

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

   107       //

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

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

   110 

   111       T xp1 = one + x;

   112       T xm1 = x - one;

   113 

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

   115       {

   116          T yy = y * y;

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

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

   119          T a = half * (r + s);

   120          T b = x / a;

   121 

   122          if(b <= b_crossover)

   123          {

   124             real = std::acos(b);

   125          }

   126          else

   127          {

   128             T apx = a + x;

   129             if(x <= one)

   130             {

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

   132             }

   133             else

   134             {

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

   136             }

   137          }

   138 

   139          if(a <= a_crossover)

   140          {

   141             T am1;

   142             if(x < one)

   143             {

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

   145             }

   146             else

   147             {

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

   149             }

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

   151          }

   152          else

   153          {

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

   155          }

   156       }

   157       else

   158       {

   159          //

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

   161          //

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

   163          {

   164             if(x < one)

   165             {

   166                real = std::acos(x);

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

   168             }

   169             else

   170             {

   171                real = 0;

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

   173                {

   174                   // xp1 * xm1 won't overflow:

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

   176                }

   177                else

   178                {

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

   180                }

   181             }

   182          }

   183          else if(y <= safe_min)

   184          {

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

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

   187             // machines where :

   188             // 

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

   190             //

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

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

   193             //

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

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

   196             // on the assumption:

   197             //

   198             BOOST_ASSERT(x == 1);

   199             real = std::sqrt(y);

   200             imag = std::sqrt(y);

   201          }

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

   203          {

   204             real = half_pi;

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

   206          }

   207          else if(x > one)

   208          {

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

   210             T xoy = x/y;

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

   212          }

   213          else

   214          {

   215             real = half_pi;

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

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

   218          }

   219       }

   220    }

   221 

   222    //

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

   224    //

   225    if(z.real() < 0)

   226       real = s_pi - real;

   227    if(z.imag() > 0)

   228       imag = -imag;

   229 

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

   231 }

   232 

   233 } } // namespaces

   234 

   235 #endif // BOOST_MATH_COMPLEX_ACOS_INCLUDED