williamr@2: // boost asinh.hpp header file williamr@2: williamr@2: // (C) Copyright Eric Ford & 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_ASINH_HPP williamr@2: #define BOOST_ASINH_HPP williamr@2: williamr@2: williamr@2: #include williamr@2: #include williamr@2: #include williamr@2: #include williamr@2: williamr@2: williamr@2: #include williamr@2: williamr@2: williamr@2: // This is the inverse of the hyperbolic sine 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: template williamr@2: inline T asinh(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: williamr@2: T const one = static_cast(1); williamr@2: T const two = static_cast(2); williamr@2: williamr@2: static T const taylor_2_bound = sqrt(numeric_limits::epsilon()); williamr@2: static T const taylor_n_bound = sqrt(taylor_2_bound); williamr@2: static T const upper_taylor_2_bound = one/taylor_2_bound; williamr@2: static T const upper_taylor_n_bound = one/taylor_n_bound; williamr@2: williamr@2: if (x >= +taylor_n_bound) williamr@2: { williamr@2: if (x > upper_taylor_n_bound) williamr@2: { williamr@2: if (x > upper_taylor_2_bound) williamr@2: { williamr@2: // approximation by laurent series in 1/x at 0+ order from -1 to 0 williamr@2: return( log( x * two) ); williamr@2: } williamr@2: else williamr@2: { williamr@2: // approximation by laurent series in 1/x at 0+ order from -1 to 1 williamr@2: return( log( x*two + (one/(x*two)) ) ); williamr@2: } williamr@2: } williamr@2: else williamr@2: { williamr@2: return( log( x + sqrt(x*x+one) ) ); williamr@2: } williamr@2: } williamr@2: else if (x <= -taylor_n_bound) williamr@2: { williamr@2: return(-asinh(-x)); 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(6); williamr@2: } williamr@2: williamr@2: return(result); williamr@2: } williamr@2: } williamr@2: } williamr@2: } williamr@2: williamr@2: #endif /* BOOST_ASINH_HPP */