ceph/src/boost/libs/math/include/boost/math/complex/asin.hpp

   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    // note that we have tweaked the value of a_crossover as per https://svn.boost.org/trac/boost/ticket/7290:
  36    //
  37    static const T one = static_cast<T>(1);
  38    //static const T two = static_cast<T>(2);
  39    static const T half = static_cast<T>(0.5L);
  40    static const T a_crossover = static_cast<T>(10);
  41    static const T b_crossover = static_cast<T>(0.6417L);
  42    static const T s_pi = boost::math::constants::pi<T>();
  43    static const T half_pi = s_pi / 2;
  44    static const T log_two = boost::math::constants::ln_two<T>();
  45    static const T quarter_pi = s_pi / 4;
  46 #ifdef BOOST_MSVC
  47 #pragma warning(push)
  48 #pragma warning(disable:4127)
  49 #endif
  50    //
  51    // Get real and imaginary parts, discard the signs as we can
  52    // figure out the sign of the result later:
  53    //
  54    T x = std::fabs(z.real());
  55    T y = std::fabs(z.imag());
  56    T real, imag;  // our results
  57
  58    //
  59    // Begin by handling the special cases for infinities and nan's
  60    // specified in C99, most of this is handled by the regular logic
  61    // below, but handling it as a special case prevents overflow/underflow
  62    // arithmetic which may trip up some machines:
  63    //
  64    if((boost::math::isnan)(x))
  65    {
  66       if((boost::math::isnan)(y))
  67          return std::complex<T>(x, x);
  68       if((boost::math::isinf)(y))
  69       {
  70          real = x;
  71          imag = std::numeric_limits<T>::infinity();
  72       }
  73       else
  74          return std::complex<T>(x, x);
  75    }
  76    else if((boost::math::isnan)(y))
  77    {
  78       if(x == 0)
  79       {
  80          real = 0;
  81          imag = y;
  82       }
  83       else if((boost::math::isinf)(x))
  84       {
  85          real = y;
  86          imag = std::numeric_limits<T>::infinity();
  87       }
  88       else
  89          return std::complex<T>(y, y);
  90    }
  91    else if((boost::math::isinf)(x))
  92    {
  93       if((boost::math::isinf)(y))
  94       {
  95          real = quarter_pi;
  96          imag = std::numeric_limits<T>::infinity();
  97       }
  98       else
  99       {
 100          real = half_pi;
 101          imag = std::numeric_limits<T>::infinity();
 102       }
 103    }
 104    else if((boost::math::isinf)(y))
 105    {
 106       real = 0;
 107       imag = std::numeric_limits<T>::infinity();
 108    }
 109    else
 110    {
 111       //
 112       // special case for real numbers:
 113       //
 114       if((y == 0) && (x <= one))
 115          return std::complex<T>(std::asin(z.real()), z.imag());
 116       //
 117       // Figure out if our input is within the "safe area" identified by Hull et al.
 118       // This would be more efficient with portable floating point exception handling;
 119       // fortunately the quantities M and u identified by Hull et al (figure 3),
 120       // match with the max and min methods of numeric_limits<T>.
 121       //
 122       T safe_max = detail::safe_max(static_cast<T>(8));
 123       T safe_min = detail::safe_min(static_cast<T>(4));
 124
 125       T xp1 = one + x;
 126       T xm1 = x - one;
 127
 128       if((x < safe_max) && (x > safe_min) && (y < safe_max) && (y > safe_min))
 129       {
 130          T yy = y * y;
 131          T r = std::sqrt(xp1*xp1 + yy);
 132          T s = std::sqrt(xm1*xm1 + yy);
 133          T a = half * (r + s);
 134          T b = x / a;
 135
 136          if(b <= b_crossover)
 137          {
 138             real = std::asin(b);
 139          }
 140          else
 141          {
 142             T apx = a + x;
 143             if(x <= one)
 144             {
 145                real = std::atan(x/std::sqrt(half * apx * (yy /(r + xp1) + (s-xm1))));
 146             }
 147             else
 148             {
 149                real = std::atan(x/(y * std::sqrt(half * (apx/(r + xp1) + apx/(s+xm1)))));
 150             }
 151          }
 152
 153          if(a <= a_crossover)
 154          {
 155             T am1;
 156             if(x < one)
 157             {
 158                am1 = half * (yy/(r + xp1) + yy/(s - xm1));
 159             }
 160             else
 161             {
 162                am1 = half * (yy/(r + xp1) + (s + xm1));
 163             }
 164             imag = boost::math::log1p(am1 + std::sqrt(am1 * (a + one)));
 165          }
 166          else
 167          {
 168             imag = std::log(a + std::sqrt(a*a - one));
 169          }
 170       }
 171       else
 172       {
 173          //
 174          // This is the Hull et al exception handling code from Fig 3 of their paper:
 175          //
 176          if(y <= (std::numeric_limits<T>::epsilon() * std::fabs(xm1)))
 177          {
 178             if(x < one)
 179             {
 180                real = std::asin(x);
 181                imag = y / std::sqrt(-xp1*xm1);
 182             }
 183             else
 184             {
 185                real = half_pi;
 186                if(((std::numeric_limits<T>::max)() / xp1) > xm1)
 187                {
 188                   // xp1 * xm1 won't overflow:
 189                   imag = boost::math::log1p(xm1 + std::sqrt(xp1*xm1));
 190                }
 191                else
 192                {
 193                   imag = log_two + std::log(x);
 194                }
 195             }
 196          }
 197          else if(y <= safe_min)
 198          {
 199             // There is an assumption in Hull et al's analysis that
 200             // if we get here then x == 1.  This is true for all "good"
 201             // machines where :
 202             //
 203             // E^2 > 8*sqrt(u); with:
 204             //
 205             // E =  std::numeric_limits<T>::epsilon()
 206             // u = (std::numeric_limits<T>::min)()
 207             //
 208             // Hull et al provide alternative code for "bad" machines
 209             // but we have no way to test that here, so for now just assert
 210             // on the assumption:
 211             //
 212             BOOST_ASSERT(x == 1);
 213             real = half_pi - std::sqrt(y);
 214             imag = std::sqrt(y);
 215          }
 216          else if(std::numeric_limits<T>::epsilon() * y - one >= x)
 217          {
 218             real = x/y; // This can underflow!
 219             imag = log_two + std::log(y);
 220          }
 221          else if(x > one)
 222          {
 223             real = std::atan(x/y);
 224             T xoy = x/y;
 225             imag = log_two + std::log(y) + half * boost::math::log1p(xoy*xoy);
 226          }
 227          else
 228          {
 229             T a = std::sqrt(one + y*y);
 230             real = x/a; // This can underflow!
 231             imag = half * boost::math::log1p(static_cast<T>(2)*y*(y+a));
 232          }
 233       }
 234    }
 235
 236    //
 237    // Finish off by working out the sign of the result:
 238    //
 239    if((boost::math::signbit)(z.real()))
 240       real = (boost::math::changesign)(real);
 241    if((boost::math::signbit)(z.imag()))
 242       imag = (boost::math::changesign)(imag);
 243
 244    return std::complex<T>(real, imag);
 245 #ifdef BOOST_MSVC
 246 #pragma warning(pop)
 247 #endif
 248 }
 249
 250 } } // namespaces
 251
 252 #endif // BOOST_MATH_COMPLEX_ASIN_INCLUDED