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