1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/os/ossrv/ossrv_pub/boost_apis/boost/math/complex/atanh.hpp Fri Jun 15 03:10:57 2012 +0200
1.3 @@ -0,0 +1,245 @@
1.4 +// (C) Copyright John Maddock 2005.
1.5 +// Use, modification and distribution are subject to the
1.6 +// Boost Software License, Version 1.0. (See accompanying file
1.7 +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
1.8 +
1.9 +#ifndef BOOST_MATH_COMPLEX_ATANH_INCLUDED
1.10 +#define BOOST_MATH_COMPLEX_ATANH_INCLUDED
1.11 +
1.12 +#ifndef BOOST_MATH_COMPLEX_DETAILS_INCLUDED
1.13 +# include <boost/math/complex/details.hpp>
1.14 +#endif
1.15 +#ifndef BOOST_MATH_LOG1P_INCLUDED
1.16 +# include <boost/math/special_functions/log1p.hpp>
1.17 +#endif
1.18 +#include <boost/assert.hpp>
1.19 +
1.20 +#ifdef BOOST_NO_STDC_NAMESPACE
1.21 +namespace std{ using ::sqrt; using ::fabs; using ::acos; using ::asin; using ::atan; using ::atan2; }
1.22 +#endif
1.23 +
1.24 +namespace boost{ namespace math{
1.25 +
1.26 +template<class T>
1.27 +std::complex<T> atanh(const std::complex<T>& z)
1.28 +{
1.29 + //
1.30 + // References:
1.31 + //
1.32 + // Eric W. Weisstein. "Inverse Hyperbolic Tangent."
1.33 + // From MathWorld--A Wolfram Web Resource.
1.34 + // http://mathworld.wolfram.com/InverseHyperbolicTangent.html
1.35 + //
1.36 + // Also: The Wolfram Functions Site,
1.37 + // http://functions.wolfram.com/ElementaryFunctions/ArcTanh/
1.38 + //
1.39 + // Also "Abramowitz and Stegun. Handbook of Mathematical Functions."
1.40 + // at : http://jove.prohosting.com/~skripty/toc.htm
1.41 + //
1.42 +
1.43 + static const T half_pi = static_cast<T>(1.57079632679489661923132169163975144L);
1.44 + static const T pi = static_cast<T>(3.141592653589793238462643383279502884197L);
1.45 + static const T one = static_cast<T>(1.0L);
1.46 + static const T two = static_cast<T>(2.0L);
1.47 + static const T four = static_cast<T>(4.0L);
1.48 + static const T zero = static_cast<T>(0);
1.49 + static const T a_crossover = static_cast<T>(0.3L);
1.50 +
1.51 + T x = std::fabs(z.real());
1.52 + T y = std::fabs(z.imag());
1.53 +
1.54 + T real, imag; // our results
1.55 +
1.56 + T safe_upper = detail::safe_max(two);
1.57 + T safe_lower = detail::safe_min(static_cast<T>(2));
1.58 +
1.59 + //
1.60 + // Begin by handling the special cases specified in C99:
1.61 + //
1.62 + if(detail::test_is_nan(x))
1.63 + {
1.64 + if(detail::test_is_nan(y))
1.65 + return std::complex<T>(x, x);
1.66 + else if(std::numeric_limits<T>::has_infinity && (y == std::numeric_limits<T>::infinity()))
1.67 + return std::complex<T>(0, ((z.imag() < 0) ? -half_pi : half_pi));
1.68 + else
1.69 + return std::complex<T>(x, x);
1.70 + }
1.71 + else if(detail::test_is_nan(y))
1.72 + {
1.73 + if(x == 0)
1.74 + return std::complex<T>(x, y);
1.75 + if(std::numeric_limits<T>::has_infinity && (x == std::numeric_limits<T>::infinity()))
1.76 + return std::complex<T>(0, y);
1.77 + else
1.78 + return std::complex<T>(y, y);
1.79 + }
1.80 + else if((x > safe_lower) && (x < safe_upper) && (y > safe_lower) && (y < safe_upper))
1.81 + {
1.82 +
1.83 + T xx = x*x;
1.84 + T yy = y*y;
1.85 + T x2 = x * two;
1.86 +
1.87 + ///
1.88 + // The real part is given by:
1.89 + //
1.90 + // real(atanh(z)) == log((1 + x^2 + y^2 + 2x) / (1 + x^2 + y^2 - 2x))
1.91 + //
1.92 + // However, when x is either large (x > 1/E) or very small
1.93 + // (x < E) then this effectively simplifies
1.94 + // to log(1), leading to wildly inaccurate results.
1.95 + // By dividing the above (top and bottom) by (1 + x^2 + y^2) we get:
1.96 + //
1.97 + // real(atanh(z)) == log((1 + (2x / (1 + x^2 + y^2))) / (1 - (-2x / (1 + x^2 + y^2))))
1.98 + //
1.99 + // which is much more sensitive to the value of x, when x is not near 1
1.100 + // (remember we can compute log(1+x) for small x very accurately).
1.101 + //
1.102 + // The cross-over from one method to the other has to be determined
1.103 + // experimentally, the value used below appears correct to within a
1.104 + // factor of 2 (and there are larger errors from other parts
1.105 + // of the input domain anyway).
1.106 + //
1.107 + T alpha = two*x / (one + xx + yy);
1.108 + if(alpha < a_crossover)
1.109 + {
1.110 + real = boost::math::log1p(alpha) - boost::math::log1p(-alpha);
1.111 + }
1.112 + else
1.113 + {
1.114 + T xm1 = x - one;
1.115 + real = boost::math::log1p(x2 + xx + yy) - std::log(xm1*xm1 + yy);
1.116 + }
1.117 + real /= four;
1.118 + if(z.real() < 0)
1.119 + real = -real;
1.120 +
1.121 + imag = std::atan2((y * two), (one - xx - yy));
1.122 + imag /= two;
1.123 + if(z.imag() < 0)
1.124 + imag = -imag;
1.125 + }
1.126 + else
1.127 + {
1.128 + //
1.129 + // This section handles exception cases that would normally cause
1.130 + // underflow or overflow in the main formulas.
1.131 + //
1.132 + // Begin by working out the real part, we need to approximate
1.133 + // alpha = 2x / (1 + x^2 + y^2)
1.134 + // without either overflow or underflow in the squared terms.
1.135 + //
1.136 + T alpha = 0;
1.137 + if(x >= safe_upper)
1.138 + {
1.139 + // this is really a test for infinity,
1.140 + // but we may not have the necessary numeric_limits support:
1.141 + if((x > (std::numeric_limits<T>::max)()) || (y > (std::numeric_limits<T>::max)()))
1.142 + {
1.143 + alpha = 0;
1.144 + }
1.145 + else if(y >= safe_upper)
1.146 + {
1.147 + // Big x and y: divide alpha through by x*y:
1.148 + alpha = (two/y) / (x/y + y/x);
1.149 + }
1.150 + else if(y > one)
1.151 + {
1.152 + // Big x: divide through by x:
1.153 + alpha = two / (x + y*y/x);
1.154 + }
1.155 + else
1.156 + {
1.157 + // Big x small y, as above but neglect y^2/x:
1.158 + alpha = two/x;
1.159 + }
1.160 + }
1.161 + else if(y >= safe_upper)
1.162 + {
1.163 + if(x > one)
1.164 + {
1.165 + // Big y, medium x, divide through by y:
1.166 + alpha = (two*x/y) / (y + x*x/y);
1.167 + }
1.168 + else
1.169 + {
1.170 + // Small x and y, whatever alpha is, it's too small to calculate:
1.171 + alpha = 0;
1.172 + }
1.173 + }
1.174 + else
1.175 + {
1.176 + // one or both of x and y are small, calculate divisor carefully:
1.177 + T div = one;
1.178 + if(x > safe_lower)
1.179 + div += x*x;
1.180 + if(y > safe_lower)
1.181 + div += y*y;
1.182 + alpha = two*x/div;
1.183 + }
1.184 + if(alpha < a_crossover)
1.185 + {
1.186 + real = boost::math::log1p(alpha) - boost::math::log1p(-alpha);
1.187 + }
1.188 + else
1.189 + {
1.190 + // We can only get here as a result of small y and medium sized x,
1.191 + // we can simply neglect the y^2 terms:
1.192 + BOOST_ASSERT(x >= safe_lower);
1.193 + BOOST_ASSERT(x <= safe_upper);
1.194 + //BOOST_ASSERT(y <= safe_lower);
1.195 + T xm1 = x - one;
1.196 + real = std::log(1 + two*x + x*x) - std::log(xm1*xm1);
1.197 + }
1.198 +
1.199 + real /= four;
1.200 + if(z.real() < 0)
1.201 + real = -real;
1.202 +
1.203 + //
1.204 + // Now handle imaginary part, this is much easier,
1.205 + // if x or y are large, then the formula:
1.206 + // atan2(2y, 1 - x^2 - y^2)
1.207 + // evaluates to +-(PI - theta) where theta is negligible compared to PI.
1.208 + //
1.209 + if((x >= safe_upper) || (y >= safe_upper))
1.210 + {
1.211 + imag = pi;
1.212 + }
1.213 + else if(x <= safe_lower)
1.214 + {
1.215 + //
1.216 + // If both x and y are small then atan(2y),
1.217 + // otherwise just x^2 is negligible in the divisor:
1.218 + //
1.219 + if(y <= safe_lower)
1.220 + imag = std::atan2(two*y, one);
1.221 + else
1.222 + {
1.223 + if((y == zero) && (x == zero))
1.224 + imag = 0;
1.225 + else
1.226 + imag = std::atan2(two*y, one - y*y);
1.227 + }
1.228 + }
1.229 + else
1.230 + {
1.231 + //
1.232 + // y^2 is negligible:
1.233 + //
1.234 + if((y == zero) && (x == one))
1.235 + imag = 0;
1.236 + else
1.237 + imag = std::atan2(two*y, 1 - x*x);
1.238 + }
1.239 + imag /= two;
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_ATANH_INCLUDED