//    boost atanh.hpp header file

 //  (C) Copyright Hubert Holin 2001.

 //  Distributed under the Boost Software License, Version 1.0. (See

 //  accompanying file LICENSE_1_0.txt or copy at

 //  http://www.boost.org/LICENSE_1_0.txt)

 // See http://www.boost.org for updates, documentation, and revision history.

 #ifndef BOOST_ATANH_HPP

 #define BOOST_ATANH_HPP

 #include <cmath>

 #include <limits>

 #include <string>

 #include <stdexcept>

 #include <boost/config.hpp>

 // This is the inverse of the hyperbolic tangent function.

 namespace boost

 {

     namespace math

     {

 #if defined(__GNUC__) && (__GNUC__ < 3)

         // gcc 2.x ignores function scope using declarations,

         // put them in the scope of the enclosing namespace instead:

         using    ::std::abs;

         using    ::std::sqrt;

         using    ::std::log;

         using    ::std::numeric_limits;

 #endif

 #if defined(BOOST_NO_TEMPLATE_PARTIAL_SPECIALIZATION)

         // This is the main fare

         template<typename T>

         inline T    atanh(const T x)

         {

             using    ::std::abs;

             using    ::std::sqrt;

             using    ::std::log;

             using    ::std::numeric_limits;

             T const            one = static_cast<T>(1);

             T const            two = static_cast<T>(2);

             static T const    taylor_2_bound = sqrt(numeric_limits<T>::epsilon());

             static T const    taylor_n_bound = sqrt(taylor_2_bound);

             if        (x < -one)

             {

                 if    (numeric_limits<T>::has_quiet_NaN)

                 {

                     return(numeric_limits<T>::quiet_NaN());

                 }

                 else

                 {

                     ::std::string        error_reporting("Argument to atanh is strictly greater than +1 or strictly smaller than -1!");

                     ::std::domain_error  bad_argument(error_reporting);

                     throw(bad_argument);

                 }

             }

             else if    (x < -one+numeric_limits<T>::epsilon())

             {

                 if    (numeric_limits<T>::has_infinity)

                 {

                     return(-numeric_limits<T>::infinity());

                 }

                 else

                 {

                     ::std::string        error_reporting("Argument to atanh is -1 (result: -Infinity)!");

                     ::std::out_of_range  bad_argument(error_reporting);

                     throw(bad_argument);

                 }

             }

             else if    (x > +one-numeric_limits<T>::epsilon())

             {

                 if    (numeric_limits<T>::has_infinity)

                 {

                     return(+numeric_limits<T>::infinity());

                 }

                 else

                 {

                     ::std::string        error_reporting("Argument to atanh is +1 (result: +Infinity)!");

                     ::std::out_of_range  bad_argument(error_reporting);

                     throw(bad_argument);

                 }

             }

             else if    (x > +one)

             {

                 if    (numeric_limits<T>::has_quiet_NaN)

                 {

                     return(numeric_limits<T>::quiet_NaN());

                 }

                 else

                 {

                     ::std::string        error_reporting("Argument to atanh is strictly greater than +1 or strictly smaller than -1!");

                     ::std::domain_error  bad_argument(error_reporting);

                     throw(bad_argument);

                 }

             }

             else if    (abs(x) >= taylor_n_bound)

             {

                 return(log( (one + x) / (one - x) ) / two);

             }

             else

             {

                 // approximation by taylor series in x at 0 up to order 2

                 T    result = x;

                 if    (abs(x) >= taylor_2_bound)

                 {

                     T    x3 = x*x*x;

                     // approximation by taylor series in x at 0 up to order 4

                     result += x3/static_cast<T>(3);

                 }

                 return(result);

             }

         }

 #else

         // These are implementation details (for main fare see below)

         namespace detail

         {

             template    <

                             typename T,

                             bool InfinitySupported

                         >

             struct    atanh_helper1_t

             {

                 static T    get_pos_infinity()

                 {

                     return(+::std::numeric_limits<T>::infinity());

                 }

                 static T    get_neg_infinity()

                 {

                     return(-::std::numeric_limits<T>::infinity());

                 }

             };    // boost::math::detail::atanh_helper1_t

             template<typename T>

             struct    atanh_helper1_t<T, false>

             {

                 static T    get_pos_infinity()

                 {

                     ::std::string        error_reporting("Argument to atanh is +1 (result: +Infinity)!");

                     ::std::out_of_range  bad_argument(error_reporting);

                     throw(bad_argument);

                 }

                 static T    get_neg_infinity()

                 {

                     ::std::string        error_reporting("Argument to atanh is -1 (result: -Infinity)!");

                     ::std::out_of_range  bad_argument(error_reporting);

                     throw(bad_argument);

                 }

             };    // boost::math::detail::atanh_helper1_t

             template    <

                             typename T,

                             bool QuietNanSupported

                         >

             struct    atanh_helper2_t

             {

                 static T    get_NaN()

                 {

                     return(::std::numeric_limits<T>::quiet_NaN());

                 }

             };    // boost::detail::atanh_helper2_t

             template<typename T>

             struct    atanh_helper2_t<T, false>

             {

                 static T    get_NaN()

                 {

                     ::std::string        error_reporting("Argument to atanh is strictly greater than +1 or strictly smaller than -1!");

                     ::std::domain_error  bad_argument(error_reporting);

                     throw(bad_argument);

                 }

             };    // boost::detail::atanh_helper2_t

         }    // boost::detail

         // This is the main fare

         template<typename T>

         inline T    atanh(const T x)

         {

             using    ::std::abs;

             using    ::std::sqrt;

             using    ::std::log;

             using    ::std::numeric_limits;

             typedef  detail::atanh_helper1_t<T, ::std::numeric_limits<T>::has_infinity>    helper1_type;

             typedef  detail::atanh_helper2_t<T, ::std::numeric_limits<T>::has_quiet_NaN>    helper2_type;

             T const           one = static_cast<T>(1);

             T const           two = static_cast<T>(2);

             static T const    taylor_2_bound = sqrt(numeric_limits<T>::epsilon());

             static T const    taylor_n_bound = sqrt(taylor_2_bound);

             if        (x < -one)

             {

                 return(helper2_type::get_NaN());

             }

             else if    (x < -one+numeric_limits<T>::epsilon())

             {

                 return(helper1_type::get_neg_infinity());

             }

             else if    (x > +one-numeric_limits<T>::epsilon())

             {

                 return(helper1_type::get_pos_infinity());

             }

             else if    (x > +one)

             {

                 return(helper2_type::get_NaN());

             }

             else if    (abs(x) >= taylor_n_bound)

             {

                 return(log( (one + x) / (one - x) ) / two);

             }

             else

             {

                 // approximation by taylor series in x at 0 up to order 2

                 T    result = x;

                 if    (abs(x) >= taylor_2_bound)

                 {

                     T    x3 = x*x*x;

                     // approximation by taylor series in x at 0 up to order 4

                     result += x3/static_cast<T>(3);

                 }

                 return(result);

             }

         }

 #endif /* defined(BOOST_NO_TEMPLATE_PARTIAL_SPECIALIZATION) */

     }

 }

 #endif /* BOOST_ATANH_HPP */