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
williamr@2
     1
//  (C) Copyright John Maddock 2005.
williamr@2
     2
//  Distributed under the Boost Software License, Version 1.0. (See accompanying
williamr@2
     3
//  file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
williamr@2
     4
williamr@2
     5
#ifndef BOOST_MATH_COMPLEX_ASIN_INCLUDED
williamr@2
     6
#define BOOST_MATH_COMPLEX_ASIN_INCLUDED
williamr@2
     7
williamr@2
     8
#ifndef BOOST_MATH_COMPLEX_DETAILS_INCLUDED
williamr@2
     9
#  include <boost/math/complex/details.hpp>
williamr@2
    10
#endif
williamr@2
    11
#ifndef BOOST_MATH_LOG1P_INCLUDED
williamr@2
    12
#  include <boost/math/special_functions/log1p.hpp>
williamr@2
    13
#endif
williamr@2
    14
#include <boost/assert.hpp>
williamr@2
    15
williamr@2
    16
#ifdef BOOST_NO_STDC_NAMESPACE
williamr@2
    17
namespace std{ using ::sqrt; using ::fabs; using ::acos; using ::asin; using ::atan; using ::atan2; }
williamr@2
    18
#endif
williamr@2
    19
williamr@2
    20
namespace boost{ namespace math{
williamr@2
    21
williamr@2
    22
template<class T> 
williamr@2
    23
inline std::complex<T> asin(const std::complex<T>& z)
williamr@2
    24
{
williamr@2
    25
   //
williamr@2
    26
   // This implementation is a transcription of the pseudo-code in:
williamr@2
    27
   //
williamr@2
    28
   // "Implementing the complex Arcsine and Arccosine Functions using Exception Handling."
williamr@2
    29
   // T E Hull, Thomas F Fairgrieve and Ping Tak Peter Tang.
williamr@2
    30
   // ACM Transactions on Mathematical Software, Vol 23, No 3, Sept 1997.
williamr@2
    31
   //
williamr@2
    32
williamr@2
    33
   //
williamr@2
    34
   // These static constants should really be in a maths constants library:
williamr@2
    35
   //
williamr@2
    36
   static const T one = static_cast<T>(1);
williamr@2
    37
   //static const T two = static_cast<T>(2);
williamr@2
    38
   static const T half = static_cast<T>(0.5L);
williamr@2
    39
   static const T a_crossover = static_cast<T>(1.5L);
williamr@2
    40
   static const T b_crossover = static_cast<T>(0.6417L);
williamr@2
    41
   //static const T pi = static_cast<T>(3.141592653589793238462643383279502884197L);
williamr@2
    42
   static const T half_pi = static_cast<T>(1.57079632679489661923132169163975144L);
williamr@2
    43
   static const T log_two = static_cast<T>(0.69314718055994530941723212145817657L);
williamr@2
    44
   static const T quarter_pi = static_cast<T>(0.78539816339744830961566084581987572L);
williamr@2
    45
   
williamr@2
    46
   //
williamr@2
    47
   // Get real and imaginary parts, discard the signs as we can 
williamr@2
    48
   // figure out the sign of the result later:
williamr@2
    49
   //
williamr@2
    50
   T x = std::fabs(z.real());
williamr@2
    51
   T y = std::fabs(z.imag());
williamr@2
    52
   T real, imag;  // our results
williamr@2
    53
williamr@2
    54
   //
williamr@2
    55
   // Begin by handling the special cases for infinities and nan's
williamr@2
    56
   // specified in C99, most of this is handled by the regular logic
williamr@2
    57
   // below, but handling it as a special case prevents overflow/underflow
williamr@2
    58
   // arithmetic which may trip up some machines:
williamr@2
    59
   //
williamr@2
    60
   if(detail::test_is_nan(x))
williamr@2
    61
   {
williamr@2
    62
      if(detail::test_is_nan(y))
williamr@2
    63
         return std::complex<T>(x, x);
williamr@2
    64
      if(std::numeric_limits<T>::has_infinity && (y == std::numeric_limits<T>::infinity()))
williamr@2
    65
      {
williamr@2
    66
         real = x;
williamr@2
    67
         imag = std::numeric_limits<T>::infinity();
williamr@2
    68
      }
williamr@2
    69
      else
williamr@2
    70
         return std::complex<T>(x, x);
williamr@2
    71
   }
williamr@2
    72
   else if(detail::test_is_nan(y))
williamr@2
    73
   {
williamr@2
    74
      if(x == 0)
williamr@2
    75
      {
williamr@2
    76
         real = 0;
williamr@2
    77
         imag = y;
williamr@2
    78
      }
williamr@2
    79
      else if(std::numeric_limits<T>::has_infinity && (x == std::numeric_limits<T>::infinity()))
williamr@2
    80
      {
williamr@2
    81
         real = y;
williamr@2
    82
         imag = std::numeric_limits<T>::infinity();
williamr@2
    83
      }
williamr@2
    84
      else
williamr@2
    85
         return std::complex<T>(y, y);
williamr@2
    86
   }
williamr@2
    87
   else if(std::numeric_limits<T>::has_infinity && (x == std::numeric_limits<T>::infinity()))
williamr@2
    88
   {
williamr@2
    89
      if(y == std::numeric_limits<T>::infinity())
williamr@2
    90
      {
williamr@2
    91
         real = quarter_pi;
williamr@2
    92
         imag = std::numeric_limits<T>::infinity();
williamr@2
    93
      }
williamr@2
    94
      else
williamr@2
    95
      {
williamr@2
    96
         real = half_pi;
williamr@2
    97
         imag = std::numeric_limits<T>::infinity();
williamr@2
    98
      }
williamr@2
    99
   }
williamr@2
   100
   else if(std::numeric_limits<T>::has_infinity && (y == std::numeric_limits<T>::infinity()))
williamr@2
   101
   {
williamr@2
   102
      real = 0;
williamr@2
   103
      imag = std::numeric_limits<T>::infinity();
williamr@2
   104
   }
williamr@2
   105
   else
williamr@2
   106
   {
williamr@2
   107
      //
williamr@2
   108
      // special case for real numbers:
williamr@2
   109
      //
williamr@2
   110
      if((y == 0) && (x <= one))
williamr@2
   111
         return std::complex<T>(std::asin(z.real()));
williamr@2
   112
      //
williamr@2
   113
      // Figure out if our input is within the "safe area" identified by Hull et al.
williamr@2
   114
      // This would be more efficient with portable floating point exception handling;
williamr@2
   115
      // fortunately the quantities M and u identified by Hull et al (figure 3), 
williamr@2
   116
      // match with the max and min methods of numeric_limits<T>.
williamr@2
   117
      //
williamr@2
   118
      T safe_max = detail::safe_max(static_cast<T>(8));
williamr@2
   119
      T safe_min = detail::safe_min(static_cast<T>(4));
williamr@2
   120
williamr@2
   121
      T xp1 = one + x;
williamr@2
   122
      T xm1 = x - one;
williamr@2
   123
williamr@2
   124
      if((x < safe_max) && (x > safe_min) && (y < safe_max) && (y > safe_min))
williamr@2
   125
      {
williamr@2
   126
         T yy = y * y;
williamr@2
   127
         T r = std::sqrt(xp1*xp1 + yy);
williamr@2
   128
         T s = std::sqrt(xm1*xm1 + yy);
williamr@2
   129
         T a = half * (r + s);
williamr@2
   130
         T b = x / a;
williamr@2
   131
williamr@2
   132
         if(b <= b_crossover)
williamr@2
   133
         {
williamr@2
   134
            real = std::asin(b);
williamr@2
   135
         }
williamr@2
   136
         else
williamr@2
   137
         {
williamr@2
   138
            T apx = a + x;
williamr@2
   139
            if(x <= one)
williamr@2
   140
            {
williamr@2
   141
               real = std::atan(x/std::sqrt(half * apx * (yy /(r + xp1) + (s-xm1))));
williamr@2
   142
            }
williamr@2
   143
            else
williamr@2
   144
            {
williamr@2
   145
               real = std::atan(x/(y * std::sqrt(half * (apx/(r + xp1) + apx/(s+xm1)))));
williamr@2
   146
            }
williamr@2
   147
         }
williamr@2
   148
williamr@2
   149
         if(a <= a_crossover)
williamr@2
   150
         {
williamr@2
   151
            T am1;
williamr@2
   152
            if(x < one)
williamr@2
   153
            {
williamr@2
   154
               am1 = half * (yy/(r + xp1) + yy/(s - xm1));
williamr@2
   155
            }
williamr@2
   156
            else
williamr@2
   157
            {
williamr@2
   158
               am1 = half * (yy/(r + xp1) + (s + xm1));
williamr@2
   159
            }
williamr@2
   160
            imag = boost::math::log1p(am1 + std::sqrt(am1 * (a + one)));
williamr@2
   161
         }
williamr@2
   162
         else
williamr@2
   163
         {
williamr@2
   164
            imag = std::log(a + std::sqrt(a*a - one));
williamr@2
   165
         }
williamr@2
   166
      }
williamr@2
   167
      else
williamr@2
   168
      {
williamr@2
   169
         //
williamr@2
   170
         // This is the Hull et al exception handling code from Fig 3 of their paper:
williamr@2
   171
         //
williamr@2
   172
         if(y <= (std::numeric_limits<T>::epsilon() * std::fabs(xm1)))
williamr@2
   173
         {
williamr@2
   174
            if(x < one)
williamr@2
   175
            {
williamr@2
   176
               real = std::asin(x);
williamr@2
   177
               imag = y / std::sqrt(xp1*xm1);
williamr@2
   178
            }
williamr@2
   179
            else
williamr@2
   180
            {
williamr@2
   181
               real = half_pi;
williamr@2
   182
               if(((std::numeric_limits<T>::max)() / xp1) > xm1)
williamr@2
   183
               {
williamr@2
   184
                  // xp1 * xm1 won't overflow:
williamr@2
   185
                  imag = boost::math::log1p(xm1 + std::sqrt(xp1*xm1));
williamr@2
   186
               }
williamr@2
   187
               else
williamr@2
   188
               {
williamr@2
   189
                  imag = log_two + std::log(x);
williamr@2
   190
               }
williamr@2
   191
            }
williamr@2
   192
         }
williamr@2
   193
         else if(y <= safe_min)
williamr@2
   194
         {
williamr@2
   195
            // There is an assumption in Hull et al's analysis that
williamr@2
   196
            // if we get here then x == 1.  This is true for all "good"
williamr@2
   197
            // machines where :
williamr@2
   198
            // 
williamr@2
   199
            // E^2 > 8*sqrt(u); with:
williamr@2
   200
            //
williamr@2
   201
            // E =  std::numeric_limits<T>::epsilon()
williamr@2
   202
            // u = (std::numeric_limits<T>::min)()
williamr@2
   203
            //
williamr@2
   204
            // Hull et al provide alternative code for "bad" machines
williamr@2
   205
            // but we have no way to test that here, so for now just assert
williamr@2
   206
            // on the assumption:
williamr@2
   207
            //
williamr@2
   208
            BOOST_ASSERT(x == 1);
williamr@2
   209
            real = half_pi - std::sqrt(y);
williamr@2
   210
            imag = std::sqrt(y);
williamr@2
   211
         }
williamr@2
   212
         else if(std::numeric_limits<T>::epsilon() * y - one >= x)
williamr@2
   213
         {
williamr@2
   214
            real = x/y; // This can underflow!
williamr@2
   215
            imag = log_two + std::log(y);
williamr@2
   216
         }
williamr@2
   217
         else if(x > one)
williamr@2
   218
         {
williamr@2
   219
            real = std::atan(x/y);
williamr@2
   220
            T xoy = x/y;
williamr@2
   221
            imag = log_two + std::log(y) + half * boost::math::log1p(xoy*xoy);
williamr@2
   222
         }
williamr@2
   223
         else
williamr@2
   224
         {
williamr@2
   225
            T a = std::sqrt(one + y*y);
williamr@2
   226
            real = x/a; // This can underflow!
williamr@2
   227
            imag = half * boost::math::log1p(static_cast<T>(2)*y*(y+a));
williamr@2
   228
         }
williamr@2
   229
      }
williamr@2
   230
   }
williamr@2
   231
williamr@2
   232
   //
williamr@2
   233
   // Finish off by working out the sign of the result:
williamr@2
   234
   //
williamr@2
   235
   if(z.real() < 0)
williamr@2
   236
      real = -real;
williamr@2
   237
   if(z.imag() < 0)
williamr@2
   238
      imag = -imag;
williamr@2
   239
williamr@2
   240
   return std::complex<T>(real, imag);
williamr@2
   241
}
williamr@2
   242
williamr@2
   243
} } // namespaces
williamr@2
   244
williamr@2
   245
#endif // BOOST_MATH_COMPLEX_ASIN_INCLUDED