math/complex/asin.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_ASIN_INCLUDED
0006 #define BOOST_MATH_COMPLEX_ASIN_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")]] inline std::complex<T> asin(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 the value of 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 #ifdef _MSC_VER
0047 #pragma warning(push)
0048 #pragma warning(disable:4127)
0049 #endif
0050    //
0051    // Get real and imaginary parts, discard the signs as we can
0052    // figure out the sign of the result later:
0053    //
0054    T x = std::fabs(z.real());
0055    T y = std::fabs(z.imag());
0056    T real, imag;  // our results
0057
0058    //
0059    // Begin by handling the special cases for infinities and nan's
0060    // specified in C99, most of this is handled by the regular logic
0061    // below, but handling it as a special case prevents overflow/underflow
0062    // arithmetic which may trip up some machines:
0063    //
0064    if((boost::math::isnan)(x))
0065    {
0066       if((boost::math::isnan)(y))
0067          return std::complex<T>(x, x);
0068       if((boost::math::isinf)(y))
0069       {
0070          real = x;
0071          imag = std::numeric_limits<T>::infinity();
0072       }
0073       else
0074          return std::complex<T>(x, x);
0075    }
0076    else if((boost::math::isnan)(y))
0077    {
0078       if(x == 0)
0079       {
0080          real = 0;
0081          imag = y;
0082       }
0083       else if((boost::math::isinf)(x))
0084       {
0085          real = y;
0086          imag = std::numeric_limits<T>::infinity();
0087       }
0088       else
0089          return std::complex<T>(y, y);
0090    }
0091    else if((boost::math::isinf)(x))
0092    {
0093       if((boost::math::isinf)(y))
0094       {
0095          real = quarter_pi;
0096          imag = std::numeric_limits<T>::infinity();
0097       }
0098       else
0099       {
0100          real = half_pi;
0101          imag = std::numeric_limits<T>::infinity();
0102       }
0103    }
0104    else if((boost::math::isinf)(y))
0105    {
0106       real = 0;
0107       imag = std::numeric_limits<T>::infinity();
0108    }
0109    else
0110    {
0111       //
0112       // special case for real numbers:
0113       //
0114       if((y == 0) && (x <= one))
0115          return std::complex<T>(std::asin(z.real()), z.imag());
0116       //
0117       // Figure out if our input is within the "safe area" identified by Hull et al.
0118       // This would be more efficient with portable floating point exception handling;
0119       // fortunately the quantities M and u identified by Hull et al (figure 3),
0120       // match with the max and min methods of numeric_limits<T>.
0121       //
0122       T safe_max = detail::safe_max(static_cast<T>(8));
0123       T safe_min = detail::safe_min(static_cast<T>(4));
0124
0125       T xp1 = one + x;
0126       T xm1 = x - one;
0127
0128       if((x < safe_max) && (x > safe_min) && (y < safe_max) && (y > safe_min))
0129       {
0130          T yy = y * y;
0131          T r = std::sqrt(xp1*xp1 + yy);
0132          T s = std::sqrt(xm1*xm1 + yy);
0133          T a = half * (r + s);
0134          T b = x / a;
0135
0136          if(b <= b_crossover)
0137          {
0138             real = std::asin(b);
0139          }
0140          else
0141          {
0142             T apx = a + x;
0143             if(x <= one)
0144             {
0145                real = std::atan(x/std::sqrt(half * apx * (yy /(r + xp1) + (s-xm1))));
0146             }
0147             else
0148             {
0149                real = std::atan(x/(y * std::sqrt(half * (apx/(r + xp1) + apx/(s+xm1)))));
0150             }
0151          }
0152
0153          if(a <= a_crossover)
0154          {
0155             T am1;
0156             if(x < one)
0157             {
0158                am1 = half * (yy/(r + xp1) + yy/(s - xm1));
0159             }
0160             else
0161             {
0162                am1 = half * (yy/(r + xp1) + (s + xm1));
0163             }
0164             imag = boost::math::log1p(am1 + std::sqrt(am1 * (a + one)));
0165          }
0166          else
0167          {
0168             imag = std::log(a + std::sqrt(a*a - one));
0169          }
0170       }
0171       else
0172       {
0173          //
0174          // This is the Hull et al exception handling code from Fig 3 of their paper:
0175          //
0176          if(y <= (std::numeric_limits<T>::epsilon() * std::fabs(xm1)))
0177          {
0178             if(x < one)
0179             {
0180                real = std::asin(x);
0181                imag = y / std::sqrt(-xp1*xm1);
0182             }
0183             else
0184             {
0185                real = half_pi;
0186                if(((std::numeric_limits<T>::max)() / xp1) > xm1)
0187                {
0188                   // xp1 * xm1 won't overflow:
0189                   imag = boost::math::log1p(xm1 + std::sqrt(xp1*xm1));
0190                }
0191                else
0192                {
0193                   imag = log_two + std::log(x);
0194                }
0195             }
0196          }
0197          else if(y <= safe_min)
0198          {
0199             // There is an assumption in Hull et al's analysis that
0200             // if we get here then x == 1.  This is true for all "good"
0201             // machines where :
0202             //
0203             // E^2 > 8*sqrt(u); with:
0204             //
0205             // E =  std::numeric_limits<T>::epsilon()
0206             // u = (std::numeric_limits<T>::min)()
0207             //
0208             // Hull et al provide alternative code for "bad" machines
0209             // but we have no way to test that here, so for now just assert
0210             // on the assumption:
0211             //
0212             BOOST_MATH_ASSERT(x == 1);
0213             real = half_pi - std::sqrt(y);
0214             imag = std::sqrt(y);
0215          }
0216          else if(std::numeric_limits<T>::epsilon() * y - one >= x)
0217          {
0218             real = x/y; // This can underflow!
0219             imag = log_two + std::log(y);
0220          }
0221          else if(x > one)
0222          {
0223             real = std::atan(x/y);
0224             T xoy = x/y;
0225             imag = log_two + std::log(y) + half * boost::math::log1p(xoy*xoy);
0226          }
0227          else
0228          {
0229             T a = std::sqrt(one + y*y);
0230             real = x/a; // This can underflow!
0231             imag = half * boost::math::log1p(static_cast<T>(2)*y*(y+a));
0232          }
0233       }
0234    }
0235
0236    //
0237    // Finish off by working out the sign of the result:
0238    //
0239    if((boost::math::signbit)(z.real()))
0240       real = (boost::math::changesign)(real);
0241    if((boost::math::signbit)(z.imag()))
0242       imag = (boost::math::changesign)(imag);
0243
0244    return std::complex<T>(real, imag);
0245 #ifdef _MSC_VER
0246 #pragma warning(pop)
0247 #endif
0248 }
0249
0250 } } // namespaces
0251
0252 #endif // BOOST_MATH_COMPLEX_ASIN_INCLUDED