--- /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 <boost/math/complex/details.hpp>
+#endif
+#ifndef BOOST_MATH_LOG1P_INCLUDED
+# include <boost/math/special_functions/log1p.hpp>
+#endif
+#include <boost/assert.hpp>
+
+#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<class T>
+inline std::complex<T> asin(const std::complex<T>& 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<T>(1);
+ //static const T two = static_cast<T>(2);
+ static const T half = static_cast<T>(0.5L);
+ static const T a_crossover = static_cast<T>(1.5L);
+ static const T b_crossover = static_cast<T>(0.6417L);
+ //static const T pi = static_cast<T>(3.141592653589793238462643383279502884197L);
+ static const T half_pi = static_cast<T>(1.57079632679489661923132169163975144L);
+ static const T log_two = static_cast<T>(0.69314718055994530941723212145817657L);
+ static const T quarter_pi = static_cast<T>(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<T>(x, x);
+ if(std::numeric_limits<T>::has_infinity && (y == std::numeric_limits<T>::infinity()))
+ {
+ real = x;
+ imag = std::numeric_limits<T>::infinity();
+ }
+ else
+ return std::complex<T>(x, x);
+ }
+ else if(detail::test_is_nan(y))
+ {
+ if(x == 0)
+ {
+ real = 0;
+ imag = y;
+ }
+ else if(std::numeric_limits<T>::has_infinity && (x == std::numeric_limits<T>::infinity()))
+ {
+ real = y;
+ imag = std::numeric_limits<T>::infinity();
+ }
+ else
+ return std::complex<T>(y, y);
+ }
+ else if(std::numeric_limits<T>::has_infinity && (x == std::numeric_limits<T>::infinity()))
+ {
+ if(y == std::numeric_limits<T>::infinity())
+ {
+ real = quarter_pi;
+ imag = std::numeric_limits<T>::infinity();
+ }
+ else
+ {
+ real = half_pi;
+ imag = std::numeric_limits<T>::infinity();
+ }
+ }
+ else if(std::numeric_limits<T>::has_infinity && (y == std::numeric_limits<T>::infinity()))
+ {
+ real = 0;
+ imag = std::numeric_limits<T>::infinity();
+ }
+ else
+ {
+ //
+ // special case for real numbers:
+ //
+ if((y == 0) && (x <= one))
+ return std::complex<T>(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>.
+ //
+ T safe_max = detail::safe_max(static_cast<T>(8));
+ T safe_min = detail::safe_min(static_cast<T>(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<T>::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<T>::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<T>::epsilon()
+ // u = (std::numeric_limits<T>::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<T>::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<T>(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<T>(real, imag);
+}
+
+} } // namespaces
+
+#endif // BOOST_MATH_COMPLEX_ASIN_INCLUDED