epoc32/include/stdapis/boost/math/complex/acos.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_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