417 lines · plain
1// (C) Copyright John Maddock 2006.2// (C) Copyright Matt Borland 2024.3// Use, modification and distribution are subject to the4// Boost Software License, Version 1.0. (See accompanying file5// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)6 7#ifndef BOOST_MATH_EXPM1_INCLUDED8#define BOOST_MATH_EXPM1_INCLUDED9 10#ifdef _MSC_VER11#pragma once12#endif13 14#include <boost/math/tools/config.hpp>15 16#ifndef BOOST_MATH_HAS_NVRTC17 18#if defined __has_include19# if ((__cplusplus > 202002L) || (defined(_MSVC_LANG) && (_MSVC_LANG > 202002L)))20# if __has_include (<stdfloat>)21# include <stdfloat>22# endif23# endif24#endif25 26#include <boost/math/tools/series.hpp>27#include <boost/math/tools/precision.hpp>28#include <boost/math/tools/big_constant.hpp>29#include <boost/math/policies/error_handling.hpp>30#include <boost/math/tools/rational.hpp>31#include <boost/math/special_functions/math_fwd.hpp>32#include <boost/math/special_functions/fpclassify.hpp>33#include <boost/math/tools/assert.hpp>34#include <boost/math/tools/numeric_limits.hpp>35#include <boost/math/tools/type_traits.hpp>36#include <boost/math/tools/cstdint.hpp>37 38#if defined(__GNUC__) && defined(BOOST_MATH_USE_FLOAT128)39//40// This is the only way we can avoid41// warning: non-standard suffix on floating constant [-Wpedantic]42// when building with -Wall -pedantic. Neither __extension__43// nor #pragma diagnostic ignored work :(44//45#pragma GCC system_header46#endif47 48namespace boost {49 namespace math {50 51 namespace detail52 {53 // Functor expm1_series returns the next term in the Taylor series54 // x^k / k!55 // each time that operator() is invoked.56 //57 // LCOV_EXCL_START multiprecision case only, excluded from coverage analysis58 template <class T>59 struct expm1_series60 {61 typedef T result_type;62 63 BOOST_MATH_GPU_ENABLED expm1_series(T x)64 : k(0), m_x(x), m_term(1) {65 }66 67 BOOST_MATH_GPU_ENABLED T operator()()68 {69 ++k;70 m_term *= m_x;71 m_term /= k;72 return m_term;73 }74 75 BOOST_MATH_GPU_ENABLED int count()const76 {77 return k;78 }79 80 private:81 int k;82 const T m_x;83 T m_term;84 expm1_series(const expm1_series&) = delete;85 expm1_series& operator=(const expm1_series&) = delete;86 };87 88 //89 // Algorithm expm1 is part of C99, but is not yet provided by many compilers.90 //91 // This version uses a Taylor series expansion for 0.5 > |x| > epsilon.92 //93 template <class T, class Policy>94 T expm1_imp(T x, const boost::math::integral_constant<int, 0>&, const Policy& pol)95 {96 BOOST_MATH_STD_USING97 98 T a = fabs(x);99 if ((boost::math::isnan)(a))100 {101 return policies::raise_domain_error<T>("boost::math::expm1<%1%>(%1%)", "expm1 requires a finite argument, but got %1%", a, pol);102 }103 if (a > T(0.5f))104 {105 if (a >= tools::log_max_value<T>())106 {107 if (x > 0)108 return policies::raise_overflow_error<T>("boost::math::expm1<%1%>(%1%)", nullptr, pol);109 return -1;110 }111 return exp(x) - T(1);112 }113 if (a < tools::epsilon<T>())114 return x;115 detail::expm1_series<T> s(x);116 boost::math::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();117 118 T result = tools::sum_series(s, policies::get_epsilon<T, Policy>(), max_iter);119 120 policies::check_series_iterations<T>("boost::math::expm1<%1%>(%1%)", max_iter, pol);121 return result;122 }123 // LCOV_EXCL_STOP124 125 template <class T, class P>126 BOOST_MATH_GPU_ENABLED T expm1_imp(T x, const boost::math::integral_constant<int, 53>&, const P& pol)127 {128 BOOST_MATH_STD_USING129 130 T a = fabs(x);131 if ((boost::math::isnan)(a))132 {133 return policies::raise_domain_error<T>("boost::math::expm1<%1%>(%1%)", "expm1 requires a finite argument, but got %1%", a, pol);134 }135 if (a > T(0.5L))136 {137 if (a >= tools::log_max_value<T>())138 {139 if (x > 0)140 return policies::raise_overflow_error<T>("boost::math::expm1<%1%>(%1%)", nullptr, pol);141 return -1;142 }143 return exp(x) - T(1);144 }145 if (a < tools::epsilon<T>())146 return x;147 148 BOOST_MATH_STATIC const float Y = 0.10281276702880859e1f;149 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) };150 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) };151 152 T result = x * Y + x * tools::evaluate_polynomial(n, x) / tools::evaluate_polynomial(d, x);153 return result;154 }155 156 template <class T, class P>157 BOOST_MATH_GPU_ENABLED T expm1_imp(T x, const boost::math::integral_constant<int, 64>&, const P& pol)158 {159 BOOST_MATH_STD_USING160 161 T a = fabs(x);162 if ((boost::math::isnan)(a))163 {164 return policies::raise_domain_error<T>("boost::math::expm1<%1%>(%1%)", "expm1 requires a finite argument, but got %1%", a, pol);165 }166 if (a > T(0.5L))167 {168 if (a >= tools::log_max_value<T>())169 {170 if (x > 0)171 return policies::raise_overflow_error<T>("boost::math::expm1<%1%>(%1%)", nullptr, pol);172 return -1;173 }174 return exp(x) - T(1);175 }176 if (a < tools::epsilon<T>())177 return x;178 179 // LCOV_EXCL_START180 BOOST_MATH_STATIC const float Y = 0.10281276702880859375e1f;181 BOOST_MATH_STATIC const T n[] = {182 BOOST_MATH_BIG_CONSTANT(T, 64, -0.281276702880859375e-1),183 BOOST_MATH_BIG_CONSTANT(T, 64, 0.512980290285154286358e0),184 BOOST_MATH_BIG_CONSTANT(T, 64, -0.667758794592881019644e-1),185 BOOST_MATH_BIG_CONSTANT(T, 64, 0.131432469658444745835e-1),186 BOOST_MATH_BIG_CONSTANT(T, 64, -0.72303795326880286965e-3),187 BOOST_MATH_BIG_CONSTANT(T, 64, 0.447441185192951335042e-4),188 BOOST_MATH_BIG_CONSTANT(T, 64, -0.714539134024984593011e-6)189 };190 BOOST_MATH_STATIC const T d[] = {191 BOOST_MATH_BIG_CONSTANT(T, 64, 1.0),192 BOOST_MATH_BIG_CONSTANT(T, 64, -0.461477618025562520389e0),193 BOOST_MATH_BIG_CONSTANT(T, 64, 0.961237488025708540713e-1),194 BOOST_MATH_BIG_CONSTANT(T, 64, -0.116483957658204450739e-1),195 BOOST_MATH_BIG_CONSTANT(T, 64, 0.873308008461557544458e-3),196 BOOST_MATH_BIG_CONSTANT(T, 64, -0.387922804997682392562e-4),197 BOOST_MATH_BIG_CONSTANT(T, 64, 0.807473180049193557294e-6)198 };199 // LCOV_EXCL_STOP200 201 T result = x * Y + x * tools::evaluate_polynomial(n, x) / tools::evaluate_polynomial(d, x);202 return result;203 }204 205 template <class T, class P>206 BOOST_MATH_GPU_ENABLED T expm1_imp(T x, const boost::math::integral_constant<int, 113>&, const P& pol)207 {208 BOOST_MATH_STD_USING209 210 T a = fabs(x);211 if ((boost::math::isnan)(a))212 {213 return policies::raise_domain_error<T>("boost::math::expm1<%1%>(%1%)", "expm1 requires a finite argument, but got %1%", a, pol);214 }215 if (a > T(0.5L))216 {217 if (a >= tools::log_max_value<T>())218 {219 if (x > 0)220 return policies::raise_overflow_error<T>("boost::math::expm1<%1%>(%1%)", nullptr, pol);221 return -1;222 }223 return exp(x) - T(1);224 }225 if (a < tools::epsilon<T>())226 return x;227 228 // LCOV_EXCL_START229 static const float Y = 0.10281276702880859375e1f;230 static const T n[] = {231 BOOST_MATH_BIG_CONSTANT(T, 113, -0.28127670288085937499999999999999999854e-1),232 BOOST_MATH_BIG_CONSTANT(T, 113, 0.51278156911210477556524452177540792214e0),233 BOOST_MATH_BIG_CONSTANT(T, 113, -0.63263178520747096729500254678819588223e-1),234 BOOST_MATH_BIG_CONSTANT(T, 113, 0.14703285606874250425508446801230572252e-1),235 BOOST_MATH_BIG_CONSTANT(T, 113, -0.8675686051689527802425310407898459386e-3),236 BOOST_MATH_BIG_CONSTANT(T, 113, 0.88126359618291165384647080266133492399e-4),237 BOOST_MATH_BIG_CONSTANT(T, 113, -0.25963087867706310844432390015463138953e-5),238 BOOST_MATH_BIG_CONSTANT(T, 113, 0.14226691087800461778631773363204081194e-6),239 BOOST_MATH_BIG_CONSTANT(T, 113, -0.15995603306536496772374181066765665596e-8),240 BOOST_MATH_BIG_CONSTANT(T, 113, 0.45261820069007790520447958280473183582e-10)241 };242 static const T d[] = {243 BOOST_MATH_BIG_CONSTANT(T, 113, 1.0),244 BOOST_MATH_BIG_CONSTANT(T, 113, -0.45441264709074310514348137469214538853e0),245 BOOST_MATH_BIG_CONSTANT(T, 113, 0.96827131936192217313133611655555298106e-1),246 BOOST_MATH_BIG_CONSTANT(T, 113, -0.12745248725908178612540554584374876219e-1),247 BOOST_MATH_BIG_CONSTANT(T, 113, 0.11473613871583259821612766907781095472e-2),248 BOOST_MATH_BIG_CONSTANT(T, 113, -0.73704168477258911962046591907690764416e-4),249 BOOST_MATH_BIG_CONSTANT(T, 113, 0.34087499397791555759285503797256103259e-5),250 BOOST_MATH_BIG_CONSTANT(T, 113, -0.11114024704296196166272091230695179724e-6),251 BOOST_MATH_BIG_CONSTANT(T, 113, 0.23987051614110848595909588343223896577e-8),252 BOOST_MATH_BIG_CONSTANT(T, 113, -0.29477341859111589208776402638429026517e-10),253 BOOST_MATH_BIG_CONSTANT(T, 113, 0.13222065991022301420255904060628100924e-12)254 };255 // LCOV_EXCL_STOP256 257 T result = x * Y + x * tools::evaluate_polynomial(n, x) / tools::evaluate_polynomial(d, x);258 return result;259 }260 261 } // namespace detail262 263 template <class T, class Policy>264 BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<T>::type expm1(T x, const Policy& /* pol */)265 {266 typedef typename tools::promote_args<T>::type result_type;267 typedef typename policies::evaluation<result_type, Policy>::type value_type;268 typedef typename policies::precision<result_type, Policy>::type precision_type;269 typedef typename policies::normalise<270 Policy,271 policies::promote_float<false>,272 policies::promote_double<false>,273 policies::discrete_quantile<>,274 policies::assert_undefined<> >::type forwarding_policy;275 276 typedef boost::math::integral_constant<int,277 precision_type::value <= 0 ? 0 :278 precision_type::value <= 53 ? 53 :279 precision_type::value <= 64 ? 64 :280 precision_type::value <= 113 ? 113 : 0281 > tag_type;282 283 return policies::checked_narrowing_cast<result_type, forwarding_policy>(detail::expm1_imp(284 static_cast<value_type>(x),285 tag_type(), forwarding_policy()), "boost::math::expm1<%1%>(%1%)");286 }287 288 //289 // Since we now live in a post C++11 world, we can always defer to std::expm1 when appropriate:290 //291 template <class Policy>292 BOOST_MATH_GPU_ENABLED inline float expm1(float x, const Policy&)293 {294 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)295 {296 if ((boost::math::isnan)(x))297 return policies::raise_domain_error<float>("boost::math::expm1<%1%>(%1%)", "expm1 requires a finite argument, but got %1%", x, Policy());298 }299 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)300 {301 if (x >= tools::log_max_value<float>())302 return policies::raise_overflow_error<float>("boost::math::expm1<%1%>(%1%)", nullptr, Policy());303 }304 return std::expm1(x);305 }306#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS307 template <class Policy>308 inline long double expm1(long double x, const Policy&)309 {310 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)311 {312 if ((boost::math::isnan)(x))313 return policies::raise_domain_error<long double>("boost::math::expm1<%1%>(%1%)", "expm1 requires a finite argument, but got %1%", x, Policy());314 }315 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)316 {317 if (x >= tools::log_max_value<long double>())318 return policies::raise_overflow_error<long double>("boost::math::expm1<%1%>(%1%)", nullptr, Policy());319 }320 return std::expm1(x);321 }322#endif323 template <class Policy>324 BOOST_MATH_GPU_ENABLED inline double expm1(double x, const Policy&)325 {326 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)327 {328 if ((boost::math::isnan)(x))329 return policies::raise_domain_error<double>("boost::math::expm1<%1%>(%1%)", "expm1 requires a finite argument, but got %1%", x, Policy());330 }331 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)332 {333 if (x >= tools::log_max_value<double>())334 return policies::raise_overflow_error<double>("boost::math::expm1<%1%>(%1%)", nullptr, Policy());335 }336 return std::expm1(x);337 }338 339 template <class T>340 BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<T>::type expm1(T x)341 {342 return expm1(x, policies::policy<>());343 }344 //345 // Specific width floating point types:346 //347#ifdef __STDCPP_FLOAT32_T__348 template <class Policy>349 BOOST_MATH_GPU_ENABLED inline std::float32_t expm1(std::float32_t x, const Policy& pol)350 {351 return boost::math::expm1(static_cast<float>(x), pol);352 }353#endif354#ifdef __STDCPP_FLOAT64_T__355 template <class Policy>356 BOOST_MATH_GPU_ENABLED inline std::float64_t expm1(std::float64_t x, const Policy& pol)357 {358 return boost::math::expm1(static_cast<double>(x), pol);359 }360#endif361#ifdef __STDCPP_FLOAT128_T__362 template <class Policy>363 BOOST_MATH_GPU_ENABLED inline std::float128_t expm1(std::float128_t x, const Policy& pol)364 {365 if constexpr (std::numeric_limits<long double>::digits == std::numeric_limits<std::float128_t>::digits)366 {367 return boost::math::expm1(static_cast<long double>(x), pol);368 }369 else370 {371 return boost::math::detail::expm1_imp(x, boost::math::integral_constant<int, 113>(), pol);372 }373 }374#endif375} // namespace math376} // namespace boost377 378#else // Special handling for NVRTC 379 380namespace boost {381namespace math {382 383template <typename T>384BOOST_MATH_GPU_ENABLED auto expm1(T x)385{386 return ::expm1(x);387}388 389template <>390BOOST_MATH_GPU_ENABLED auto expm1(float x)391{392 return ::expm1f(x);393}394 395template <typename T, typename Policy>396BOOST_MATH_GPU_ENABLED auto expm1(T x, const Policy&)397{398 return ::expm1(x);399}400 401template <typename Policy>402BOOST_MATH_GPU_ENABLED auto expm1(float x, const Policy&)403{404 return ::expm1f(x);405}406 407} // Namespace math408} // Namespace boost409 410#endif // BOOST_MATH_HAS_NVRTC411 412#endif // BOOST_MATH_HYPOT_INCLUDED413 414 415 416 417