epoc32/include/stdapis/boost/math/complex/atanh.hpp
branchSymbian3
changeset 4 837f303aceeb
parent 3 e1b950c65cb4
     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