williamr@2: // (C) Copyright John Maddock 2005. williamr@2: // Distributed under the Boost Software License, Version 1.0. (See accompanying williamr@2: // file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) williamr@2: williamr@2: #ifndef BOOST_MATH_COMPLEX_ACOS_INCLUDED williamr@2: #define BOOST_MATH_COMPLEX_ACOS_INCLUDED williamr@2: williamr@2: #ifndef BOOST_MATH_COMPLEX_DETAILS_INCLUDED williamr@2: # include williamr@2: #endif williamr@2: #ifndef BOOST_MATH_LOG1P_INCLUDED williamr@2: # include williamr@2: #endif williamr@2: #include williamr@2: williamr@2: #ifdef BOOST_NO_STDC_NAMESPACE williamr@2: namespace std{ using ::sqrt; using ::fabs; using ::acos; using ::asin; using ::atan; using ::atan2; } williamr@2: #endif williamr@2: williamr@2: namespace boost{ namespace math{ williamr@2: williamr@2: template williamr@2: std::complex acos(const std::complex& z) williamr@2: { williamr@2: // williamr@2: // This implementation is a transcription of the pseudo-code in: williamr@2: // williamr@2: // "Implementing the Complex Arcsine and Arccosine Functions using Exception Handling." williamr@2: // T E Hull, Thomas F Fairgrieve and Ping Tak Peter Tang. williamr@2: // ACM Transactions on Mathematical Software, Vol 23, No 3, Sept 1997. williamr@2: // williamr@2: williamr@2: // williamr@2: // These static constants should really be in a maths constants library: williamr@2: // williamr@2: static const T one = static_cast(1); williamr@2: //static const T two = static_cast(2); williamr@2: static const T half = static_cast(0.5L); williamr@2: static const T a_crossover = static_cast(1.5L); williamr@2: static const T b_crossover = static_cast(0.6417L); williamr@2: static const T s_pi = static_cast(3.141592653589793238462643383279502884197L); williamr@2: static const T half_pi = static_cast(1.57079632679489661923132169163975144L); williamr@2: static const T log_two = static_cast(0.69314718055994530941723212145817657L); williamr@2: static const T quarter_pi = static_cast(0.78539816339744830961566084581987572L); williamr@2: williamr@2: // williamr@2: // Get real and imaginary parts, discard the signs as we can williamr@2: // figure out the sign of the result later: williamr@2: // williamr@2: T x = std::fabs(z.real()); williamr@2: T y = std::fabs(z.imag()); williamr@2: williamr@2: T real, imag; // these hold our result williamr@2: williamr@2: // williamr@2: // Handle special cases specified by the C99 standard, williamr@2: // many of these special cases aren't really needed here, williamr@2: // but doing it this way prevents overflow/underflow arithmetic williamr@2: // in the main body of the logic, which may trip up some machines: williamr@2: // williamr@2: if(std::numeric_limits::has_infinity && (x == std::numeric_limits::infinity())) williamr@2: { williamr@2: if(y == std::numeric_limits::infinity()) williamr@2: { williamr@2: real = quarter_pi; williamr@2: imag = std::numeric_limits::infinity(); williamr@2: } williamr@2: else if(detail::test_is_nan(y)) williamr@2: { williamr@2: return std::complex(y, -std::numeric_limits::infinity()); williamr@2: } williamr@2: else williamr@2: { williamr@2: // y is not infinity or nan: williamr@2: real = 0; williamr@2: imag = std::numeric_limits::infinity(); williamr@2: } williamr@2: } williamr@2: else if(detail::test_is_nan(x)) williamr@2: { williamr@2: if(y == std::numeric_limits::infinity()) williamr@2: return std::complex(x, (z.imag() < 0) ? std::numeric_limits::infinity() : -std::numeric_limits::infinity()); williamr@2: return std::complex(x, x); williamr@2: } williamr@2: else if(std::numeric_limits::has_infinity && (y == std::numeric_limits::infinity())) williamr@2: { williamr@2: real = half_pi; williamr@2: imag = std::numeric_limits::infinity(); williamr@2: } williamr@2: else if(detail::test_is_nan(y)) williamr@2: { williamr@2: return std::complex((x == 0) ? half_pi : y, y); williamr@2: } williamr@2: else williamr@2: { williamr@2: // williamr@2: // What follows is the regular Hull et al code, williamr@2: // begin with the special case for real numbers: williamr@2: // williamr@2: if((y == 0) && (x <= one)) williamr@2: return std::complex((x == 0) ? half_pi : std::acos(z.real())); williamr@2: // williamr@2: // Figure out if our input is within the "safe area" identified by Hull et al. williamr@2: // This would be more efficient with portable floating point exception handling; williamr@2: // fortunately the quantities M and u identified by Hull et al (figure 3), williamr@2: // match with the max and min methods of numeric_limits. williamr@2: // williamr@2: T safe_max = detail::safe_max(static_cast(8)); williamr@2: T safe_min = detail::safe_min(static_cast(4)); williamr@2: williamr@2: T xp1 = one + x; williamr@2: T xm1 = x - one; williamr@2: williamr@2: if((x < safe_max) && (x > safe_min) && (y < safe_max) && (y > safe_min)) williamr@2: { williamr@2: T yy = y * y; williamr@2: T r = std::sqrt(xp1*xp1 + yy); williamr@2: T s = std::sqrt(xm1*xm1 + yy); williamr@2: T a = half * (r + s); williamr@2: T b = x / a; williamr@2: williamr@2: if(b <= b_crossover) williamr@2: { williamr@2: real = std::acos(b); williamr@2: } williamr@2: else williamr@2: { williamr@2: T apx = a + x; williamr@2: if(x <= one) williamr@2: { williamr@2: real = std::atan(std::sqrt(half * apx * (yy /(r + xp1) + (s-xm1)))/x); williamr@2: } williamr@2: else williamr@2: { williamr@2: real = std::atan((y * std::sqrt(half * (apx/(r + xp1) + apx/(s+xm1))))/x); williamr@2: } williamr@2: } williamr@2: williamr@2: if(a <= a_crossover) williamr@2: { williamr@2: T am1; williamr@2: if(x < one) williamr@2: { williamr@2: am1 = half * (yy/(r + xp1) + yy/(s - xm1)); williamr@2: } williamr@2: else williamr@2: { williamr@2: am1 = half * (yy/(r + xp1) + (s + xm1)); williamr@2: } williamr@2: imag = boost::math::log1p(am1 + std::sqrt(am1 * (a + one))); williamr@2: } williamr@2: else williamr@2: { williamr@2: imag = std::log(a + std::sqrt(a*a - one)); williamr@2: } williamr@2: } williamr@2: else williamr@2: { williamr@2: // williamr@2: // This is the Hull et al exception handling code from Fig 6 of their paper: williamr@2: // williamr@2: if(y <= (std::numeric_limits::epsilon() * std::fabs(xm1))) williamr@2: { williamr@2: if(x < one) williamr@2: { williamr@2: real = std::acos(x); williamr@2: imag = y / std::sqrt(xp1*(one-x)); williamr@2: } williamr@2: else williamr@2: { williamr@2: real = 0; williamr@2: if(((std::numeric_limits::max)() / xp1) > xm1) williamr@2: { williamr@2: // xp1 * xm1 won't overflow: williamr@2: imag = boost::math::log1p(xm1 + std::sqrt(xp1*xm1)); williamr@2: } williamr@2: else williamr@2: { williamr@2: imag = log_two + std::log(x); williamr@2: } williamr@2: } williamr@2: } williamr@2: else if(y <= safe_min) williamr@2: { williamr@2: // There is an assumption in Hull et al's analysis that williamr@2: // if we get here then x == 1. This is true for all "good" williamr@2: // machines where : williamr@2: // williamr@2: // E^2 > 8*sqrt(u); with: williamr@2: // williamr@2: // E = std::numeric_limits::epsilon() williamr@2: // u = (std::numeric_limits::min)() williamr@2: // williamr@2: // Hull et al provide alternative code for "bad" machines williamr@2: // but we have no way to test that here, so for now just assert williamr@2: // on the assumption: williamr@2: // williamr@2: BOOST_ASSERT(x == 1); williamr@2: real = std::sqrt(y); williamr@2: imag = std::sqrt(y); williamr@2: } williamr@2: else if(std::numeric_limits::epsilon() * y - one >= x) williamr@2: { williamr@2: real = half_pi; williamr@2: imag = log_two + std::log(y); williamr@2: } williamr@2: else if(x > one) williamr@2: { williamr@2: real = std::atan(y/x); williamr@2: T xoy = x/y; williamr@2: imag = log_two + std::log(y) + half * boost::math::log1p(xoy*xoy); williamr@2: } williamr@2: else williamr@2: { williamr@2: real = half_pi; williamr@2: T a = std::sqrt(one + y*y); williamr@2: imag = half * boost::math::log1p(static_cast(2)*y*(y+a)); williamr@2: } williamr@2: } williamr@2: } williamr@2: williamr@2: // williamr@2: // Finish off by working out the sign of the result: williamr@2: // williamr@2: if(z.real() < 0) williamr@2: real = s_pi - real; williamr@2: if(z.imag() > 0) williamr@2: imag = -imag; williamr@2: williamr@2: return std::complex(real, imag); williamr@2: } williamr@2: williamr@2: } } // namespaces williamr@2: williamr@2: #endif // BOOST_MATH_COMPLEX_ACOS_INCLUDED