math/special_functions/expm1.hpp

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