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_ASIN_INCLUDED

     6 #define BOOST_MATH_COMPLEX_ASIN_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 inline std::complex<T> asin(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 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    T real, imag;  // our results

    53 

    54    //

    55    // Begin by handling the special cases for infinities and nan's

    56    // specified in C99, most of this is handled by the regular logic

    57    // below, but handling it as a special case prevents overflow/underflow

    58    // arithmetic which may trip up some machines:

    59    //

    60    if(detail::test_is_nan(x))

    61    {

    62       if(detail::test_is_nan(y))

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

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

    65       {

    66          real = x;

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

    68       }

    69       else

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

    71    }

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

    73    {

    74       if(x == 0)

    75       {

    76          real = 0;

    77          imag = y;

    78       }

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

    80       {

    81          real = y;

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

    83       }

    84       else

    85          return std::complex<T>(y, y);

    86    }

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

    88    {

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

    90       {

    91          real = quarter_pi;

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

    93       }

    94       else

    95       {

    96          real = half_pi;

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

    98       }

    99    }

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

   101    {

   102       real = 0;

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

   104    }

   105    else

   106    {

   107       //

   108       // special case for real numbers:

   109       //

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

   111          return std::complex<T>(std::asin(z.real()));

   112       //

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

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

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

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

   117       //

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

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

   120 

   121       T xp1 = one + x;

   122       T xm1 = x - one;

   123 

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

   125       {

   126          T yy = y * y;

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

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

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

   130          T b = x / a;

   131 

   132          if(b <= b_crossover)

   133          {

   134             real = std::asin(b);

   135          }

   136          else

   137          {

   138             T apx = a + x;

   139             if(x <= one)

   140             {

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

   142             }

   143             else

   144             {

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

   146             }

   147          }

   148 

   149          if(a <= a_crossover)

   150          {

   151             T am1;

   152             if(x < one)

   153             {

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

   155             }

   156             else

   157             {

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

   159             }

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

   161          }

   162          else

   163          {

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

   165          }

   166       }

   167       else

   168       {

   169          //

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

   171          //

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

   173          {

   174             if(x < one)

   175             {

   176                real = std::asin(x);

   177                imag = y / std::sqrt(xp1*xm1);

   178             }

   179             else

   180             {

   181                real = half_pi;

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

   183                {

   184                   // xp1 * xm1 won't overflow:

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

   186                }

   187                else

   188                {

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

   190                }

   191             }

   192          }

   193          else if(y <= safe_min)

   194          {

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

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

   197             // machines where :

   198             // 

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

   200             //

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

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

   203             //

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

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

   206             // on the assumption:

   207             //

   208             BOOST_ASSERT(x == 1);

   209             real = half_pi - std::sqrt(y);

   210             imag = std::sqrt(y);

   211          }

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

   213          {

   214             real = x/y; // This can underflow!

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

   216          }

   217          else if(x > one)

   218          {

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

   220             T xoy = x/y;

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

   222          }

   223          else

   224          {

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

   226             real = x/a; // This can underflow!

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

   228          }

   229       }

   230    }

   231 

   232    //

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

   234    //

   235    if(z.real() < 0)

   236       real = -real;

   237    if(z.imag() < 0)

   238       imag = -imag;

   239 

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

   241 }

   242 

   243 } } // namespaces

   244 

   245 #endif // BOOST_MATH_COMPLEX_ASIN_INCLUDED