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