1.1 --- a/epoc32/include/stdapis/boost/math/complex/atanh.hpp Wed Mar 31 12:27:01 2010 +0100
1.2 +++ b/epoc32/include/stdapis/boost/math/complex/atanh.hpp Wed Mar 31 12:33:34 2010 +0100
1.3 @@ -1,267 +1,245 @@
1.4 -// boost atanh.hpp header file
1.5 +// (C) Copyright John Maddock 2005.
1.6 +// Use, modification and distribution are subject to the
1.7 +// Boost Software License, Version 1.0. (See accompanying file
1.8 +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
1.9
1.10 -// (C) Copyright Hubert Holin 2001.
1.11 -// Distributed under the Boost Software License, Version 1.0. (See
1.12 -// accompanying file LICENSE_1_0.txt or copy at
1.13 -// http://www.boost.org/LICENSE_1_0.txt)
1.14 +#ifndef BOOST_MATH_COMPLEX_ATANH_INCLUDED
1.15 +#define BOOST_MATH_COMPLEX_ATANH_INCLUDED
1.16
1.17 -// See http://www.boost.org for updates, documentation, and revision history.
1.18 +#ifndef BOOST_MATH_COMPLEX_DETAILS_INCLUDED
1.19 +# include <boost/math/complex/details.hpp>
1.20 +#endif
1.21 +#ifndef BOOST_MATH_LOG1P_INCLUDED
1.22 +# include <boost/math/special_functions/log1p.hpp>
1.23 +#endif
1.24 +#include <boost/assert.hpp>
1.25
1.26 -#ifndef BOOST_ATANH_HPP
1.27 -#define BOOST_ATANH_HPP
1.28 +#ifdef BOOST_NO_STDC_NAMESPACE
1.29 +namespace std{ using ::sqrt; using ::fabs; using ::acos; using ::asin; using ::atan; using ::atan2; }
1.30 +#endif
1.31
1.32 +namespace boost{ namespace math{
1.33
1.34 -#include <cmath>
1.35 -#include <limits>
1.36 -#include <string>
1.37 -#include <stdexcept>
1.38 +template<class T>
1.39 +std::complex<T> atanh(const std::complex<T>& z)
1.40 +{
1.41 + //
1.42 + // References:
1.43 + //
1.44 + // Eric W. Weisstein. "Inverse Hyperbolic Tangent."
1.45 + // From MathWorld--A Wolfram Web Resource.
1.46 + // http://mathworld.wolfram.com/InverseHyperbolicTangent.html
1.47 + //
1.48 + // Also: The Wolfram Functions Site,
1.49 + // http://functions.wolfram.com/ElementaryFunctions/ArcTanh/
1.50 + //
1.51 + // Also "Abramowitz and Stegun. Handbook of Mathematical Functions."
1.52 + // at : http://jove.prohosting.com/~skripty/toc.htm
1.53 + //
1.54 +
1.55 + static const T half_pi = static_cast<T>(1.57079632679489661923132169163975144L);
1.56 + static const T pi = static_cast<T>(3.141592653589793238462643383279502884197L);
1.57 + static const T one = static_cast<T>(1.0L);
1.58 + static const T two = static_cast<T>(2.0L);
1.59 + static const T four = static_cast<T>(4.0L);
1.60 + static const T zero = static_cast<T>(0);
1.61 + static const T a_crossover = static_cast<T>(0.3L);
1.62
1.63 + T x = std::fabs(z.real());
1.64 + T y = std::fabs(z.imag());
1.65
1.66 -#include <boost/config.hpp>
1.67 + T real, imag; // our results
1.68
1.69 + T safe_upper = detail::safe_max(two);
1.70 + T safe_lower = detail::safe_min(static_cast<T>(2));
1.71
1.72 -// This is the inverse of the hyperbolic tangent function.
1.73 + //
1.74 + // Begin by handling the special cases specified in C99:
1.75 + //
1.76 + if(detail::test_is_nan(x))
1.77 + {
1.78 + if(detail::test_is_nan(y))
1.79 + return std::complex<T>(x, x);
1.80 + else if(std::numeric_limits<T>::has_infinity && (y == std::numeric_limits<T>::infinity()))
1.81 + return std::complex<T>(0, ((z.imag() < 0) ? -half_pi : half_pi));
1.82 + else
1.83 + return std::complex<T>(x, x);
1.84 + }
1.85 + else if(detail::test_is_nan(y))
1.86 + {
1.87 + if(x == 0)
1.88 + return std::complex<T>(x, y);
1.89 + if(std::numeric_limits<T>::has_infinity && (x == std::numeric_limits<T>::infinity()))
1.90 + return std::complex<T>(0, y);
1.91 + else
1.92 + return std::complex<T>(y, y);
1.93 + }
1.94 + else if((x > safe_lower) && (x < safe_upper) && (y > safe_lower) && (y < safe_upper))
1.95 + {
1.96
1.97 -namespace boost
1.98 -{
1.99 - namespace math
1.100 - {
1.101 -#if defined(__GNUC__) && (__GNUC__ < 3)
1.102 - // gcc 2.x ignores function scope using declarations,
1.103 - // put them in the scope of the enclosing namespace instead:
1.104 -
1.105 - using ::std::abs;
1.106 - using ::std::sqrt;
1.107 - using ::std::log;
1.108 -
1.109 - using ::std::numeric_limits;
1.110 -#endif
1.111 -
1.112 -#if defined(BOOST_NO_TEMPLATE_PARTIAL_SPECIALIZATION)
1.113 - // This is the main fare
1.114 -
1.115 - template<typename T>
1.116 - inline T atanh(const T x)
1.117 - {
1.118 - using ::std::abs;
1.119 - using ::std::sqrt;
1.120 - using ::std::log;
1.121 -
1.122 - using ::std::numeric_limits;
1.123 -
1.124 - T const one = static_cast<T>(1);
1.125 - T const two = static_cast<T>(2);
1.126 -
1.127 - static T const taylor_2_bound = sqrt(numeric_limits<T>::epsilon());
1.128 - static T const taylor_n_bound = sqrt(taylor_2_bound);
1.129 -
1.130 - if (x < -one)
1.131 - {
1.132 - if (numeric_limits<T>::has_quiet_NaN)
1.133 - {
1.134 - return(numeric_limits<T>::quiet_NaN());
1.135 - }
1.136 - else
1.137 - {
1.138 - ::std::string error_reporting("Argument to atanh is strictly greater than +1 or strictly smaller than -1!");
1.139 - ::std::domain_error bad_argument(error_reporting);
1.140 -
1.141 - throw(bad_argument);
1.142 - }
1.143 - }
1.144 - else if (x < -one+numeric_limits<T>::epsilon())
1.145 - {
1.146 - if (numeric_limits<T>::has_infinity)
1.147 - {
1.148 - return(-numeric_limits<T>::infinity());
1.149 - }
1.150 - else
1.151 - {
1.152 - ::std::string error_reporting("Argument to atanh is -1 (result: -Infinity)!");
1.153 - ::std::out_of_range bad_argument(error_reporting);
1.154 -
1.155 - throw(bad_argument);
1.156 - }
1.157 - }
1.158 - else if (x > +one-numeric_limits<T>::epsilon())
1.159 - {
1.160 - if (numeric_limits<T>::has_infinity)
1.161 - {
1.162 - return(+numeric_limits<T>::infinity());
1.163 - }
1.164 - else
1.165 - {
1.166 - ::std::string error_reporting("Argument to atanh is +1 (result: +Infinity)!");
1.167 - ::std::out_of_range bad_argument(error_reporting);
1.168 -
1.169 - throw(bad_argument);
1.170 - }
1.171 - }
1.172 - else if (x > +one)
1.173 - {
1.174 - if (numeric_limits<T>::has_quiet_NaN)
1.175 - {
1.176 - return(numeric_limits<T>::quiet_NaN());
1.177 - }
1.178 - else
1.179 - {
1.180 - ::std::string error_reporting("Argument to atanh is strictly greater than +1 or strictly smaller than -1!");
1.181 - ::std::domain_error bad_argument(error_reporting);
1.182 -
1.183 - throw(bad_argument);
1.184 - }
1.185 - }
1.186 - else if (abs(x) >= taylor_n_bound)
1.187 - {
1.188 - return(log( (one + x) / (one - x) ) / two);
1.189 - }
1.190 + T xx = x*x;
1.191 + T yy = y*y;
1.192 + T x2 = x * two;
1.193 +
1.194 + ///
1.195 + // The real part is given by:
1.196 + //
1.197 + // real(atanh(z)) == log((1 + x^2 + y^2 + 2x) / (1 + x^2 + y^2 - 2x))
1.198 + //
1.199 + // However, when x is either large (x > 1/E) or very small
1.200 + // (x < E) then this effectively simplifies
1.201 + // to log(1), leading to wildly inaccurate results.
1.202 + // By dividing the above (top and bottom) by (1 + x^2 + y^2) we get:
1.203 + //
1.204 + // real(atanh(z)) == log((1 + (2x / (1 + x^2 + y^2))) / (1 - (-2x / (1 + x^2 + y^2))))
1.205 + //
1.206 + // which is much more sensitive to the value of x, when x is not near 1
1.207 + // (remember we can compute log(1+x) for small x very accurately).
1.208 + //
1.209 + // The cross-over from one method to the other has to be determined
1.210 + // experimentally, the value used below appears correct to within a
1.211 + // factor of 2 (and there are larger errors from other parts
1.212 + // of the input domain anyway).
1.213 + //
1.214 + T alpha = two*x / (one + xx + yy);
1.215 + if(alpha < a_crossover)
1.216 + {
1.217 + real = boost::math::log1p(alpha) - boost::math::log1p(-alpha);
1.218 + }
1.219 + else
1.220 + {
1.221 + T xm1 = x - one;
1.222 + real = boost::math::log1p(x2 + xx + yy) - std::log(xm1*xm1 + yy);
1.223 + }
1.224 + real /= four;
1.225 + if(z.real() < 0)
1.226 + real = -real;
1.227 +
1.228 + imag = std::atan2((y * two), (one - xx - yy));
1.229 + imag /= two;
1.230 + if(z.imag() < 0)
1.231 + imag = -imag;
1.232 + }
1.233 + else
1.234 + {
1.235 + //
1.236 + // This section handles exception cases that would normally cause
1.237 + // underflow or overflow in the main formulas.
1.238 + //
1.239 + // Begin by working out the real part, we need to approximate
1.240 + // alpha = 2x / (1 + x^2 + y^2)
1.241 + // without either overflow or underflow in the squared terms.
1.242 + //
1.243 + T alpha = 0;
1.244 + if(x >= safe_upper)
1.245 + {
1.246 + // this is really a test for infinity,
1.247 + // but we may not have the necessary numeric_limits support:
1.248 + if((x > (std::numeric_limits<T>::max)()) || (y > (std::numeric_limits<T>::max)()))
1.249 + {
1.250 + alpha = 0;
1.251 + }
1.252 + else if(y >= safe_upper)
1.253 + {
1.254 + // Big x and y: divide alpha through by x*y:
1.255 + alpha = (two/y) / (x/y + y/x);
1.256 + }
1.257 + else if(y > one)
1.258 + {
1.259 + // Big x: divide through by x:
1.260 + alpha = two / (x + y*y/x);
1.261 + }
1.262 + else
1.263 + {
1.264 + // Big x small y, as above but neglect y^2/x:
1.265 + alpha = two/x;
1.266 + }
1.267 + }
1.268 + else if(y >= safe_upper)
1.269 + {
1.270 + if(x > one)
1.271 + {
1.272 + // Big y, medium x, divide through by y:
1.273 + alpha = (two*x/y) / (y + x*x/y);
1.274 + }
1.275 + else
1.276 + {
1.277 + // Small x and y, whatever alpha is, it's too small to calculate:
1.278 + alpha = 0;
1.279 + }
1.280 + }
1.281 + else
1.282 + {
1.283 + // one or both of x and y are small, calculate divisor carefully:
1.284 + T div = one;
1.285 + if(x > safe_lower)
1.286 + div += x*x;
1.287 + if(y > safe_lower)
1.288 + div += y*y;
1.289 + alpha = two*x/div;
1.290 + }
1.291 + if(alpha < a_crossover)
1.292 + {
1.293 + real = boost::math::log1p(alpha) - boost::math::log1p(-alpha);
1.294 + }
1.295 + else
1.296 + {
1.297 + // We can only get here as a result of small y and medium sized x,
1.298 + // we can simply neglect the y^2 terms:
1.299 + BOOST_ASSERT(x >= safe_lower);
1.300 + BOOST_ASSERT(x <= safe_upper);
1.301 + //BOOST_ASSERT(y <= safe_lower);
1.302 + T xm1 = x - one;
1.303 + real = std::log(1 + two*x + x*x) - std::log(xm1*xm1);
1.304 + }
1.305 +
1.306 + real /= four;
1.307 + if(z.real() < 0)
1.308 + real = -real;
1.309 +
1.310 + //
1.311 + // Now handle imaginary part, this is much easier,
1.312 + // if x or y are large, then the formula:
1.313 + // atan2(2y, 1 - x^2 - y^2)
1.314 + // evaluates to +-(PI - theta) where theta is negligible compared to PI.
1.315 + //
1.316 + if((x >= safe_upper) || (y >= safe_upper))
1.317 + {
1.318 + imag = pi;
1.319 + }
1.320 + else if(x <= safe_lower)
1.321 + {
1.322 + //
1.323 + // If both x and y are small then atan(2y),
1.324 + // otherwise just x^2 is negligible in the divisor:
1.325 + //
1.326 + if(y <= safe_lower)
1.327 + imag = std::atan2(two*y, one);
1.328 + else
1.329 + {
1.330 + if((y == zero) && (x == zero))
1.331 + imag = 0;
1.332 else
1.333 - {
1.334 - // approximation by taylor series in x at 0 up to order 2
1.335 - T result = x;
1.336 -
1.337 - if (abs(x) >= taylor_2_bound)
1.338 - {
1.339 - T x3 = x*x*x;
1.340 -
1.341 - // approximation by taylor series in x at 0 up to order 4
1.342 - result += x3/static_cast<T>(3);
1.343 - }
1.344 -
1.345 - return(result);
1.346 - }
1.347 - }
1.348 -#else
1.349 - // These are implementation details (for main fare see below)
1.350 -
1.351 - namespace detail
1.352 - {
1.353 - template <
1.354 - typename T,
1.355 - bool InfinitySupported
1.356 - >
1.357 - struct atanh_helper1_t
1.358 - {
1.359 - static T get_pos_infinity()
1.360 - {
1.361 - return(+::std::numeric_limits<T>::infinity());
1.362 - }
1.363 -
1.364 - static T get_neg_infinity()
1.365 - {
1.366 - return(-::std::numeric_limits<T>::infinity());
1.367 - }
1.368 - }; // boost::math::detail::atanh_helper1_t
1.369 -
1.370 -
1.371 - template<typename T>
1.372 - struct atanh_helper1_t<T, false>
1.373 - {
1.374 - static T get_pos_infinity()
1.375 - {
1.376 - ::std::string error_reporting("Argument to atanh is +1 (result: +Infinity)!");
1.377 - ::std::out_of_range bad_argument(error_reporting);
1.378 -
1.379 - throw(bad_argument);
1.380 - }
1.381 -
1.382 - static T get_neg_infinity()
1.383 - {
1.384 - ::std::string error_reporting("Argument to atanh is -1 (result: -Infinity)!");
1.385 - ::std::out_of_range bad_argument(error_reporting);
1.386 -
1.387 - throw(bad_argument);
1.388 - }
1.389 - }; // boost::math::detail::atanh_helper1_t
1.390 -
1.391 -
1.392 - template <
1.393 - typename T,
1.394 - bool QuietNanSupported
1.395 - >
1.396 - struct atanh_helper2_t
1.397 - {
1.398 - static T get_NaN()
1.399 - {
1.400 - return(::std::numeric_limits<T>::quiet_NaN());
1.401 - }
1.402 - }; // boost::detail::atanh_helper2_t
1.403 -
1.404 -
1.405 - template<typename T>
1.406 - struct atanh_helper2_t<T, false>
1.407 - {
1.408 - static T get_NaN()
1.409 - {
1.410 - ::std::string error_reporting("Argument to atanh is strictly greater than +1 or strictly smaller than -1!");
1.411 - ::std::domain_error bad_argument(error_reporting);
1.412 -
1.413 - throw(bad_argument);
1.414 - }
1.415 - }; // boost::detail::atanh_helper2_t
1.416 - } // boost::detail
1.417 -
1.418 -
1.419 - // This is the main fare
1.420 -
1.421 - template<typename T>
1.422 - inline T atanh(const T x)
1.423 - {
1.424 - using ::std::abs;
1.425 - using ::std::sqrt;
1.426 - using ::std::log;
1.427 -
1.428 - using ::std::numeric_limits;
1.429 -
1.430 - typedef detail::atanh_helper1_t<T, ::std::numeric_limits<T>::has_infinity> helper1_type;
1.431 - typedef detail::atanh_helper2_t<T, ::std::numeric_limits<T>::has_quiet_NaN> helper2_type;
1.432 -
1.433 -
1.434 - T const one = static_cast<T>(1);
1.435 - T const two = static_cast<T>(2);
1.436 -
1.437 - static T const taylor_2_bound = sqrt(numeric_limits<T>::epsilon());
1.438 - static T const taylor_n_bound = sqrt(taylor_2_bound);
1.439 -
1.440 - if (x < -one)
1.441 - {
1.442 - return(helper2_type::get_NaN());
1.443 - }
1.444 - else if (x < -one+numeric_limits<T>::epsilon())
1.445 - {
1.446 - return(helper1_type::get_neg_infinity());
1.447 - }
1.448 - else if (x > +one-numeric_limits<T>::epsilon())
1.449 - {
1.450 - return(helper1_type::get_pos_infinity());
1.451 - }
1.452 - else if (x > +one)
1.453 - {
1.454 - return(helper2_type::get_NaN());
1.455 - }
1.456 - else if (abs(x) >= taylor_n_bound)
1.457 - {
1.458 - return(log( (one + x) / (one - x) ) / two);
1.459 - }
1.460 - else
1.461 - {
1.462 - // approximation by taylor series in x at 0 up to order 2
1.463 - T result = x;
1.464 -
1.465 - if (abs(x) >= taylor_2_bound)
1.466 - {
1.467 - T x3 = x*x*x;
1.468 -
1.469 - // approximation by taylor series in x at 0 up to order 4
1.470 - result += x3/static_cast<T>(3);
1.471 - }
1.472 -
1.473 - return(result);
1.474 - }
1.475 - }
1.476 -#endif /* defined(BOOST_NO_TEMPLATE_PARTIAL_SPECIALIZATION) */
1.477 - }
1.478 + imag = std::atan2(two*y, one - y*y);
1.479 + }
1.480 + }
1.481 + else
1.482 + {
1.483 + //
1.484 + // y^2 is negligible:
1.485 + //
1.486 + if((y == zero) && (x == one))
1.487 + imag = 0;
1.488 + else
1.489 + imag = std::atan2(two*y, 1 - x*x);
1.490 + }
1.491 + imag /= two;
1.492 + if(z.imag() < 0)
1.493 + imag = -imag;
1.494 + }
1.495 + return std::complex<T>(real, imag);
1.496 }
1.497
1.498 -#endif /* BOOST_ATANH_HPP */
1.499 +} } // namespaces
1.500
1.501 +#endif // BOOST_MATH_COMPLEX_ATANH_INCLUDED