253 lines · plain
1// (C) Copyright John Maddock 2005.2// Distributed under the Boost Software License, Version 1.0. (See accompanying3// file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)4 5#ifndef BOOST_MATH_COMPLEX_ASIN_INCLUDED6#define BOOST_MATH_COMPLEX_ASIN_INCLUDED7 8#ifndef BOOST_MATH_COMPLEX_DETAILS_INCLUDED9# include <boost/math/complex/details.hpp>10#endif11#ifndef BOOST_MATH_LOG1P_INCLUDED12# include <boost/math/special_functions/log1p.hpp>13#endif14#include <boost/math/tools/assert.hpp>15 16#ifdef BOOST_NO_STDC_NAMESPACE17namespace std{ using ::sqrt; using ::fabs; using ::acos; using ::asin; using ::atan; using ::atan2; }18#endif19 20namespace boost{ namespace math{21 22template<class T> 23[[deprecated("Replaced by C++11")]] 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 _MSC_VER47#pragma warning(push)48#pragma warning(disable:4127)49#endif50 //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 results57 58 //59 // Begin by handling the special cases for infinities and nan's60 // specified in C99, most of this is handled by the regular logic61 // below, but handling it as a special case prevents overflow/underflow62 // 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 else74 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 else89 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 else99 {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 else110 {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 else141 {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 else148 {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 else161 {162 am1 = half * (yy/(r + xp1) + (s + xm1));163 }164 imag = boost::math::log1p(am1 + std::sqrt(am1 * (a + one)));165 }166 else167 {168 imag = std::log(a + std::sqrt(a*a - one));169 }170 }171 else172 {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 else184 {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 else192 {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 that200 // 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" machines209 // but we have no way to test that here, so for now just assert210 // on the assumption:211 //212 BOOST_MATH_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 else228 {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 _MSC_VER246#pragma warning(pop)247#endif248}249 250} } // namespaces251 252#endif // BOOST_MATH_COMPLEX_ASIN_INCLUDED253