1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/epoc32/include/stdapis/boost/math/complex/asin.hpp Tue Mar 16 16:12:26 2010 +0000
1.3 @@ -0,0 +1,245 @@
1.4 +// (C) Copyright John Maddock 2005.
1.5 +// Distributed under the Boost Software License, Version 1.0. (See accompanying
1.6 +// file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
1.7 +
1.8 +#ifndef BOOST_MATH_COMPLEX_ASIN_INCLUDED
1.9 +#define BOOST_MATH_COMPLEX_ASIN_INCLUDED
1.10 +
1.11 +#ifndef BOOST_MATH_COMPLEX_DETAILS_INCLUDED
1.12 +# include <boost/math/complex/details.hpp>
1.13 +#endif
1.14 +#ifndef BOOST_MATH_LOG1P_INCLUDED
1.15 +# include <boost/math/special_functions/log1p.hpp>
1.16 +#endif
1.17 +#include <boost/assert.hpp>
1.18 +
1.19 +#ifdef BOOST_NO_STDC_NAMESPACE
1.20 +namespace std{ using ::sqrt; using ::fabs; using ::acos; using ::asin; using ::atan; using ::atan2; }
1.21 +#endif
1.22 +
1.23 +namespace boost{ namespace math{
1.24 +
1.25 +template<class T>
1.26 +inline std::complex<T> asin(const std::complex<T>& z)
1.27 +{
1.28 + //
1.29 + // This implementation is a transcription of the pseudo-code in:
1.30 + //
1.31 + // "Implementing the complex Arcsine and Arccosine Functions using Exception Handling."
1.32 + // T E Hull, Thomas F Fairgrieve and Ping Tak Peter Tang.
1.33 + // ACM Transactions on Mathematical Software, Vol 23, No 3, Sept 1997.
1.34 + //
1.35 +
1.36 + //
1.37 + // These static constants should really be in a maths constants library:
1.38 + //
1.39 + static const T one = static_cast<T>(1);
1.40 + //static const T two = static_cast<T>(2);
1.41 + static const T half = static_cast<T>(0.5L);
1.42 + static const T a_crossover = static_cast<T>(1.5L);
1.43 + static const T b_crossover = static_cast<T>(0.6417L);
1.44 + //static const T pi = static_cast<T>(3.141592653589793238462643383279502884197L);
1.45 + static const T half_pi = static_cast<T>(1.57079632679489661923132169163975144L);
1.46 + static const T log_two = static_cast<T>(0.69314718055994530941723212145817657L);
1.47 + static const T quarter_pi = static_cast<T>(0.78539816339744830961566084581987572L);
1.48 +
1.49 + //
1.50 + // Get real and imaginary parts, discard the signs as we can
1.51 + // figure out the sign of the result later:
1.52 + //
1.53 + T x = std::fabs(z.real());
1.54 + T y = std::fabs(z.imag());
1.55 + T real, imag; // our results
1.56 +
1.57 + //
1.58 + // Begin by handling the special cases for infinities and nan's
1.59 + // specified in C99, most of this is handled by the regular logic
1.60 + // below, but handling it as a special case prevents overflow/underflow
1.61 + // arithmetic which may trip up some machines:
1.62 + //
1.63 + if(detail::test_is_nan(x))
1.64 + {
1.65 + if(detail::test_is_nan(y))
1.66 + return std::complex<T>(x, x);
1.67 + if(std::numeric_limits<T>::has_infinity && (y == std::numeric_limits<T>::infinity()))
1.68 + {
1.69 + real = x;
1.70 + imag = std::numeric_limits<T>::infinity();
1.71 + }
1.72 + else
1.73 + return std::complex<T>(x, x);
1.74 + }
1.75 + else if(detail::test_is_nan(y))
1.76 + {
1.77 + if(x == 0)
1.78 + {
1.79 + real = 0;
1.80 + imag = y;
1.81 + }
1.82 + else if(std::numeric_limits<T>::has_infinity && (x == std::numeric_limits<T>::infinity()))
1.83 + {
1.84 + real = y;
1.85 + imag = std::numeric_limits<T>::infinity();
1.86 + }
1.87 + else
1.88 + return std::complex<T>(y, y);
1.89 + }
1.90 + else if(std::numeric_limits<T>::has_infinity && (x == std::numeric_limits<T>::infinity()))
1.91 + {
1.92 + if(y == std::numeric_limits<T>::infinity())
1.93 + {
1.94 + real = quarter_pi;
1.95 + imag = std::numeric_limits<T>::infinity();
1.96 + }
1.97 + else
1.98 + {
1.99 + real = half_pi;
1.100 + imag = std::numeric_limits<T>::infinity();
1.101 + }
1.102 + }
1.103 + else if(std::numeric_limits<T>::has_infinity && (y == std::numeric_limits<T>::infinity()))
1.104 + {
1.105 + real = 0;
1.106 + imag = std::numeric_limits<T>::infinity();
1.107 + }
1.108 + else
1.109 + {
1.110 + //
1.111 + // special case for real numbers:
1.112 + //
1.113 + if((y == 0) && (x <= one))
1.114 + return std::complex<T>(std::asin(z.real()));
1.115 + //
1.116 + // Figure out if our input is within the "safe area" identified by Hull et al.
1.117 + // This would be more efficient with portable floating point exception handling;
1.118 + // fortunately the quantities M and u identified by Hull et al (figure 3),
1.119 + // match with the max and min methods of numeric_limits<T>.
1.120 + //
1.121 + T safe_max = detail::safe_max(static_cast<T>(8));
1.122 + T safe_min = detail::safe_min(static_cast<T>(4));
1.123 +
1.124 + T xp1 = one + x;
1.125 + T xm1 = x - one;
1.126 +
1.127 + if((x < safe_max) && (x > safe_min) && (y < safe_max) && (y > safe_min))
1.128 + {
1.129 + T yy = y * y;
1.130 + T r = std::sqrt(xp1*xp1 + yy);
1.131 + T s = std::sqrt(xm1*xm1 + yy);
1.132 + T a = half * (r + s);
1.133 + T b = x / a;
1.134 +
1.135 + if(b <= b_crossover)
1.136 + {
1.137 + real = std::asin(b);
1.138 + }
1.139 + else
1.140 + {
1.141 + T apx = a + x;
1.142 + if(x <= one)
1.143 + {
1.144 + real = std::atan(x/std::sqrt(half * apx * (yy /(r + xp1) + (s-xm1))));
1.145 + }
1.146 + else
1.147 + {
1.148 + real = std::atan(x/(y * std::sqrt(half * (apx/(r + xp1) + apx/(s+xm1)))));
1.149 + }
1.150 + }
1.151 +
1.152 + if(a <= a_crossover)
1.153 + {
1.154 + T am1;
1.155 + if(x < one)
1.156 + {
1.157 + am1 = half * (yy/(r + xp1) + yy/(s - xm1));
1.158 + }
1.159 + else
1.160 + {
1.161 + am1 = half * (yy/(r + xp1) + (s + xm1));
1.162 + }
1.163 + imag = boost::math::log1p(am1 + std::sqrt(am1 * (a + one)));
1.164 + }
1.165 + else
1.166 + {
1.167 + imag = std::log(a + std::sqrt(a*a - one));
1.168 + }
1.169 + }
1.170 + else
1.171 + {
1.172 + //
1.173 + // This is the Hull et al exception handling code from Fig 3 of their paper:
1.174 + //
1.175 + if(y <= (std::numeric_limits<T>::epsilon() * std::fabs(xm1)))
1.176 + {
1.177 + if(x < one)
1.178 + {
1.179 + real = std::asin(x);
1.180 + imag = y / std::sqrt(xp1*xm1);
1.181 + }
1.182 + else
1.183 + {
1.184 + real = half_pi;
1.185 + if(((std::numeric_limits<T>::max)() / xp1) > xm1)
1.186 + {
1.187 + // xp1 * xm1 won't overflow:
1.188 + imag = boost::math::log1p(xm1 + std::sqrt(xp1*xm1));
1.189 + }
1.190 + else
1.191 + {
1.192 + imag = log_two + std::log(x);
1.193 + }
1.194 + }
1.195 + }
1.196 + else if(y <= safe_min)
1.197 + {
1.198 + // There is an assumption in Hull et al's analysis that
1.199 + // if we get here then x == 1. This is true for all "good"
1.200 + // machines where :
1.201 + //
1.202 + // E^2 > 8*sqrt(u); with:
1.203 + //
1.204 + // E = std::numeric_limits<T>::epsilon()
1.205 + // u = (std::numeric_limits<T>::min)()
1.206 + //
1.207 + // Hull et al provide alternative code for "bad" machines
1.208 + // but we have no way to test that here, so for now just assert
1.209 + // on the assumption:
1.210 + //
1.211 + BOOST_ASSERT(x == 1);
1.212 + real = half_pi - std::sqrt(y);
1.213 + imag = std::sqrt(y);
1.214 + }
1.215 + else if(std::numeric_limits<T>::epsilon() * y - one >= x)
1.216 + {
1.217 + real = x/y; // This can underflow!
1.218 + imag = log_two + std::log(y);
1.219 + }
1.220 + else if(x > one)
1.221 + {
1.222 + real = std::atan(x/y);
1.223 + T xoy = x/y;
1.224 + imag = log_two + std::log(y) + half * boost::math::log1p(xoy*xoy);
1.225 + }
1.226 + else
1.227 + {
1.228 + T a = std::sqrt(one + y*y);
1.229 + real = x/a; // This can underflow!
1.230 + imag = half * boost::math::log1p(static_cast<T>(2)*y*(y+a));
1.231 + }
1.232 + }
1.233 + }
1.234 +
1.235 + //
1.236 + // Finish off by working out the sign of the result:
1.237 + //
1.238 + if(z.real() < 0)
1.239 + real = -real;
1.240 + if(z.imag() < 0)
1.241 + imag = -imag;
1.242 +
1.243 + return std::complex<T>(real, imag);
1.244 +}
1.245 +
1.246 +} } // namespaces
1.247 +
1.248 +#endif // BOOST_MATH_COMPLEX_ASIN_INCLUDED