math/special_functions/expm1.hpp

0001 //  (C) Copyright John Maddock 2006.
0002 //  Use, modification and distribution are subject to the
0003 //  Boost Software License, Version 1.0. (See accompanying file
0004 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
0005
0006 #ifndef BOOST_MATH_EXPM1_INCLUDED
0007 #define BOOST_MATH_EXPM1_INCLUDED
0008
0009 #ifdef _MSC_VER
0010 #pragma once
0011 #endif
0012
0013 #include <cmath>
0014 #include <cstdint>
0015 #include <limits>
0016 #include <boost/math/tools/config.hpp>
0017 #include <boost/math/tools/series.hpp>
0018 #include <boost/math/tools/precision.hpp>
0019 #include <boost/math/tools/big_constant.hpp>
0020 #include <boost/math/policies/error_handling.hpp>
0021 #include <boost/math/tools/rational.hpp>
0022 #include <boost/math/special_functions/math_fwd.hpp>
0023 #include <boost/math/tools/assert.hpp>
0024
0025 #if defined(__GNUC__) && defined(BOOST_MATH_USE_FLOAT128)
0026 //
0027 // This is the only way we can avoid
0028 // warning: non-standard suffix on floating constant [-Wpedantic]
0029 // when building with -Wall -pedantic.  Neither __extension__
0030 // nor #pragma diagnostic ignored work :(
0031 //
0032 #pragma GCC system_header
0033 #endif
0034
0035 namespace boost{ namespace math{
0036
0037 namespace detail
0038 {
0039   // Functor expm1_series returns the next term in the Taylor series
0040   // x^k / k!
0041   // each time that operator() is invoked.
0042   //
0043   template <class T>
0044   struct expm1_series
0045   {
0046      typedef T result_type;
0047
0048      expm1_series(T x)
0049         : k(0), m_x(x), m_term(1) {}
0050
0051      T operator()()
0052      {
0053         ++k;
0054         m_term *= m_x;
0055         m_term /= k;
0056         return m_term;
0057      }
0058
0059      int count()const
0060      {
0061         return k;
0062      }
0063
0064   private:
0065      int k;
0066      const T m_x;
0067      T m_term;
0068      expm1_series(const expm1_series&) = delete;
0069      expm1_series& operator=(const expm1_series&) = delete;
0070   };
0071
0072 template <class T, class Policy, class tag>
0073 struct expm1_initializer
0074 {
0075    struct init
0076    {
0077       init()
0078       {
0079          do_init(tag());
0080       }
0081       template <int N>
0082       static void do_init(const std::integral_constant<int, N>&){}
0083       static void do_init(const std::integral_constant<int, 64>&)
0084       {
0085          expm1(T(0.5));
0086       }
0087       static void do_init(const std::integral_constant<int, 113>&)
0088       {
0089          expm1(T(0.5));
0090       }
0091       void force_instantiate()const{}
0092    };
0093    static const init initializer;
0094    static void force_instantiate()
0095    {
0096       initializer.force_instantiate();
0097    }
0098 };
0099
0100 template <class T, class Policy, class tag>
0101 const typename expm1_initializer<T, Policy, tag>::init expm1_initializer<T, Policy, tag>::initializer;
0102
0103 //
0104 // Algorithm expm1 is part of C99, but is not yet provided by many compilers.
0105 //
0106 // This version uses a Taylor series expansion for 0.5 > |x| > epsilon.
0107 //
0108 template <class T, class Policy>
0109 T expm1_imp(T x, const std::integral_constant<int, 0>&, const Policy& pol)
0110 {
0111    BOOST_MATH_STD_USING
0112
0113    T a = fabs(x);
0114    if((boost::math::isnan)(a))
0115    {
0116       return policies::raise_domain_error<T>("boost::math::expm1<%1%>(%1%)", "expm1 requires a finite argument, but got %1%", a, pol);
0117    }
0118    if(a > T(0.5f))
0119    {
0120       if(a >= tools::log_max_value<T>())
0121       {
0122          if(x > 0)
0123             return policies::raise_overflow_error<T>("boost::math::expm1<%1%>(%1%)", nullptr, pol);
0124          return -1;
0125       }
0126       return exp(x) - T(1);
0127    }
0128    if(a < tools::epsilon<T>())
0129       return x;
0130    detail::expm1_series<T> s(x);
0131    std::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
0132
0133    T result = tools::sum_series(s, policies::get_epsilon<T, Policy>(), max_iter);
0134
0135    policies::check_series_iterations<T>("boost::math::expm1<%1%>(%1%)", max_iter, pol);
0136    return result;
0137 }
0138
0139 template <class T, class P>
0140 T expm1_imp(T x, const std::integral_constant<int, 53>&, const P& pol)
0141 {
0142    BOOST_MATH_STD_USING
0143
0144    T a = fabs(x);
0145    if(a > T(0.5L))
0146    {
0147       if(a >= tools::log_max_value<T>())
0148       {
0149          if(x > 0)
0150             return policies::raise_overflow_error<T>("boost::math::expm1<%1%>(%1%)", nullptr, pol);
0151          return -1;
0152       }
0153       return exp(x) - T(1);
0154    }
0155    if(a < tools::epsilon<T>())
0156       return x;
0157
0158    static const float Y = 0.10281276702880859e1f;
0159    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) };
0160    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) };
0161
0162    T result = x * Y + x * tools::evaluate_polynomial(n, x) / tools::evaluate_polynomial(d, x);
0163    return result;
0164 }
0165
0166 template <class T, class P>
0167 T expm1_imp(T x, const std::integral_constant<int, 64>&, const P& pol)
0168 {
0169    BOOST_MATH_STD_USING
0170
0171    T a = fabs(x);
0172    if(a > T(0.5L))
0173    {
0174       if(a >= tools::log_max_value<T>())
0175       {
0176          if(x > 0)
0177             return policies::raise_overflow_error<T>("boost::math::expm1<%1%>(%1%)", nullptr, pol);
0178          return -1;
0179       }
0180       return exp(x) - T(1);
0181    }
0182    if(a < tools::epsilon<T>())
0183       return x;
0184
0185    static const float Y = 0.10281276702880859375e1f;
0186    static const T n[] = {
0187       BOOST_MATH_BIG_CONSTANT(T, 64, -0.281276702880859375e-1),
0188        BOOST_MATH_BIG_CONSTANT(T, 64, 0.512980290285154286358e0),
0189        BOOST_MATH_BIG_CONSTANT(T, 64, -0.667758794592881019644e-1),
0190        BOOST_MATH_BIG_CONSTANT(T, 64, 0.131432469658444745835e-1),
0191        BOOST_MATH_BIG_CONSTANT(T, 64, -0.72303795326880286965e-3),
0192        BOOST_MATH_BIG_CONSTANT(T, 64, 0.447441185192951335042e-4),
0193        BOOST_MATH_BIG_CONSTANT(T, 64, -0.714539134024984593011e-6)
0194    };
0195    static const T d[] = {
0196       BOOST_MATH_BIG_CONSTANT(T, 64, 1.0),
0197       BOOST_MATH_BIG_CONSTANT(T, 64, -0.461477618025562520389e0),
0198       BOOST_MATH_BIG_CONSTANT(T, 64, 0.961237488025708540713e-1),
0199       BOOST_MATH_BIG_CONSTANT(T, 64, -0.116483957658204450739e-1),
0200       BOOST_MATH_BIG_CONSTANT(T, 64, 0.873308008461557544458e-3),
0201       BOOST_MATH_BIG_CONSTANT(T, 64, -0.387922804997682392562e-4),
0202       BOOST_MATH_BIG_CONSTANT(T, 64, 0.807473180049193557294e-6)
0203    };
0204
0205    T result = x * Y + x * tools::evaluate_polynomial(n, x) / tools::evaluate_polynomial(d, x);
0206    return result;
0207 }
0208
0209 template <class T, class P>
0210 T expm1_imp(T x, const std::integral_constant<int, 113>&, const P& pol)
0211 {
0212    BOOST_MATH_STD_USING
0213
0214    T a = fabs(x);
0215    if(a > T(0.5L))
0216    {
0217       if(a >= tools::log_max_value<T>())
0218       {
0219          if(x > 0)
0220             return policies::raise_overflow_error<T>("boost::math::expm1<%1%>(%1%)", nullptr, pol);
0221          return -1;
0222       }
0223       return exp(x) - T(1);
0224    }
0225    if(a < tools::epsilon<T>())
0226       return x;
0227
0228    static const float Y = 0.10281276702880859375e1f;
0229    static const T n[] = {
0230       BOOST_MATH_BIG_CONSTANT(T, 113, -0.28127670288085937499999999999999999854e-1),
0231       BOOST_MATH_BIG_CONSTANT(T, 113, 0.51278156911210477556524452177540792214e0),
0232       BOOST_MATH_BIG_CONSTANT(T, 113, -0.63263178520747096729500254678819588223e-1),
0233       BOOST_MATH_BIG_CONSTANT(T, 113, 0.14703285606874250425508446801230572252e-1),
0234       BOOST_MATH_BIG_CONSTANT(T, 113, -0.8675686051689527802425310407898459386e-3),
0235       BOOST_MATH_BIG_CONSTANT(T, 113, 0.88126359618291165384647080266133492399e-4),
0236       BOOST_MATH_BIG_CONSTANT(T, 113, -0.25963087867706310844432390015463138953e-5),
0237       BOOST_MATH_BIG_CONSTANT(T, 113, 0.14226691087800461778631773363204081194e-6),
0238       BOOST_MATH_BIG_CONSTANT(T, 113, -0.15995603306536496772374181066765665596e-8),
0239       BOOST_MATH_BIG_CONSTANT(T, 113, 0.45261820069007790520447958280473183582e-10)
0240    };
0241    static const T d[] = {
0242       BOOST_MATH_BIG_CONSTANT(T, 113, 1.0),
0243       BOOST_MATH_BIG_CONSTANT(T, 113, -0.45441264709074310514348137469214538853e0),
0244       BOOST_MATH_BIG_CONSTANT(T, 113, 0.96827131936192217313133611655555298106e-1),
0245       BOOST_MATH_BIG_CONSTANT(T, 113, -0.12745248725908178612540554584374876219e-1),
0246       BOOST_MATH_BIG_CONSTANT(T, 113, 0.11473613871583259821612766907781095472e-2),
0247       BOOST_MATH_BIG_CONSTANT(T, 113, -0.73704168477258911962046591907690764416e-4),
0248       BOOST_MATH_BIG_CONSTANT(T, 113, 0.34087499397791555759285503797256103259e-5),
0249       BOOST_MATH_BIG_CONSTANT(T, 113, -0.11114024704296196166272091230695179724e-6),
0250       BOOST_MATH_BIG_CONSTANT(T, 113, 0.23987051614110848595909588343223896577e-8),
0251       BOOST_MATH_BIG_CONSTANT(T, 113, -0.29477341859111589208776402638429026517e-10),
0252       BOOST_MATH_BIG_CONSTANT(T, 113, 0.13222065991022301420255904060628100924e-12)
0253    };
0254
0255    T result = x * Y + x * tools::evaluate_polynomial(n, x) / tools::evaluate_polynomial(d, x);
0256    return result;
0257 }
0258
0259 } // namespace detail
0260
0261 template <class T, class Policy>
0262 inline typename tools::promote_args<T>::type expm1(T x, const Policy& /* pol */)
0263 {
0264    typedef typename tools::promote_args<T>::type result_type;
0265    typedef typename policies::evaluation<result_type, Policy>::type value_type;
0266    typedef typename policies::precision<result_type, Policy>::type precision_type;
0267    typedef typename policies::normalise<
0268       Policy,
0269       policies::promote_float<false>,
0270       policies::promote_double<false>,
0271       policies::discrete_quantile<>,
0272       policies::assert_undefined<> >::type forwarding_policy;
0273
0274    typedef std::integral_constant<int,
0275       precision_type::value <= 0 ? 0 :
0276       precision_type::value <= 53 ? 53 :
0277       precision_type::value <= 64 ? 64 :
0278       precision_type::value <= 113 ? 113 : 0
0279    > tag_type;
0280
0281    detail::expm1_initializer<value_type, forwarding_policy, tag_type>::force_instantiate();
0282
0283    return policies::checked_narrowing_cast<result_type, forwarding_policy>(detail::expm1_imp(
0284       static_cast<value_type>(x),
0285       tag_type(), forwarding_policy()), "boost::math::expm1<%1%>(%1%)");
0286 }
0287
0288 #ifdef expm1
0289 #  ifndef BOOST_HAS_expm1
0290 #     define BOOST_HAS_expm1
0291 #  endif
0292 #  undef expm1
0293 #endif
0294
0295 #if defined(BOOST_HAS_EXPM1) && !(defined(__osf__) && defined(__DECCXX_VER))
0296 #  ifdef BOOST_MATH_USE_C99
0297 inline float expm1(float x, const policies::policy<>&){ return ::expm1f(x); }
0298 #     ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
0299 inline long double expm1(long double x, const policies::policy<>&){ return ::expm1l(x); }
0300 #     endif
0301 #  else
0302 inline float expm1(float x, const policies::policy<>&){ return static_cast<float>(::expm1(x)); }
0303 #  endif
0304 inline double expm1(double x, const policies::policy<>&){ return ::expm1(x); }
0305 #endif
0306
0307 template <class T>
0308 inline typename tools::promote_args<T>::type expm1(T x)
0309 {
0310    return expm1(x, policies::policy<>());
0311 }
0312
0313 } // namespace math
0314 } // namespace boost
0315
0316 #endif // BOOST_MATH_HYPOT_INCLUDED
0317
0318
0319
0320