Chris@16: //  (C) Copyright John Maddock 2006.
Chris@16: //  Use, modification and distribution are subject to the
Chris@16: //  Boost Software License, Version 1.0. (See accompanying file
Chris@16: //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
Chris@16: 
Chris@16: #ifndef BOOST_MATH_EXPM1_INCLUDED
Chris@16: #define BOOST_MATH_EXPM1_INCLUDED
Chris@16: 
Chris@16: #ifdef _MSC_VER
Chris@16: #pragma once
Chris@16: #endif
Chris@16: 
Chris@16: #include <boost/config/no_tr1/cmath.hpp>
Chris@16: #include <math.h> // platform's ::expm1
Chris@16: #include <boost/limits.hpp>
Chris@16: #include <boost/math/tools/config.hpp>
Chris@16: #include <boost/math/tools/series.hpp>
Chris@16: #include <boost/math/tools/precision.hpp>
Chris@16: #include <boost/math/tools/big_constant.hpp>
Chris@16: #include <boost/math/policies/error_handling.hpp>
Chris@16: #include <boost/math/tools/rational.hpp>
Chris@16: #include <boost/math/special_functions/math_fwd.hpp>
Chris@16: #include <boost/mpl/less_equal.hpp>
Chris@16: 
Chris@16: #ifndef BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS
Chris@16: #  include <boost/static_assert.hpp>
Chris@16: #else
Chris@16: #  include <boost/assert.hpp>
Chris@16: #endif
Chris@16: 
Chris@16: namespace boost{ namespace math{
Chris@16: 
Chris@16: namespace detail
Chris@16: {
Chris@16:   // Functor expm1_series returns the next term in the Taylor series
Chris@16:   // x^k / k!
Chris@16:   // each time that operator() is invoked.
Chris@16:   //
Chris@16:   template <class T>
Chris@16:   struct expm1_series
Chris@16:   {
Chris@16:      typedef T result_type;
Chris@16: 
Chris@16:      expm1_series(T x)
Chris@16:         : k(0), m_x(x), m_term(1) {}
Chris@16: 
Chris@16:      T operator()()
Chris@16:      {
Chris@16:         ++k;
Chris@16:         m_term *= m_x;
Chris@16:         m_term /= k;
Chris@16:         return m_term;
Chris@16:      }
Chris@16: 
Chris@16:      int count()const
Chris@16:      {
Chris@16:         return k;
Chris@16:      }
Chris@16: 
Chris@16:   private:
Chris@16:      int k;
Chris@16:      const T m_x;
Chris@16:      T m_term;
Chris@16:      expm1_series(const expm1_series&);
Chris@16:      expm1_series& operator=(const expm1_series&);
Chris@16:   };
Chris@16: 
Chris@16: template <class T, class Policy, class tag>
Chris@16: struct expm1_initializer
Chris@16: {
Chris@16:    struct init
Chris@16:    {
Chris@16:       init()
Chris@16:       {
Chris@16:          do_init(tag());
Chris@16:       }
Chris@16:       template <int N>
Chris@16:       static void do_init(const mpl::int_<N>&){}
Chris@16:       static void do_init(const mpl::int_<64>&)
Chris@16:       {
Chris@16:          expm1(T(0.5));
Chris@16:       }
Chris@16:       static void do_init(const mpl::int_<113>&)
Chris@16:       {
Chris@16:          expm1(T(0.5));
Chris@16:       }
Chris@16:       void force_instantiate()const{}
Chris@16:    };
Chris@16:    static const init initializer;
Chris@16:    static void force_instantiate()
Chris@16:    {
Chris@16:       initializer.force_instantiate();
Chris@16:    }
Chris@16: };
Chris@16: 
Chris@16: template <class T, class Policy, class tag>
Chris@16: const typename expm1_initializer<T, Policy, tag>::init expm1_initializer<T, Policy, tag>::initializer;
Chris@16: 
Chris@16: //
Chris@16: // Algorithm expm1 is part of C99, but is not yet provided by many compilers.
Chris@16: //
Chris@16: // This version uses a Taylor series expansion for 0.5 > |x| > epsilon.
Chris@16: //
Chris@16: template <class T, class Policy>
Chris@16: T expm1_imp(T x, const mpl::int_<0>&, const Policy& pol)
Chris@16: {
Chris@16:    BOOST_MATH_STD_USING
Chris@16: 
Chris@16:    T a = fabs(x);
Chris@16:    if(a > T(0.5f))
Chris@16:    {
Chris@16:       if(a >= tools::log_max_value<T>())
Chris@16:       {
Chris@16:          if(x > 0)
Chris@16:             return policies::raise_overflow_error<T>("boost::math::expm1<%1%>(%1%)", 0, pol);
Chris@16:          return -1;
Chris@16:       }
Chris@16:       return exp(x) - T(1);
Chris@16:    }
Chris@16:    if(a < tools::epsilon<T>())
Chris@16:       return x;
Chris@16:    detail::expm1_series<T> s(x);
Chris@16:    boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
Chris@16: #if !BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x582)) && !BOOST_WORKAROUND(__EDG_VERSION__, <= 245)
Chris@16:    T result = tools::sum_series(s, policies::get_epsilon<T, Policy>(), max_iter);
Chris@16: #else
Chris@16:    T zero = 0;
Chris@16:    T result = tools::sum_series(s, policies::get_epsilon<T, Policy>(), max_iter, zero);
Chris@16: #endif
Chris@16:    policies::check_series_iterations<T>("boost::math::expm1<%1%>(%1%)", max_iter, pol);
Chris@16:    return result;
Chris@16: }
Chris@16: 
Chris@16: template <class T, class P>
Chris@16: T expm1_imp(T x, const mpl::int_<53>&, const P& pol)
Chris@16: {
Chris@16:    BOOST_MATH_STD_USING
Chris@16: 
Chris@16:    T a = fabs(x);
Chris@16:    if(a > T(0.5L))
Chris@16:    {
Chris@16:       if(a >= tools::log_max_value<T>())
Chris@16:       {
Chris@16:          if(x > 0)
Chris@16:             return policies::raise_overflow_error<T>("boost::math::expm1<%1%>(%1%)", 0, pol);
Chris@16:          return -1;
Chris@16:       }
Chris@16:       return exp(x) - T(1);
Chris@16:    }
Chris@16:    if(a < tools::epsilon<T>())
Chris@16:       return x;
Chris@16: 
Chris@16:    static const float Y = 0.10281276702880859e1f;
Chris@16:    static const T n[] = { static_cast<T>(-0.28127670288085937e-1), static_cast<T>(0.51278186299064534e0), static_cast<T>(-0.6310029069350198e-1), static_cast<T>(0.11638457975729296e-1), static_cast<T>(-0.52143390687521003e-3), static_cast<T>(0.21491399776965688e-4) };
Chris@16:    static const T d[] = { 1, static_cast<T>(-0.45442309511354755e0), static_cast<T>(0.90850389570911714e-1), static_cast<T>(-0.10088963629815502e-1), static_cast<T>(0.63003407478692265e-3), static_cast<T>(-0.17976570003654402e-4) };
Chris@16: 
Chris@16:    T result = x * Y + x * tools::evaluate_polynomial(n, x) / tools::evaluate_polynomial(d, x);
Chris@16:    return result;
Chris@16: }
Chris@16: 
Chris@16: template <class T, class P>
Chris@16: T expm1_imp(T x, const mpl::int_<64>&, const P& pol)
Chris@16: {
Chris@16:    BOOST_MATH_STD_USING
Chris@16: 
Chris@16:    T a = fabs(x);
Chris@16:    if(a > T(0.5L))
Chris@16:    {
Chris@16:       if(a >= tools::log_max_value<T>())
Chris@16:       {
Chris@16:          if(x > 0)
Chris@16:             return policies::raise_overflow_error<T>("boost::math::expm1<%1%>(%1%)", 0, pol);
Chris@16:          return -1;
Chris@16:       }
Chris@16:       return exp(x) - T(1);
Chris@16:    }
Chris@16:    if(a < tools::epsilon<T>())
Chris@16:       return x;
Chris@16: 
Chris@16:    static const float Y = 0.10281276702880859375e1f;
Chris@16:    static const T n[] = { 
Chris@16:       BOOST_MATH_BIG_CONSTANT(T, 64, -0.281276702880859375e-1), 
Chris@16:        BOOST_MATH_BIG_CONSTANT(T, 64, 0.512980290285154286358e0), 
Chris@16:        BOOST_MATH_BIG_CONSTANT(T, 64, -0.667758794592881019644e-1),
Chris@16:        BOOST_MATH_BIG_CONSTANT(T, 64, 0.131432469658444745835e-1),
Chris@16:        BOOST_MATH_BIG_CONSTANT(T, 64, -0.72303795326880286965e-3),
Chris@16:        BOOST_MATH_BIG_CONSTANT(T, 64, 0.447441185192951335042e-4),
Chris@16:        BOOST_MATH_BIG_CONSTANT(T, 64, -0.714539134024984593011e-6)
Chris@16:    };
Chris@16:    static const T d[] = { 
Chris@16:       BOOST_MATH_BIG_CONSTANT(T, 64, 1.0),
Chris@16:       BOOST_MATH_BIG_CONSTANT(T, 64, -0.461477618025562520389e0),
Chris@16:       BOOST_MATH_BIG_CONSTANT(T, 64, 0.961237488025708540713e-1),
Chris@16:       BOOST_MATH_BIG_CONSTANT(T, 64, -0.116483957658204450739e-1),
Chris@16:       BOOST_MATH_BIG_CONSTANT(T, 64, 0.873308008461557544458e-3),
Chris@16:       BOOST_MATH_BIG_CONSTANT(T, 64, -0.387922804997682392562e-4),
Chris@16:       BOOST_MATH_BIG_CONSTANT(T, 64, 0.807473180049193557294e-6)
Chris@16:    };
Chris@16: 
Chris@16:    T result = x * Y + x * tools::evaluate_polynomial(n, x) / tools::evaluate_polynomial(d, x);
Chris@16:    return result;
Chris@16: }
Chris@16: 
Chris@16: template <class T, class P>
Chris@16: T expm1_imp(T x, const mpl::int_<113>&, const P& pol)
Chris@16: {
Chris@16:    BOOST_MATH_STD_USING
Chris@16: 
Chris@16:    T a = fabs(x);
Chris@16:    if(a > T(0.5L))
Chris@16:    {
Chris@16:       if(a >= tools::log_max_value<T>())
Chris@16:       {
Chris@16:          if(x > 0)
Chris@16:             return policies::raise_overflow_error<T>("boost::math::expm1<%1%>(%1%)", 0, pol);
Chris@16:          return -1;
Chris@16:       }
Chris@16:       return exp(x) - T(1);
Chris@16:    }
Chris@16:    if(a < tools::epsilon<T>())
Chris@16:       return x;
Chris@16: 
Chris@16:    static const float Y = 0.10281276702880859375e1f;
Chris@16:    static const T n[] = { 
Chris@16:       BOOST_MATH_BIG_CONSTANT(T, 113, -0.28127670288085937499999999999999999854e-1),
Chris@16:       BOOST_MATH_BIG_CONSTANT(T, 113, 0.51278156911210477556524452177540792214e0),
Chris@16:       BOOST_MATH_BIG_CONSTANT(T, 113, -0.63263178520747096729500254678819588223e-1),
Chris@16:       BOOST_MATH_BIG_CONSTANT(T, 113, 0.14703285606874250425508446801230572252e-1),
Chris@16:       BOOST_MATH_BIG_CONSTANT(T, 113, -0.8675686051689527802425310407898459386e-3),
Chris@16:       BOOST_MATH_BIG_CONSTANT(T, 113, 0.88126359618291165384647080266133492399e-4),
Chris@16:       BOOST_MATH_BIG_CONSTANT(T, 113, -0.25963087867706310844432390015463138953e-5),
Chris@16:       BOOST_MATH_BIG_CONSTANT(T, 113, 0.14226691087800461778631773363204081194e-6),
Chris@16:       BOOST_MATH_BIG_CONSTANT(T, 113, -0.15995603306536496772374181066765665596e-8),
Chris@16:       BOOST_MATH_BIG_CONSTANT(T, 113, 0.45261820069007790520447958280473183582e-10)
Chris@16:    };
Chris@16:    static const T d[] = { 
Chris@16:       BOOST_MATH_BIG_CONSTANT(T, 113, 1.0),
Chris@16:       BOOST_MATH_BIG_CONSTANT(T, 113, -0.45441264709074310514348137469214538853e0),
Chris@16:       BOOST_MATH_BIG_CONSTANT(T, 113, 0.96827131936192217313133611655555298106e-1),
Chris@16:       BOOST_MATH_BIG_CONSTANT(T, 113, -0.12745248725908178612540554584374876219e-1),
Chris@16:       BOOST_MATH_BIG_CONSTANT(T, 113, 0.11473613871583259821612766907781095472e-2),
Chris@16:       BOOST_MATH_BIG_CONSTANT(T, 113, -0.73704168477258911962046591907690764416e-4),
Chris@16:       BOOST_MATH_BIG_CONSTANT(T, 113, 0.34087499397791555759285503797256103259e-5),
Chris@16:       BOOST_MATH_BIG_CONSTANT(T, 113, -0.11114024704296196166272091230695179724e-6),
Chris@16:       BOOST_MATH_BIG_CONSTANT(T, 113, 0.23987051614110848595909588343223896577e-8),
Chris@16:       BOOST_MATH_BIG_CONSTANT(T, 113, -0.29477341859111589208776402638429026517e-10),
Chris@16:       BOOST_MATH_BIG_CONSTANT(T, 113, 0.13222065991022301420255904060628100924e-12)
Chris@16:    };
Chris@16: 
Chris@16:    T result = x * Y + x * tools::evaluate_polynomial(n, x) / tools::evaluate_polynomial(d, x);
Chris@16:    return result;
Chris@16: }
Chris@16: 
Chris@16: } // namespace detail
Chris@16: 
Chris@16: template <class T, class Policy>
Chris@16: inline typename tools::promote_args<T>::type expm1(T x, const Policy& /* pol */)
Chris@16: {
Chris@16:    typedef typename tools::promote_args<T>::type result_type;
Chris@16:    typedef typename policies::evaluation<result_type, Policy>::type value_type;
Chris@16:    typedef typename policies::precision<result_type, Policy>::type precision_type;
Chris@16:    typedef typename policies::normalise<
Chris@16:       Policy, 
Chris@16:       policies::promote_float<false>, 
Chris@16:       policies::promote_double<false>, 
Chris@16:       policies::discrete_quantile<>,
Chris@16:       policies::assert_undefined<> >::type forwarding_policy;
Chris@16: 
Chris@16:    typedef typename mpl::if_c<
Chris@16:       ::std::numeric_limits<result_type>::is_specialized == 0,
Chris@16:       mpl::int_<0>,  // no numeric_limits, use generic solution
Chris@16:       typename mpl::if_<
Chris@16:          typename mpl::less_equal<precision_type, mpl::int_<53> >::type,
Chris@16:          mpl::int_<53>,  // double
Chris@16:          typename mpl::if_<
Chris@16:             typename mpl::less_equal<precision_type, mpl::int_<64> >::type,
Chris@16:             mpl::int_<64>, // 80-bit long double
Chris@16:             typename mpl::if_<
Chris@16:                typename mpl::less_equal<precision_type, mpl::int_<113> >::type,
Chris@16:                mpl::int_<113>, // 128-bit long double
Chris@16:                mpl::int_<0> // too many bits, use generic version.
Chris@16:             >::type
Chris@16:          >::type
Chris@16:       >::type
Chris@16:    >::type tag_type;
Chris@16: 
Chris@16:    detail::expm1_initializer<value_type, forwarding_policy, tag_type>::force_instantiate();
Chris@16:    
Chris@16:    return policies::checked_narrowing_cast<result_type, forwarding_policy>(detail::expm1_imp(
Chris@16:       static_cast<value_type>(x),
Chris@16:       tag_type(), forwarding_policy()), "boost::math::expm1<%1%>(%1%)");
Chris@16: }
Chris@16: 
Chris@16: #ifdef expm1
Chris@16: #  ifndef BOOST_HAS_expm1
Chris@16: #     define BOOST_HAS_expm1
Chris@16: #  endif
Chris@16: #  undef expm1
Chris@16: #endif
Chris@16: 
Chris@16: #if defined(BOOST_HAS_EXPM1) && !(defined(__osf__) && defined(__DECCXX_VER))
Chris@16: #  ifdef BOOST_MATH_USE_C99
Chris@16: inline float expm1(float x, const policies::policy<>&){ return ::expm1f(x); }
Chris@16: #     ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
Chris@16: inline long double expm1(long double x, const policies::policy<>&){ return ::expm1l(x); }
Chris@16: #     endif
Chris@16: #  else
Chris@16: inline float expm1(float x, const policies::policy<>&){ return static_cast<float>(::expm1(x)); }
Chris@16: #  endif
Chris@16: inline double expm1(double x, const policies::policy<>&){ return ::expm1(x); }
Chris@16: #endif
Chris@16: 
Chris@16: template <class T>
Chris@16: inline typename tools::promote_args<T>::type expm1(T x)
Chris@16: {
Chris@16:    return expm1(x, policies::policy<>());
Chris@16: }
Chris@16: 
Chris@16: #if BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x564))
Chris@16: inline float expm1(float z)
Chris@16: {
Chris@16:    return expm1<float>(z);
Chris@16: }
Chris@16: inline double expm1(double z)
Chris@16: {
Chris@16:    return expm1<double>(z);
Chris@16: }
Chris@16: #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
Chris@16: inline long double expm1(long double z)
Chris@16: {
Chris@16:    return expm1<long double>(z);
Chris@16: }
Chris@16: #endif
Chris@16: #endif
Chris@16: 
Chris@16: } // namespace math
Chris@16: } // namespace boost
Chris@16: 
Chris@16: #endif // BOOST_MATH_HYPOT_INCLUDED
Chris@16: 
Chris@16: 
Chris@16: 
Chris@16: