Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2026-08-14 08:52:14 

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 #if defined __has_include
0019 #  if ((__cplusplus > 202002L) || (defined(_MSVC_LANG) && (_MSVC_LANG > 202002L)))
0020 #    if __has_include (<stdfloat>)
0021 #    include <stdfloat>
0022 #    endif
0023 #  endif
0024 #endif
0025 
0026 #include <boost/math/tools/series.hpp>
0027 #include <boost/math/tools/precision.hpp>
0028 #include <boost/math/tools/big_constant.hpp>
0029 #include <boost/math/policies/error_handling.hpp>
0030 #include <boost/math/tools/rational.hpp>
0031 #include <boost/math/special_functions/math_fwd.hpp>
0032 #include <boost/math/special_functions/fpclassify.hpp>
0033 #include <boost/math/tools/assert.hpp>
0034 #include <boost/math/tools/numeric_limits.hpp>
0035 #include <boost/math/tools/type_traits.hpp>
0036 #include <boost/math/tools/cstdint.hpp>
0037 
0038 #if defined(__GNUC__) && defined(BOOST_MATH_USE_FLOAT128)
0039 //
0040 // This is the only way we can avoid
0041 // warning: non-standard suffix on floating constant [-Wpedantic]
0042 // when building with -Wall -pedantic.  Neither __extension__
0043 // nor #pragma diagnostic ignored work :(
0044 //
0045 #pragma GCC system_header
0046 #endif
0047 
0048 namespace boost {
0049    namespace math {
0050 
0051       namespace detail
0052       {
0053          // Functor expm1_series returns the next term in the Taylor series
0054          // x^k / k!
0055          // each time that operator() is invoked.
0056          //
0057          // LCOV_EXCL_START multiprecision case only, excluded from coverage analysis
0058          template <class T>
0059          struct expm1_series
0060          {
0061             typedef T result_type;
0062 
0063             BOOST_MATH_GPU_ENABLED expm1_series(T x)
0064                : k(0), m_x(x), m_term(1) {
0065             }
0066 
0067             BOOST_MATH_GPU_ENABLED T operator()()
0068             {
0069                ++k;
0070                m_term *= m_x;
0071                m_term /= k;
0072                return m_term;
0073             }
0074 
0075             BOOST_MATH_GPU_ENABLED int count()const
0076             {
0077                return k;
0078             }
0079 
0080          private:
0081             int k;
0082             const T m_x;
0083             T m_term;
0084             expm1_series(const expm1_series&) = delete;
0085             expm1_series& operator=(const expm1_series&) = delete;
0086          };
0087 
0088          //
0089          // Algorithm expm1 is part of C99, but is not yet provided by many compilers.
0090          //
0091          // This version uses a Taylor series expansion for 0.5 > |x| > epsilon.
0092          //
0093          template <class T, class Policy>
0094          T expm1_imp(T x, const boost::math::integral_constant<int, 0>&, const Policy& pol)
0095          {
0096             BOOST_MATH_STD_USING
0097 
0098                T a = fabs(x);
0099             if ((boost::math::isnan)(a))
0100             {
0101                return policies::raise_domain_error<T>("boost::math::expm1<%1%>(%1%)", "expm1 requires a finite argument, but got %1%", a, pol);
0102             }
0103             if (a > T(0.5f))
0104             {
0105                if (a >= tools::log_max_value<T>())
0106                {
0107                   if (x > 0)
0108                      return policies::raise_overflow_error<T>("boost::math::expm1<%1%>(%1%)", nullptr, pol);
0109                   return -1;
0110                }
0111                return exp(x) - T(1);
0112             }
0113             if (a < tools::epsilon<T>())
0114                return x;
0115             detail::expm1_series<T> s(x);
0116             boost::math::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
0117 
0118             T result = tools::sum_series(s, policies::get_epsilon<T, Policy>(), max_iter);
0119 
0120             policies::check_series_iterations<T>("boost::math::expm1<%1%>(%1%)", max_iter, pol);
0121             return result;
0122          }
0123          // LCOV_EXCL_STOP
0124 
0125          template <class T, class P>
0126          BOOST_MATH_GPU_ENABLED T expm1_imp(T x, const boost::math::integral_constant<int, 53>&, const P& pol)
0127          {
0128             BOOST_MATH_STD_USING
0129 
0130                T a = fabs(x);
0131             if ((boost::math::isnan)(a))
0132             {
0133                return policies::raise_domain_error<T>("boost::math::expm1<%1%>(%1%)", "expm1 requires a finite argument, but got %1%", a, pol);
0134             }
0135             if (a > T(0.5L))
0136             {
0137                if (a >= tools::log_max_value<T>())
0138                {
0139                   if (x > 0)
0140                      return policies::raise_overflow_error<T>("boost::math::expm1<%1%>(%1%)", nullptr, pol);
0141                   return -1;
0142                }
0143                return exp(x) - T(1);
0144             }
0145             if (a < tools::epsilon<T>())
0146                return x;
0147 
0148             BOOST_MATH_STATIC const float Y = 0.10281276702880859e1f;
0149             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) };
0150             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) };
0151 
0152             T result = x * Y + x * tools::evaluate_polynomial(n, x) / tools::evaluate_polynomial(d, x);
0153             return result;
0154          }
0155 
0156          template <class T, class P>
0157          BOOST_MATH_GPU_ENABLED T expm1_imp(T x, const boost::math::integral_constant<int, 64>&, const P& pol)
0158          {
0159             BOOST_MATH_STD_USING
0160 
0161                T a = fabs(x);
0162             if ((boost::math::isnan)(a))
0163             {
0164                return policies::raise_domain_error<T>("boost::math::expm1<%1%>(%1%)", "expm1 requires a finite argument, but got %1%", a, pol);
0165             }
0166             if (a > T(0.5L))
0167             {
0168                if (a >= tools::log_max_value<T>())
0169                {
0170                   if (x > 0)
0171                      return policies::raise_overflow_error<T>("boost::math::expm1<%1%>(%1%)", nullptr, pol);
0172                   return -1;
0173                }
0174                return exp(x) - T(1);
0175             }
0176             if (a < tools::epsilon<T>())
0177                return x;
0178 
0179             // LCOV_EXCL_START
0180             BOOST_MATH_STATIC const float Y = 0.10281276702880859375e1f;
0181             BOOST_MATH_STATIC const T n[] = {
0182                BOOST_MATH_BIG_CONSTANT(T, 64, -0.281276702880859375e-1),
0183                 BOOST_MATH_BIG_CONSTANT(T, 64, 0.512980290285154286358e0),
0184                 BOOST_MATH_BIG_CONSTANT(T, 64, -0.667758794592881019644e-1),
0185                 BOOST_MATH_BIG_CONSTANT(T, 64, 0.131432469658444745835e-1),
0186                 BOOST_MATH_BIG_CONSTANT(T, 64, -0.72303795326880286965e-3),
0187                 BOOST_MATH_BIG_CONSTANT(T, 64, 0.447441185192951335042e-4),
0188                 BOOST_MATH_BIG_CONSTANT(T, 64, -0.714539134024984593011e-6)
0189             };
0190             BOOST_MATH_STATIC const T d[] = {
0191                BOOST_MATH_BIG_CONSTANT(T, 64, 1.0),
0192                BOOST_MATH_BIG_CONSTANT(T, 64, -0.461477618025562520389e0),
0193                BOOST_MATH_BIG_CONSTANT(T, 64, 0.961237488025708540713e-1),
0194                BOOST_MATH_BIG_CONSTANT(T, 64, -0.116483957658204450739e-1),
0195                BOOST_MATH_BIG_CONSTANT(T, 64, 0.873308008461557544458e-3),
0196                BOOST_MATH_BIG_CONSTANT(T, 64, -0.387922804997682392562e-4),
0197                BOOST_MATH_BIG_CONSTANT(T, 64, 0.807473180049193557294e-6)
0198             };
0199             // LCOV_EXCL_STOP
0200 
0201             T result = x * Y + x * tools::evaluate_polynomial(n, x) / tools::evaluate_polynomial(d, x);
0202             return result;
0203          }
0204 
0205          template <class T, class P>
0206          BOOST_MATH_GPU_ENABLED T expm1_imp(T x, const boost::math::integral_constant<int, 113>&, const P& pol)
0207          {
0208             BOOST_MATH_STD_USING
0209 
0210                T a = fabs(x);
0211             if ((boost::math::isnan)(a))
0212             {
0213                return policies::raise_domain_error<T>("boost::math::expm1<%1%>(%1%)", "expm1 requires a finite argument, but got %1%", a, pol);
0214             }
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             // LCOV_EXCL_START
0229             static const float Y = 0.10281276702880859375e1f;
0230             static const T n[] = {
0231                BOOST_MATH_BIG_CONSTANT(T, 113, -0.28127670288085937499999999999999999854e-1),
0232                BOOST_MATH_BIG_CONSTANT(T, 113, 0.51278156911210477556524452177540792214e0),
0233                BOOST_MATH_BIG_CONSTANT(T, 113, -0.63263178520747096729500254678819588223e-1),
0234                BOOST_MATH_BIG_CONSTANT(T, 113, 0.14703285606874250425508446801230572252e-1),
0235                BOOST_MATH_BIG_CONSTANT(T, 113, -0.8675686051689527802425310407898459386e-3),
0236                BOOST_MATH_BIG_CONSTANT(T, 113, 0.88126359618291165384647080266133492399e-4),
0237                BOOST_MATH_BIG_CONSTANT(T, 113, -0.25963087867706310844432390015463138953e-5),
0238                BOOST_MATH_BIG_CONSTANT(T, 113, 0.14226691087800461778631773363204081194e-6),
0239                BOOST_MATH_BIG_CONSTANT(T, 113, -0.15995603306536496772374181066765665596e-8),
0240                BOOST_MATH_BIG_CONSTANT(T, 113, 0.45261820069007790520447958280473183582e-10)
0241             };
0242             static const T d[] = {
0243                BOOST_MATH_BIG_CONSTANT(T, 113, 1.0),
0244                BOOST_MATH_BIG_CONSTANT(T, 113, -0.45441264709074310514348137469214538853e0),
0245                BOOST_MATH_BIG_CONSTANT(T, 113, 0.96827131936192217313133611655555298106e-1),
0246                BOOST_MATH_BIG_CONSTANT(T, 113, -0.12745248725908178612540554584374876219e-1),
0247                BOOST_MATH_BIG_CONSTANT(T, 113, 0.11473613871583259821612766907781095472e-2),
0248                BOOST_MATH_BIG_CONSTANT(T, 113, -0.73704168477258911962046591907690764416e-4),
0249                BOOST_MATH_BIG_CONSTANT(T, 113, 0.34087499397791555759285503797256103259e-5),
0250                BOOST_MATH_BIG_CONSTANT(T, 113, -0.11114024704296196166272091230695179724e-6),
0251                BOOST_MATH_BIG_CONSTANT(T, 113, 0.23987051614110848595909588343223896577e-8),
0252                BOOST_MATH_BIG_CONSTANT(T, 113, -0.29477341859111589208776402638429026517e-10),
0253                BOOST_MATH_BIG_CONSTANT(T, 113, 0.13222065991022301420255904060628100924e-12)
0254             };
0255             // LCOV_EXCL_STOP
0256 
0257             T result = x * Y + x * tools::evaluate_polynomial(n, x) / tools::evaluate_polynomial(d, x);
0258             return result;
0259          }
0260 
0261       } // namespace detail
0262 
0263       template <class T, class Policy>
0264       BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<T>::type expm1(T x, const Policy& /* pol */)
0265       {
0266          typedef typename tools::promote_args<T>::type result_type;
0267          typedef typename policies::evaluation<result_type, Policy>::type value_type;
0268          typedef typename policies::precision<result_type, Policy>::type precision_type;
0269          typedef typename policies::normalise<
0270             Policy,
0271             policies::promote_float<false>,
0272             policies::promote_double<false>,
0273             policies::discrete_quantile<>,
0274             policies::assert_undefined<> >::type forwarding_policy;
0275 
0276          typedef boost::math::integral_constant<int,
0277             precision_type::value <= 0 ? 0 :
0278             precision_type::value <= 53 ? 53 :
0279             precision_type::value <= 64 ? 64 :
0280             precision_type::value <= 113 ? 113 : 0
0281          > tag_type;
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       //
0289       // Since we now live in a post C++11 world, we can always defer to std::expm1 when appropriate:
0290       //
0291       template <class Policy>
0292       BOOST_MATH_GPU_ENABLED inline float expm1(float x, const Policy&)
0293       {
0294          BOOST_MATH_IF_CONSTEXPR(Policy::domain_error_type::value != boost::math::policies::ignore_error && Policy::domain_error_type::value != boost::math::policies::errno_on_error)
0295          {
0296             if ((boost::math::isnan)(x))
0297                return policies::raise_domain_error<float>("boost::math::expm1<%1%>(%1%)", "expm1 requires a finite argument, but got %1%", x, Policy());
0298          }
0299          BOOST_MATH_IF_CONSTEXPR(Policy::overflow_error_type::value != boost::math::policies::ignore_error && Policy::overflow_error_type::value != boost::math::policies::errno_on_error)
0300          {
0301             if (x >= tools::log_max_value<float>())
0302                return policies::raise_overflow_error<float>("boost::math::expm1<%1%>(%1%)", nullptr, Policy());
0303          }
0304          return std::expm1(x);
0305       }
0306 #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
0307       template <class Policy>
0308       inline long double expm1(long double x, const Policy&)
0309       {
0310          BOOST_MATH_IF_CONSTEXPR(Policy::domain_error_type::value != boost::math::policies::ignore_error && Policy::domain_error_type::value != boost::math::policies::errno_on_error)
0311          {
0312             if ((boost::math::isnan)(x))
0313                return policies::raise_domain_error<long double>("boost::math::expm1<%1%>(%1%)", "expm1 requires a finite argument, but got %1%", x, Policy());
0314          }
0315          BOOST_MATH_IF_CONSTEXPR(Policy::overflow_error_type::value != boost::math::policies::ignore_error && Policy::overflow_error_type::value != boost::math::policies::errno_on_error)
0316          {
0317             if (x >= tools::log_max_value<long double>())
0318                return policies::raise_overflow_error<long double>("boost::math::expm1<%1%>(%1%)", nullptr, Policy());
0319          }
0320          return std::expm1(x);
0321       }
0322 #endif
0323       template <class Policy>
0324       BOOST_MATH_GPU_ENABLED inline double expm1(double x, const Policy&)
0325       {
0326          BOOST_MATH_IF_CONSTEXPR(Policy::domain_error_type::value != boost::math::policies::ignore_error && Policy::domain_error_type::value != boost::math::policies::errno_on_error)
0327          {
0328             if ((boost::math::isnan)(x))
0329                return policies::raise_domain_error<double>("boost::math::expm1<%1%>(%1%)", "expm1 requires a finite argument, but got %1%", x, Policy());
0330          }
0331          BOOST_MATH_IF_CONSTEXPR(Policy::overflow_error_type::value != boost::math::policies::ignore_error && Policy::overflow_error_type::value != boost::math::policies::errno_on_error)
0332          {
0333             if (x >= tools::log_max_value<double>())
0334                return policies::raise_overflow_error<double>("boost::math::expm1<%1%>(%1%)", nullptr, Policy());
0335          }
0336          return std::expm1(x);
0337       }
0338 
0339       template <class T>
0340       BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<T>::type expm1(T x)
0341       {
0342          return expm1(x, policies::policy<>());
0343       }
0344       //
0345       // Specific width floating point types:
0346       //
0347 #ifdef __STDCPP_FLOAT32_T__
0348       template <class Policy>
0349       BOOST_MATH_GPU_ENABLED inline std::float32_t expm1(std::float32_t x, const Policy& pol)
0350       {
0351          return boost::math::expm1(static_cast<float>(x), pol);
0352       }
0353 #endif
0354 #ifdef __STDCPP_FLOAT64_T__
0355       template <class Policy>
0356       BOOST_MATH_GPU_ENABLED inline std::float64_t expm1(std::float64_t x, const Policy& pol)
0357       {
0358          return boost::math::expm1(static_cast<double>(x), pol);
0359       }
0360 #endif
0361 #ifdef __STDCPP_FLOAT128_T__
0362       template <class Policy>
0363       BOOST_MATH_GPU_ENABLED inline std::float128_t expm1(std::float128_t x, const Policy& pol)
0364       {
0365          if constexpr (std::numeric_limits<long double>::digits == std::numeric_limits<std::float128_t>::digits)
0366          {
0367             return boost::math::expm1(static_cast<long double>(x), pol);
0368          }
0369          else
0370          {
0371             return boost::math::detail::expm1_imp(x, boost::math::integral_constant<int, 113>(), pol);
0372          }
0373       }
0374 #endif
0375 } // namespace math
0376 } // namespace boost
0377 
0378 #else // Special handling for NVRTC 
0379 
0380 namespace boost {
0381 namespace math {
0382 
0383 template <typename T>
0384 BOOST_MATH_GPU_ENABLED auto expm1(T x)
0385 {
0386    return ::expm1(x);
0387 }
0388 
0389 template <>
0390 BOOST_MATH_GPU_ENABLED auto expm1(float x)
0391 {
0392    return ::expm1f(x);
0393 }
0394 
0395 template <typename T, typename Policy>
0396 BOOST_MATH_GPU_ENABLED auto expm1(T x, const Policy&)
0397 {
0398    return ::expm1(x);
0399 }
0400 
0401 template <typename Policy>
0402 BOOST_MATH_GPU_ENABLED auto expm1(float x, const Policy&)
0403 {
0404    return ::expm1f(x);
0405 }
0406 
0407 } // Namespace math
0408 } // Namespace boost
0409 
0410 #endif // BOOST_MATH_HAS_NVRTC
0411 
0412 #endif // BOOST_MATH_HYPOT_INCLUDED
0413 
0414 
0415 
0416
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!