williamr@2: //    boost atanh.hpp header file
williamr@2: 
williamr@2: //  (C) Copyright Hubert Holin 2001.
williamr@2: //  Distributed under the Boost Software License, Version 1.0. (See
williamr@2: //  accompanying file LICENSE_1_0.txt or copy at
williamr@2: //  http://www.boost.org/LICENSE_1_0.txt)
williamr@2: 
williamr@2: // See http://www.boost.org for updates, documentation, and revision history.
williamr@2: 
williamr@2: #ifndef BOOST_ATANH_HPP
williamr@2: #define BOOST_ATANH_HPP
williamr@2: 
williamr@2: 
williamr@2: #include <cmath>
williamr@2: #include <limits>
williamr@2: #include <string>
williamr@2: #include <stdexcept>
williamr@2: 
williamr@2: 
williamr@2: #include <boost/config.hpp>
williamr@2: 
williamr@2: 
williamr@2: // This is the inverse of the hyperbolic tangent function.
williamr@2: 
williamr@2: namespace boost
williamr@2: {
williamr@2:     namespace math
williamr@2:     {
williamr@2: #if defined(__GNUC__) && (__GNUC__ < 3)
williamr@2:         // gcc 2.x ignores function scope using declarations,
williamr@2:         // put them in the scope of the enclosing namespace instead:
williamr@2:         
williamr@2:         using    ::std::abs;
williamr@2:         using    ::std::sqrt;
williamr@2:         using    ::std::log;
williamr@2:         
williamr@2:         using    ::std::numeric_limits;
williamr@2: #endif
williamr@2:         
williamr@2: #if defined(BOOST_NO_TEMPLATE_PARTIAL_SPECIALIZATION)
williamr@2:         // This is the main fare
williamr@2:         
williamr@2:         template<typename T>
williamr@2:         inline T    atanh(const T x)
williamr@2:         {
williamr@2:             using    ::std::abs;
williamr@2:             using    ::std::sqrt;
williamr@2:             using    ::std::log;
williamr@2:             
williamr@2:             using    ::std::numeric_limits;
williamr@2:             
williamr@2:             T const            one = static_cast<T>(1);
williamr@2:             T const            two = static_cast<T>(2);
williamr@2:             
williamr@2:             static T const    taylor_2_bound = sqrt(numeric_limits<T>::epsilon());
williamr@2:             static T const    taylor_n_bound = sqrt(taylor_2_bound);
williamr@2:             
williamr@2:             if        (x < -one)
williamr@2:             {
williamr@2:                 if    (numeric_limits<T>::has_quiet_NaN)
williamr@2:                 {
williamr@2:                     return(numeric_limits<T>::quiet_NaN());
williamr@2:                 }
williamr@2:                 else
williamr@2:                 {
williamr@2:                     ::std::string        error_reporting("Argument to atanh is strictly greater than +1 or strictly smaller than -1!");
williamr@2:                     ::std::domain_error  bad_argument(error_reporting);
williamr@2:                     
williamr@2:                     throw(bad_argument);
williamr@2:                 }
williamr@2:             }
williamr@2:             else if    (x < -one+numeric_limits<T>::epsilon())
williamr@2:             {
williamr@2:                 if    (numeric_limits<T>::has_infinity)
williamr@2:                 {
williamr@2:                     return(-numeric_limits<T>::infinity());
williamr@2:                 }
williamr@2:                 else
williamr@2:                 {
williamr@2:                     ::std::string        error_reporting("Argument to atanh is -1 (result: -Infinity)!");
williamr@2:                     ::std::out_of_range  bad_argument(error_reporting);
williamr@2:                     
williamr@2:                     throw(bad_argument);
williamr@2:                 }
williamr@2:             }
williamr@2:             else if    (x > +one-numeric_limits<T>::epsilon())
williamr@2:             {
williamr@2:                 if    (numeric_limits<T>::has_infinity)
williamr@2:                 {
williamr@2:                     return(+numeric_limits<T>::infinity());
williamr@2:                 }
williamr@2:                 else
williamr@2:                 {
williamr@2:                     ::std::string        error_reporting("Argument to atanh is +1 (result: +Infinity)!");
williamr@2:                     ::std::out_of_range  bad_argument(error_reporting);
williamr@2:                     
williamr@2:                     throw(bad_argument);
williamr@2:                 }
williamr@2:             }
williamr@2:             else if    (x > +one)
williamr@2:             {
williamr@2:                 if    (numeric_limits<T>::has_quiet_NaN)
williamr@2:                 {
williamr@2:                     return(numeric_limits<T>::quiet_NaN());
williamr@2:                 }
williamr@2:                 else
williamr@2:                 {
williamr@2:                     ::std::string        error_reporting("Argument to atanh is strictly greater than +1 or strictly smaller than -1!");
williamr@2:                     ::std::domain_error  bad_argument(error_reporting);
williamr@2:                     
williamr@2:                     throw(bad_argument);
williamr@2:                 }
williamr@2:             }
williamr@2:             else if    (abs(x) >= taylor_n_bound)
williamr@2:             {
williamr@2:                 return(log( (one + x) / (one - x) ) / two);
williamr@2:             }
williamr@2:             else
williamr@2:             {
williamr@2:                 // approximation by taylor series in x at 0 up to order 2
williamr@2:                 T    result = x;
williamr@2:                 
williamr@2:                 if    (abs(x) >= taylor_2_bound)
williamr@2:                 {
williamr@2:                     T    x3 = x*x*x;
williamr@2:                     
williamr@2:                     // approximation by taylor series in x at 0 up to order 4
williamr@2:                     result += x3/static_cast<T>(3);
williamr@2:                 }
williamr@2:                 
williamr@2:                 return(result);
williamr@2:             }
williamr@2:         }
williamr@2: #else
williamr@2:         // These are implementation details (for main fare see below)
williamr@2:         
williamr@2:         namespace detail
williamr@2:         {
williamr@2:             template    <
williamr@2:                             typename T,
williamr@2:                             bool InfinitySupported
williamr@2:                         >
williamr@2:             struct    atanh_helper1_t
williamr@2:             {
williamr@2:                 static T    get_pos_infinity()
williamr@2:                 {
williamr@2:                     return(+::std::numeric_limits<T>::infinity());
williamr@2:                 }
williamr@2:                 
williamr@2:                 static T    get_neg_infinity()
williamr@2:                 {
williamr@2:                     return(-::std::numeric_limits<T>::infinity());
williamr@2:                 }
williamr@2:             };    // boost::math::detail::atanh_helper1_t
williamr@2:             
williamr@2:             
williamr@2:             template<typename T>
williamr@2:             struct    atanh_helper1_t<T, false>
williamr@2:             {
williamr@2:                 static T    get_pos_infinity()
williamr@2:                 {
williamr@2:                     ::std::string        error_reporting("Argument to atanh is +1 (result: +Infinity)!");
williamr@2:                     ::std::out_of_range  bad_argument(error_reporting);
williamr@2:                     
williamr@2:                     throw(bad_argument);
williamr@2:                 }
williamr@2:                 
williamr@2:                 static T    get_neg_infinity()
williamr@2:                 {
williamr@2:                     ::std::string        error_reporting("Argument to atanh is -1 (result: -Infinity)!");
williamr@2:                     ::std::out_of_range  bad_argument(error_reporting);
williamr@2:                     
williamr@2:                     throw(bad_argument);
williamr@2:                 }
williamr@2:             };    // boost::math::detail::atanh_helper1_t
williamr@2:             
williamr@2:             
williamr@2:             template    <
williamr@2:                             typename T,
williamr@2:                             bool QuietNanSupported
williamr@2:                         >
williamr@2:             struct    atanh_helper2_t
williamr@2:             {
williamr@2:                 static T    get_NaN()
williamr@2:                 {
williamr@2:                     return(::std::numeric_limits<T>::quiet_NaN());
williamr@2:                 }
williamr@2:             };    // boost::detail::atanh_helper2_t
williamr@2:             
williamr@2:             
williamr@2:             template<typename T>
williamr@2:             struct    atanh_helper2_t<T, false>
williamr@2:             {
williamr@2:                 static T    get_NaN()
williamr@2:                 {
williamr@2:                     ::std::string        error_reporting("Argument to atanh is strictly greater than +1 or strictly smaller than -1!");
williamr@2:                     ::std::domain_error  bad_argument(error_reporting);
williamr@2:                     
williamr@2:                     throw(bad_argument);
williamr@2:                 }
williamr@2:             };    // boost::detail::atanh_helper2_t
williamr@2:         }    // boost::detail
williamr@2:         
williamr@2:         
williamr@2:         // This is the main fare
williamr@2:         
williamr@2:         template<typename T>
williamr@2:         inline T    atanh(const T x)
williamr@2:         {
williamr@2:             using    ::std::abs;
williamr@2:             using    ::std::sqrt;
williamr@2:             using    ::std::log;
williamr@2:             
williamr@2:             using    ::std::numeric_limits;
williamr@2:             
williamr@2:             typedef  detail::atanh_helper1_t<T, ::std::numeric_limits<T>::has_infinity>    helper1_type;
williamr@2:             typedef  detail::atanh_helper2_t<T, ::std::numeric_limits<T>::has_quiet_NaN>    helper2_type;
williamr@2:             
williamr@2:             
williamr@2:             T const           one = static_cast<T>(1);
williamr@2:             T const           two = static_cast<T>(2);
williamr@2:             
williamr@2:             static T const    taylor_2_bound = sqrt(numeric_limits<T>::epsilon());
williamr@2:             static T const    taylor_n_bound = sqrt(taylor_2_bound);
williamr@2:             
williamr@2:             if        (x < -one)
williamr@2:             {
williamr@2:                 return(helper2_type::get_NaN());
williamr@2:             }
williamr@2:             else if    (x < -one+numeric_limits<T>::epsilon())
williamr@2:             {
williamr@2:                 return(helper1_type::get_neg_infinity());
williamr@2:             }
williamr@2:             else if    (x > +one-numeric_limits<T>::epsilon())
williamr@2:             {
williamr@2:                 return(helper1_type::get_pos_infinity());
williamr@2:             }
williamr@2:             else if    (x > +one)
williamr@2:             {
williamr@2:                 return(helper2_type::get_NaN());
williamr@2:             }
williamr@2:             else if    (abs(x) >= taylor_n_bound)
williamr@2:             {
williamr@2:                 return(log( (one + x) / (one - x) ) / two);
williamr@2:             }
williamr@2:             else
williamr@2:             {
williamr@2:                 // approximation by taylor series in x at 0 up to order 2
williamr@2:                 T    result = x;
williamr@2:                 
williamr@2:                 if    (abs(x) >= taylor_2_bound)
williamr@2:                 {
williamr@2:                     T    x3 = x*x*x;
williamr@2:                     
williamr@2:                     // approximation by taylor series in x at 0 up to order 4
williamr@2:                     result += x3/static_cast<T>(3);
williamr@2:                 }
williamr@2:                 
williamr@2:                 return(result);
williamr@2:             }
williamr@2:         }
williamr@2: #endif /* defined(BOOST_NO_TEMPLATE_PARTIAL_SPECIALIZATION) */
williamr@2:     }
williamr@2: }
williamr@2: 
williamr@2: #endif /* BOOST_ATANH_HPP */
williamr@2: