epoc32/include/stdapis/boost/math/complex/asin.hpp
author William Roberts <williamr@symbian.org>
Wed, 31 Mar 2010 12:27:01 +0100
branchSymbian2
changeset 3 e1b950c65cb4
permissions -rw-r--r--
Attempt to represent the S^2->S^3 header reorganisation as a series of "hg rename" operations
     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