246 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_ACOS_INCLUDED6#define BOOST_MATH_COMPLEX_ACOS_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")]] 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/729036 //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 _MSC_VER48#pragma warning(push)49#pragma warning(disable:4127)50#endif51 //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 result59 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 arithmetic64 // 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 else78 {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 else100 {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 else132 {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 else139 {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 else152 {153 am1 = half * (yy/(r + xp1) + (s + xm1));154 }155 imag = boost::math::log1p(am1 + std::sqrt(am1 * (a + one)));156 }157 else158 {159 imag = std::log(a + std::sqrt(a*a - one));160 }161 }162 else163 {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 else175 {176 // This deviates from Hull et al's paper as per https://svn.boost.org/trac/boost/ticket/7290177 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 else184 {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 that193 // 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" machines202 // but we have no way to test that here, so for now just assert203 // on the assumption:204 //205 BOOST_MATH_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 else221 {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 _MSC_VER239#pragma warning(pop)240#endif241}242 243} } // namespaces244 245#endif // BOOST_MATH_COMPLEX_ACOS_INCLUDED246