//    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 */