math/complex/acos.hpp

0001 //  (C) Copyright John Maddock 2005.
0002 //  Distributed under the Boost Software License, Version 1.0. (See accompanying
0003 //  file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
0004
0005 #ifndef BOOST_MATH_COMPLEX_ACOS_INCLUDED
0006 #define BOOST_MATH_COMPLEX_ACOS_INCLUDED
0007
0008 #ifndef BOOST_MATH_COMPLEX_DETAILS_INCLUDED
0009 #  include <boost/math/complex/details.hpp>
0010 #endif
0011 #ifndef BOOST_MATH_LOG1P_INCLUDED
0012 #  include <boost/math/special_functions/log1p.hpp>
0013 #endif
0014 #include <boost/math/tools/assert.hpp>
0015
0016 #ifdef BOOST_NO_STDC_NAMESPACE
0017 namespace std{ using ::sqrt; using ::fabs; using ::acos; using ::asin; using ::atan; using ::atan2; }
0018 #endif
0019
0020 namespace boost{ namespace math{
0021
0022 template<class T>
0023 [[deprecated("Replaced by C++11")]] std::complex<T> acos(const std::complex<T>& z)
0024 {
0025    //
0026    // This implementation is a transcription of the pseudo-code in:
0027    //
0028    // "Implementing the Complex Arcsine and Arccosine Functions using Exception Handling."
0029    // T E Hull, Thomas F Fairgrieve and Ping Tak Peter Tang.
0030    // ACM Transactions on Mathematical Software, Vol 23, No 3, Sept 1997.
0031    //
0032
0033    //
0034    // These static constants should really be in a maths constants library,
0035    // note that we have tweaked a_crossover as per: https://svn.boost.org/trac/boost/ticket/7290
0036    //
0037    static const T one = static_cast<T>(1);
0038    //static const T two = static_cast<T>(2);
0039    static const T half = static_cast<T>(0.5L);
0040    static const T a_crossover = static_cast<T>(10);
0041    static const T b_crossover = static_cast<T>(0.6417L);
0042    static const T s_pi = boost::math::constants::pi<T>();
0043    static const T half_pi = s_pi / 2;
0044    static const T log_two = boost::math::constants::ln_two<T>();
0045    static const T quarter_pi = s_pi / 4;
0046
0047 #ifdef _MSC_VER
0048 #pragma warning(push)
0049 #pragma warning(disable:4127)
0050 #endif
0051    //
0052    // Get real and imaginary parts, discard the signs as we can
0053    // figure out the sign of the result later:
0054    //
0055    T x = std::fabs(z.real());
0056    T y = std::fabs(z.imag());
0057
0058    T real, imag; // these hold our result
0059
0060    //
0061    // Handle special cases specified by the C99 standard,
0062    // many of these special cases aren't really needed here,
0063    // but doing it this way prevents overflow/underflow arithmetic
0064    // in the main body of the logic, which may trip up some machines:
0065    //
0066    if((boost::math::isinf)(x))
0067    {
0068       if((boost::math::isinf)(y))
0069       {
0070          real = quarter_pi;
0071          imag = std::numeric_limits<T>::infinity();
0072       }
0073       else if((boost::math::isnan)(y))
0074       {
0075          return std::complex<T>(y, -std::numeric_limits<T>::infinity());
0076       }
0077       else
0078       {
0079          // y is not infinity or nan:
0080          real = 0;
0081          imag = std::numeric_limits<T>::infinity();
0082       }
0083    }
0084    else if((boost::math::isnan)(x))
0085    {
0086       if((boost::math::isinf)(y))
0087          return std::complex<T>(x, ((boost::math::signbit)(z.imag())) ? std::numeric_limits<T>::infinity() :  -std::numeric_limits<T>::infinity());
0088       return std::complex<T>(x, x);
0089    }
0090    else if((boost::math::isinf)(y))
0091    {
0092       real = half_pi;
0093       imag = std::numeric_limits<T>::infinity();
0094    }
0095    else if((boost::math::isnan)(y))
0096    {
0097       return std::complex<T>((x == 0) ? half_pi : y, y);
0098    }
0099    else
0100    {
0101       //
0102       // What follows is the regular Hull et al code,
0103       // begin with the special case for real numbers:
0104       //
0105       if((y == 0) && (x <= one))
0106          return std::complex<T>((x == 0) ? half_pi : std::acos(z.real()), (boost::math::changesign)(z.imag()));
0107       //
0108       // Figure out if our input is within the "safe area" identified by Hull et al.
0109       // This would be more efficient with portable floating point exception handling;
0110       // fortunately the quantities M and u identified by Hull et al (figure 3),
0111       // match with the max and min methods of numeric_limits<T>.
0112       //
0113       T safe_max = detail::safe_max(static_cast<T>(8));
0114       T safe_min = detail::safe_min(static_cast<T>(4));
0115
0116       T xp1 = one + x;
0117       T xm1 = x - one;
0118
0119       if((x < safe_max) && (x > safe_min) && (y < safe_max) && (y > safe_min))
0120       {
0121          T yy = y * y;
0122          T r = std::sqrt(xp1*xp1 + yy);
0123          T s = std::sqrt(xm1*xm1 + yy);
0124          T a = half * (r + s);
0125          T b = x / a;
0126
0127          if(b <= b_crossover)
0128          {
0129             real = std::acos(b);
0130          }
0131          else
0132          {
0133             T apx = a + x;
0134             if(x <= one)
0135             {
0136                real = std::atan(std::sqrt(half * apx * (yy /(r + xp1) + (s-xm1)))/x);
0137             }
0138             else
0139             {
0140                real = std::atan((y * std::sqrt(half * (apx/(r + xp1) + apx/(s+xm1))))/x);
0141             }
0142          }
0143
0144          if(a <= a_crossover)
0145          {
0146             T am1;
0147             if(x < one)
0148             {
0149                am1 = half * (yy/(r + xp1) + yy/(s - xm1));
0150             }
0151             else
0152             {
0153                am1 = half * (yy/(r + xp1) + (s + xm1));
0154             }
0155             imag = boost::math::log1p(am1 + std::sqrt(am1 * (a + one)));
0156          }
0157          else
0158          {
0159             imag = std::log(a + std::sqrt(a*a - one));
0160          }
0161       }
0162       else
0163       {
0164          //
0165          // This is the Hull et al exception handling code from Fig 6 of their paper:
0166          //
0167          if(y <= (std::numeric_limits<T>::epsilon() * std::fabs(xm1)))
0168          {
0169             if(x < one)
0170             {
0171                real = std::acos(x);
0172                imag = y / std::sqrt(xp1*(one-x));
0173             }
0174             else
0175             {
0176                // This deviates from Hull et al's paper as per https://svn.boost.org/trac/boost/ticket/7290
0177                if(((std::numeric_limits<T>::max)() / xp1) > xm1)
0178                {
0179                   // xp1 * xm1 won't overflow:
0180                   real = y / std::sqrt(xm1*xp1);
0181                   imag = boost::math::log1p(xm1 + std::sqrt(xp1*xm1));
0182                }
0183                else
0184                {
0185                   real = y / x;
0186                   imag = log_two + std::log(x);
0187                }
0188             }
0189          }
0190          else if(y <= safe_min)
0191          {
0192             // There is an assumption in Hull et al's analysis that
0193             // if we get here then x == 1.  This is true for all "good"
0194             // machines where :
0195             //
0196             // E^2 > 8*sqrt(u); with:
0197             //
0198             // E =  std::numeric_limits<T>::epsilon()
0199             // u = (std::numeric_limits<T>::min)()
0200             //
0201             // Hull et al provide alternative code for "bad" machines
0202             // but we have no way to test that here, so for now just assert
0203             // on the assumption:
0204             //
0205             BOOST_MATH_ASSERT(x == 1);
0206             real = std::sqrt(y);
0207             imag = std::sqrt(y);
0208          }
0209          else if(std::numeric_limits<T>::epsilon() * y - one >= x)
0210          {
0211             real = half_pi;
0212             imag = log_two + std::log(y);
0213          }
0214          else if(x > one)
0215          {
0216             real = std::atan(y/x);
0217             T xoy = x/y;
0218             imag = log_two + std::log(y) + half * boost::math::log1p(xoy*xoy);
0219          }
0220          else
0221          {
0222             real = half_pi;
0223             T a = std::sqrt(one + y*y);
0224             imag = half * boost::math::log1p(static_cast<T>(2)*y*(y+a));
0225          }
0226       }
0227    }
0228
0229    //
0230    // Finish off by working out the sign of the result:
0231    //
0232    if((boost::math::signbit)(z.real()))
0233       real = s_pi - real;
0234    if(!(boost::math::signbit)(z.imag()))
0235       imag = (boost::math::changesign)(imag);
0236
0237    return std::complex<T>(real, imag);
0238 #ifdef _MSC_VER
0239 #pragma warning(pop)
0240 #endif
0241 }
0242
0243 } } // namespaces
0244
0245 #endif // BOOST_MATH_COMPLEX_ACOS_INCLUDED