boost/math
2
0
Fork 0
mirror of https://github.com/boostorg/math.git synced 2026-01-19 04:22:09 +00:00
Code Issues Projects Releases Wiki Activity

Modernise log1p/expm1 to mostly use their std:: equivalents.

Browse Source 
Add support for <cstdfloat> types as well.
Extend tests for better coverage.
This commit is contained in:
jzmaddock
2025-04-19 13:35:23 +01:00
parent b0200e41d1
commit 54b057bb51
6 changed files with 523 additions and 392 deletions
Download Patch File Download Diff File Expand all files Collapse all files

634
include/boost/math/special_functions/expm1.hpp
View File

@@ -15,6 +15,14 @@
#ifndef BOOST_MATH_HAS_NVRTC #ifndef BOOST_MATH_HAS_NVRTC
#if defined __has_include
# if __cplusplus > 202002L || _MSVC_LANG > 202002L
# if __has_include (<stdfloat>)
# include <stdfloat>
# endif
# endif
#endif
#include <boost/math/tools/series.hpp> #include <boost/math/tools/series.hpp>
#include <boost/math/tools/precision.hpp> #include <boost/math/tools/precision.hpp>
#include <boost/math/tools/big_constant.hpp> #include <boost/math/tools/big_constant.hpp>
@@ -37,339 +45,331 @@
#pragma GCC system_header #pragma GCC system_header
#endif #endif
namespace boost{ namespace math{ namespace boost {
namespace math {
namespace detail namespace detail
{
// Functor expm1_series returns the next term in the Taylor series
// x^k / k!
// each time that operator() is invoked.
//
template <class T>
struct expm1_series
{
typedef T result_type;
BOOST_MATH_GPU_ENABLED expm1_series(T x)
: k(0), m_x(x), m_term(1) {}
BOOST_MATH_GPU_ENABLED T operator()()
{
++k;
m_term *= m_x;
m_term /= k;
return m_term;
}
BOOST_MATH_GPU_ENABLED int count()const
{
return k;
}
private:
int k;
const T m_x;
T m_term;
expm1_series(const expm1_series&) = delete;
expm1_series& operator=(const expm1_series&) = delete;
};
//
// Algorithm expm1 is part of C99, but is not yet provided by many compilers.
//
// This version uses a Taylor series expansion for 0.5 > |x| > epsilon.
//
template <class T, class Policy>
T expm1_imp(T x, const boost::math::integral_constant<int, 0>&, const Policy& pol)
{
BOOST_MATH_STD_USING
T a = fabs(x);
if((boost::math::isnan)(a))
{
return policies::raise_domain_error<T>("boost::math::expm1<%1%>(%1%)", "expm1 requires a finite argument, but got %1%", a, pol);
}
if(a > T(0.5f))
{
if(a >= tools::log_max_value<T>())
{ {
if(x > 0) // Functor expm1_series returns the next term in the Taylor series
return policies::raise_overflow_error<T>("boost::math::expm1<%1%>(%1%)", nullptr, pol); // x^k / k!
return -1; // each time that operator() is invoked.
} //
return exp(x) - T(1); template <class T>
} struct expm1_series
if(a < tools::epsilon<T>()) {
return x; typedef T result_type;
detail::expm1_series<T> s(x);
boost::math::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
T result = tools::sum_series(s, policies::get_epsilon<T, Policy>(), max_iter); BOOST_MATH_GPU_ENABLED expm1_series(T x)
: k(0), m_x(x), m_term(1) {
}
policies::check_series_iterations<T>("boost::math::expm1<%1%>(%1%)", max_iter, pol); BOOST_MATH_GPU_ENABLED T operator()()
return result; {
} ++k;
m_term *= m_x;
m_term /= k;
return m_term;
}
template <class T, class P> BOOST_MATH_GPU_ENABLED int count()const
BOOST_MATH_GPU_ENABLED T expm1_imp(T x, const boost::math::integral_constant<int, 53>&, const P& pol) {
{ return k;
BOOST_MATH_STD_USING }
T a = fabs(x); private:
if ((boost::math::isnan)(a)) int k;
{ const T m_x;
return policies::raise_domain_error<T>("boost::math::expm1<%1%>(%1%)", "expm1 requires a finite argument, but got %1%", a, pol); T m_term;
} expm1_series(const expm1_series&) = delete;
if(a > T(0.5L)) expm1_series& operator=(const expm1_series&) = delete;
{ };
if(a >= tools::log_max_value<T>())
//
// Algorithm expm1 is part of C99, but is not yet provided by many compilers.
//
// This version uses a Taylor series expansion for 0.5 > |x| > epsilon.
//
template <class T, class Policy>
T expm1_imp(T x, const boost::math::integral_constant<int, 0>&, const Policy& pol)
{
BOOST_MATH_STD_USING
T a = fabs(x);
if ((boost::math::isnan)(a))
{
return policies::raise_domain_error<T>("boost::math::expm1<%1%>(%1%)", "expm1 requires a finite argument, but got %1%", a, pol);
}
if (a > T(0.5f))
{
if (a >= tools::log_max_value<T>())
{
if (x > 0)
return policies::raise_overflow_error<T>("boost::math::expm1<%1%>(%1%)", nullptr, pol);
return -1;
}
return exp(x) - T(1);
}
if (a < tools::epsilon<T>())
return x;
detail::expm1_series<T> s(x);
boost::math::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
T result = tools::sum_series(s, policies::get_epsilon<T, Policy>(), max_iter);
policies::check_series_iterations<T>("boost::math::expm1<%1%>(%1%)", max_iter, pol);
return result;
}
template <class T, class P>
BOOST_MATH_GPU_ENABLED T expm1_imp(T x, const boost::math::integral_constant<int, 53>&, const P& pol)
{
BOOST_MATH_STD_USING
T a = fabs(x);
if ((boost::math::isnan)(a))
{
return policies::raise_domain_error<T>("boost::math::expm1<%1%>(%1%)", "expm1 requires a finite argument, but got %1%", a, pol);
}
if (a > T(0.5L))
{
if (a >= tools::log_max_value<T>())
{
if (x > 0)
return policies::raise_overflow_error<T>("boost::math::expm1<%1%>(%1%)", nullptr, pol);
return -1;
}
return exp(x) - T(1);
}
if (a < tools::epsilon<T>())
return x;
BOOST_MATH_STATIC const float Y = 0.10281276702880859e1f;
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) };
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) };
T result = x * Y + x * tools::evaluate_polynomial(n, x) / tools::evaluate_polynomial(d, x);
return result;
}
template <class T, class P>
BOOST_MATH_GPU_ENABLED T expm1_imp(T x, const boost::math::integral_constant<int, 64>&, const P& pol)
{
BOOST_MATH_STD_USING
T a = fabs(x);
if ((boost::math::isnan)(a))
{
return policies::raise_domain_error<T>("boost::math::expm1<%1%>(%1%)", "expm1 requires a finite argument, but got %1%", a, pol);
}
if (a > T(0.5L))
{
if (a >= tools::log_max_value<T>())
{
if (x > 0)
return policies::raise_overflow_error<T>("boost::math::expm1<%1%>(%1%)", nullptr, pol);
return -1;
}
return exp(x) - T(1);
}
if (a < tools::epsilon<T>())
return x;
// LCOV_EXCL_START
BOOST_MATH_STATIC const float Y = 0.10281276702880859375e1f;
BOOST_MATH_STATIC const T n[] = {
BOOST_MATH_BIG_CONSTANT(T, 64, -0.281276702880859375e-1),
BOOST_MATH_BIG_CONSTANT(T, 64, 0.512980290285154286358e0),
BOOST_MATH_BIG_CONSTANT(T, 64, -0.667758794592881019644e-1),
BOOST_MATH_BIG_CONSTANT(T, 64, 0.131432469658444745835e-1),
BOOST_MATH_BIG_CONSTANT(T, 64, -0.72303795326880286965e-3),
BOOST_MATH_BIG_CONSTANT(T, 64, 0.447441185192951335042e-4),
BOOST_MATH_BIG_CONSTANT(T, 64, -0.714539134024984593011e-6)
};
BOOST_MATH_STATIC const T d[] = {
BOOST_MATH_BIG_CONSTANT(T, 64, 1.0),
BOOST_MATH_BIG_CONSTANT(T, 64, -0.461477618025562520389e0),
BOOST_MATH_BIG_CONSTANT(T, 64, 0.961237488025708540713e-1),
BOOST_MATH_BIG_CONSTANT(T, 64, -0.116483957658204450739e-1),
BOOST_MATH_BIG_CONSTANT(T, 64, 0.873308008461557544458e-3),
BOOST_MATH_BIG_CONSTANT(T, 64, -0.387922804997682392562e-4),
BOOST_MATH_BIG_CONSTANT(T, 64, 0.807473180049193557294e-6)
};
// LCOV_EXCL_STOP
T result = x * Y + x * tools::evaluate_polynomial(n, x) / tools::evaluate_polynomial(d, x);
return result;
}
template <class T, class P>
BOOST_MATH_GPU_ENABLED T expm1_imp(T x, const boost::math::integral_constant<int, 113>&, const P& pol)
{
BOOST_MATH_STD_USING
T a = fabs(x);
if ((boost::math::isnan)(a))
{
return policies::raise_domain_error<T>("boost::math::expm1<%1%>(%1%)", "expm1 requires a finite argument, but got %1%", a, pol);
}
if (a > T(0.5L))
{
if (a >= tools::log_max_value<T>())
{
if (x > 0)
return policies::raise_overflow_error<T>("boost::math::expm1<%1%>(%1%)", nullptr, pol);
return -1;
}
return exp(x) - T(1);
}
if (a < tools::epsilon<T>())
return x;
// LCOV_EXCL_START
static const float Y = 0.10281276702880859375e1f;
static const T n[] = {
BOOST_MATH_BIG_CONSTANT(T, 113, -0.28127670288085937499999999999999999854e-1),
BOOST_MATH_BIG_CONSTANT(T, 113, 0.51278156911210477556524452177540792214e0),
BOOST_MATH_BIG_CONSTANT(T, 113, -0.63263178520747096729500254678819588223e-1),
BOOST_MATH_BIG_CONSTANT(T, 113, 0.14703285606874250425508446801230572252e-1),
BOOST_MATH_BIG_CONSTANT(T, 113, -0.8675686051689527802425310407898459386e-3),
BOOST_MATH_BIG_CONSTANT(T, 113, 0.88126359618291165384647080266133492399e-4),
BOOST_MATH_BIG_CONSTANT(T, 113, -0.25963087867706310844432390015463138953e-5),
BOOST_MATH_BIG_CONSTANT(T, 113, 0.14226691087800461778631773363204081194e-6),
BOOST_MATH_BIG_CONSTANT(T, 113, -0.15995603306536496772374181066765665596e-8),
BOOST_MATH_BIG_CONSTANT(T, 113, 0.45261820069007790520447958280473183582e-10)
};
static const T d[] = {
BOOST_MATH_BIG_CONSTANT(T, 113, 1.0),
BOOST_MATH_BIG_CONSTANT(T, 113, -0.45441264709074310514348137469214538853e0),
BOOST_MATH_BIG_CONSTANT(T, 113, 0.96827131936192217313133611655555298106e-1),
BOOST_MATH_BIG_CONSTANT(T, 113, -0.12745248725908178612540554584374876219e-1),
BOOST_MATH_BIG_CONSTANT(T, 113, 0.11473613871583259821612766907781095472e-2),
BOOST_MATH_BIG_CONSTANT(T, 113, -0.73704168477258911962046591907690764416e-4),
BOOST_MATH_BIG_CONSTANT(T, 113, 0.34087499397791555759285503797256103259e-5),
BOOST_MATH_BIG_CONSTANT(T, 113, -0.11114024704296196166272091230695179724e-6),
BOOST_MATH_BIG_CONSTANT(T, 113, 0.23987051614110848595909588343223896577e-8),
BOOST_MATH_BIG_CONSTANT(T, 113, -0.29477341859111589208776402638429026517e-10),
BOOST_MATH_BIG_CONSTANT(T, 113, 0.13222065991022301420255904060628100924e-12)
};
// LCOV_EXCL_STOP
T result = x * Y + x * tools::evaluate_polynomial(n, x) / tools::evaluate_polynomial(d, x);
return result;
}
} // namespace detail
template <class T, class Policy>
BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<T>::type expm1(T x, const Policy& /* pol */)
{ {
if(x > 0) typedef typename tools::promote_args<T>::type result_type;
return policies::raise_overflow_error<T>("boost::math::expm1<%1%>(%1%)", nullptr, pol); typedef typename policies::evaluation<result_type, Policy>::type value_type;
return -1; typedef typename policies::precision<result_type, Policy>::type precision_type;
typedef typename policies::normalise<
Policy,
policies::promote_float<false>,
policies::promote_double<false>,
policies::discrete_quantile<>,
policies::assert_undefined<> >::type forwarding_policy;
typedef boost::math::integral_constant<int,
precision_type::value <= 0 ? 0 :
precision_type::value <= 53 ? 53 :
precision_type::value <= 64 ? 64 :
precision_type::value <= 113 ? 113 : 0
> tag_type;
return policies::checked_narrowing_cast<result_type, forwarding_policy>(detail::expm1_imp(
static_cast<value_type>(x),
tag_type(), forwarding_policy()), "boost::math::expm1<%1%>(%1%)");
} }
return exp(x) - T(1);
}
if(a < tools::epsilon<T>())
return x;
BOOST_MATH_STATIC const float Y = 0.10281276702880859e1f; //
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) }; // Since we now live in a post C++11 world, we can always defer to std::expm1 when appropriate:
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) }; //
template <class Policy>
T result = x * Y + x * tools::evaluate_polynomial(n, x) / tools::evaluate_polynomial(d, x); BOOST_MATH_GPU_ENABLED inline float expm1(float x, const Policy&)
return result;
}
template <class T, class P>
BOOST_MATH_GPU_ENABLED T expm1_imp(T x, const boost::math::integral_constant<int, 64>&, const P& pol)
{
BOOST_MATH_STD_USING
T a = fabs(x);
if ((boost::math::isnan)(a))
{
return policies::raise_domain_error<T>("boost::math::expm1<%1%>(%1%)", "expm1 requires a finite argument, but got %1%", a, pol);
}
if(a > T(0.5L))
{
if(a >= tools::log_max_value<T>())
{ {
if(x > 0) 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)
return policies::raise_overflow_error<T>("boost::math::expm1<%1%>(%1%)", nullptr, pol); {
return -1; if ((boost::math::isnan)(x))
return policies::raise_domain_error<float>("boost::math::expm1<%1%>(%1%)", "expm1 requires a finite argument, but got %1%", x, Policy());
}
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)
{
if (x >= tools::log_max_value<float>())
return policies::raise_overflow_error<float>("boost::math::expm1<%1%>(%1%)", nullptr, Policy());
}
return std::expm1(x);
} }
return exp(x) - T(1); #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
} template <class Policy>
if(a < tools::epsilon<T>()) inline long double expm1(long double x, const Policy&)
return x;
// LCOV_EXCL_START
BOOST_MATH_STATIC const float Y = 0.10281276702880859375e1f;
BOOST_MATH_STATIC const T n[] = {
BOOST_MATH_BIG_CONSTANT(T, 64, -0.281276702880859375e-1),
BOOST_MATH_BIG_CONSTANT(T, 64, 0.512980290285154286358e0),
BOOST_MATH_BIG_CONSTANT(T, 64, -0.667758794592881019644e-1),
BOOST_MATH_BIG_CONSTANT(T, 64, 0.131432469658444745835e-1),
BOOST_MATH_BIG_CONSTANT(T, 64, -0.72303795326880286965e-3),
BOOST_MATH_BIG_CONSTANT(T, 64, 0.447441185192951335042e-4),
BOOST_MATH_BIG_CONSTANT(T, 64, -0.714539134024984593011e-6)
};
BOOST_MATH_STATIC const T d[] = {
BOOST_MATH_BIG_CONSTANT(T, 64, 1.0),
BOOST_MATH_BIG_CONSTANT(T, 64, -0.461477618025562520389e0),
BOOST_MATH_BIG_CONSTANT(T, 64, 0.961237488025708540713e-1),
BOOST_MATH_BIG_CONSTANT(T, 64, -0.116483957658204450739e-1),
BOOST_MATH_BIG_CONSTANT(T, 64, 0.873308008461557544458e-3),
BOOST_MATH_BIG_CONSTANT(T, 64, -0.387922804997682392562e-4),
BOOST_MATH_BIG_CONSTANT(T, 64, 0.807473180049193557294e-6)
};
// LCOV_EXCL_STOP
T result = x * Y + x * tools::evaluate_polynomial(n, x) / tools::evaluate_polynomial(d, x);
return result;
}
template <class T, class P>
BOOST_MATH_GPU_ENABLED T expm1_imp(T x, const boost::math::integral_constant<int, 113>&, const P& pol)
{
BOOST_MATH_STD_USING
T a = fabs(x);
if ((boost::math::isnan)(a))
{
return policies::raise_domain_error<T>("boost::math::expm1<%1%>(%1%)", "expm1 requires a finite argument, but got %1%", a, pol);
}
if(a > T(0.5L))
{
if(a >= tools::log_max_value<T>())
{ {
if(x > 0) 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)
return policies::raise_overflow_error<T>("boost::math::expm1<%1%>(%1%)", nullptr, pol); {
return -1; if ((boost::math::isnan)(x))
return policies::raise_domain_error<long double>("boost::math::expm1<%1%>(%1%)", "expm1 requires a finite argument, but got %1%", x, Policy());
}
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)
{
if (x >= tools::log_max_value<long double>())
return policies::raise_overflow_error<long double>("boost::math::expm1<%1%>(%1%)", nullptr, Policy());
}
return std::expm1(x);
} }
return exp(x) - T(1);
}
if(a < tools::epsilon<T>())
return x;
// LCOV_EXCL_START
static const float Y = 0.10281276702880859375e1f;
static const T n[] = {
BOOST_MATH_BIG_CONSTANT(T, 113, -0.28127670288085937499999999999999999854e-1),
BOOST_MATH_BIG_CONSTANT(T, 113, 0.51278156911210477556524452177540792214e0),
BOOST_MATH_BIG_CONSTANT(T, 113, -0.63263178520747096729500254678819588223e-1),
BOOST_MATH_BIG_CONSTANT(T, 113, 0.14703285606874250425508446801230572252e-1),
BOOST_MATH_BIG_CONSTANT(T, 113, -0.8675686051689527802425310407898459386e-3),
BOOST_MATH_BIG_CONSTANT(T, 113, 0.88126359618291165384647080266133492399e-4),
BOOST_MATH_BIG_CONSTANT(T, 113, -0.25963087867706310844432390015463138953e-5),
BOOST_MATH_BIG_CONSTANT(T, 113, 0.14226691087800461778631773363204081194e-6),
BOOST_MATH_BIG_CONSTANT(T, 113, -0.15995603306536496772374181066765665596e-8),
BOOST_MATH_BIG_CONSTANT(T, 113, 0.45261820069007790520447958280473183582e-10)
};
static const T d[] = {
BOOST_MATH_BIG_CONSTANT(T, 113, 1.0),
BOOST_MATH_BIG_CONSTANT(T, 113, -0.45441264709074310514348137469214538853e0),
BOOST_MATH_BIG_CONSTANT(T, 113, 0.96827131936192217313133611655555298106e-1),
BOOST_MATH_BIG_CONSTANT(T, 113, -0.12745248725908178612540554584374876219e-1),
BOOST_MATH_BIG_CONSTANT(T, 113, 0.11473613871583259821612766907781095472e-2),
BOOST_MATH_BIG_CONSTANT(T, 113, -0.73704168477258911962046591907690764416e-4),
BOOST_MATH_BIG_CONSTANT(T, 113, 0.34087499397791555759285503797256103259e-5),
BOOST_MATH_BIG_CONSTANT(T, 113, -0.11114024704296196166272091230695179724e-6),
BOOST_MATH_BIG_CONSTANT(T, 113, 0.23987051614110848595909588343223896577e-8),
BOOST_MATH_BIG_CONSTANT(T, 113, -0.29477341859111589208776402638429026517e-10),
BOOST_MATH_BIG_CONSTANT(T, 113, 0.13222065991022301420255904060628100924e-12)
};
// LCOV_EXCL_STOP
T result = x * Y + x * tools::evaluate_polynomial(n, x) / tools::evaluate_polynomial(d, x);
return result;
}
} // namespace detail
template <class T, class Policy>
BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<T>::type expm1(T x, const Policy& /* pol */)
{
typedef typename tools::promote_args<T>::type result_type;
typedef typename policies::evaluation<result_type, Policy>::type value_type;
typedef typename policies::precision<result_type, Policy>::type precision_type;
typedef typename policies::normalise<
Policy,
policies::promote_float<false>,
policies::promote_double<false>,
policies::discrete_quantile<>,
policies::assert_undefined<> >::type forwarding_policy;
typedef boost::math::integral_constant<int,
precision_type::value <= 0 ? 0 :
precision_type::value <= 53 ? 53 :
precision_type::value <= 64 ? 64 :
precision_type::value <= 113 ? 113 : 0
> tag_type;
return policies::checked_narrowing_cast<result_type, forwarding_policy>(detail::expm1_imp(
static_cast<value_type>(x),
tag_type(), forwarding_policy()), "boost::math::expm1<%1%>(%1%)");
}
#ifdef expm1
# ifndef BOOST_HAS_expm1
# define BOOST_HAS_expm1
# endif
# undef expm1
#endif #endif
template <class Policy>
BOOST_MATH_GPU_ENABLED inline double expm1(double x, const Policy&)
{
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)
{
if ((boost::math::isnan)(x))
return policies::raise_domain_error<double>("boost::math::expm1<%1%>(%1%)", "expm1 requires a finite argument, but got %1%", x, Policy());
}
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)
{
if (x >= tools::log_max_value<double>())
return policies::raise_overflow_error<double>("boost::math::expm1<%1%>(%1%)", nullptr, Policy());
}
return std::expm1(x);
}
#if defined(BOOST_HAS_EXPM1) && !(defined(__osf__) && defined(__DECCXX_VER)) template <class T>
# ifdef BOOST_MATH_USE_C99 BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<T>::type expm1(T x)
template <class Policy> {
BOOST_MATH_GPU_ENABLED inline float expm1(float x, const Policy&) return expm1(x, policies::policy<>());
{ }
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) //
{ // Specific width floating point types:
if((boost::math::isnan)(x)) //
return policies::raise_domain_error<float>("boost::math::expm1<%1%>(%1%)", "expm1 requires a finite argument, but got %1%", x, Policy()); #ifdef __STDCPP_FLOAT32_T__
} template <class Policy>
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) BOOST_MATH_GPU_ENABLED inline std::float32_t expm1(std::float32_t x, const Policy& pol)
{ {
if (x >= tools::log_max_value<float>()) return boost::math::expm1(static_cast<float>(x), pol);
return policies::raise_overflow_error<float>("boost::math::expm1<%1%>(%1%)", nullptr, Policy()); }
} #endif
return ::expm1f(x); #ifdef __STDCPP_FLOAT64_T__
} template <class Policy>
# ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS BOOST_MATH_GPU_ENABLED inline std::float64_t expm1(std::float64_t x, const Policy& pol)
template <class Policy> {
inline long double expm1(long double x, const Policy&) return boost::math::expm1(static_cast<double>(x), pol);
{ }
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) #endif
{ #ifdef __STDCPP_FLOAT128_T__
if ((boost::math::isnan)(x)) template <class Policy>
return policies::raise_domain_error<long double>("boost::math::expm1<%1%>(%1%)", "expm1 requires a finite argument, but got %1%", x, Policy()); BOOST_MATH_GPU_ENABLED inline std::float128_t expm1(std::float128_t x, const Policy& pol)
} {
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) if constexpr (std::numeric_limits<long double>::digits == std::numeric_limits<std::float128_t>::digits)
{ {
if (x >= tools::log_max_value<long double>()) return boost::math::expm1(static_cast<long double>(x), pol);
return policies::raise_overflow_error<long double>("boost::math::expm1<%1%>(%1%)", nullptr, Policy()); }
} else
return ::expm1l(x); {
} return boost::math::detail::expm1_imp(x, boost::math::integral_constant<int, 113>(), pol);
# endif }
# else }
template <class Policy>
inline float expm1(float x, const Policy&)
{
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)
{
if ((boost::math::isnan)(x))
return policies::raise_domain_error<float>("boost::math::expm1<%1%>(%1%)", "expm1 requires a finite argument, but got %1%", x, Policy());
}
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)
{
if (x >= tools::log_max_value<float>())
return policies::raise_overflow_error<float>("boost::math::expm1<%1%>(%1%)", nullptr, Policy());
}
return static_cast<float>(::expm1(x));
}
template <class Policy>
inline typename std::enable_if<sizeof(double) == sizeof(long double), long double>::type expm1(long double x, const Policy&)
{
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)
{
if ((boost::math::isnan)(x))
return policies::raise_domain_error<long double>("boost::math::expm1<%1%>(%1%)", "expm1 requires a finite argument, but got %1%", x, Policy());
}
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)
{
if (x >= tools::log_max_value<float>())
return policies::raise_overflow_error<long double>("boost::math::expm1<%1%>(%1%)", nullptr, Policy());
}
return ::expm1(x);
}
# endif
template <class Policy>
BOOST_MATH_GPU_ENABLED inline double expm1(double x, const Policy&)
{
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)
{
if ((boost::math::isnan)(x))
return policies::raise_domain_error<double>("boost::math::expm1<%1%>(%1%)", "expm1 requires a finite argument, but got %1%", x, Policy());
}
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)
{
if (x >= tools::log_max_value<double>())
return policies::raise_overflow_error<double>("boost::math::expm1<%1%>(%1%)", nullptr, Policy());
}
return ::expm1(x);
}
#endif #endif
template <class T>
BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<T>::type expm1(T x)
{
return expm1(x, policies::policy<>());
}
} // namespace math } // namespace math
} // namespace boost } // namespace boost

115
include/boost/math/special_functions/log1p.hpp
View File

@@ -12,6 +12,14 @@
#pragma warning(disable:4702) // Unreachable code (release mode only warning) #pragma warning(disable:4702) // Unreachable code (release mode only warning)
#endif #endif
#if defined __has_include
# if __cplusplus > 202002L || _MSVC_LANG > 202002L
# if __has_include (<stdfloat>)
# include <stdfloat>
# endif
# endif
#endif
#include <boost/math/tools/config.hpp> #include <boost/math/tools/config.hpp>
#include <boost/math/tools/series.hpp> #include <boost/math/tools/series.hpp>
#include <boost/math/tools/rational.hpp> #include <boost/math/tools/rational.hpp>
@@ -289,15 +297,6 @@ BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<T>::type log1p(T x, c
detail::log1p_imp(static_cast<value_type>(x), forwarding_policy(), tag_type()), "boost::math::log1p<%1%>(%1%)"); detail::log1p_imp(static_cast<value_type>(x), forwarding_policy(), tag_type()), "boost::math::log1p<%1%>(%1%)");
} }
#ifdef log1p
# ifndef BOOST_HAS_LOG1P
# define BOOST_HAS_LOG1P
# endif
# undef log1p
#endif
#if defined(BOOST_HAS_LOG1P) && !(defined(__osf__) && defined(__DECCXX_VER))
# ifdef BOOST_MATH_USE_C99
template <class Policy> template <class Policy>
BOOST_MATH_GPU_ENABLED inline float log1p(float x, const Policy& pol) BOOST_MATH_GPU_ENABLED inline float log1p(float x, const Policy& pol)
{ {
@@ -307,7 +306,7 @@ BOOST_MATH_GPU_ENABLED inline float log1p(float x, const Policy& pol)
if(x == -1) if(x == -1)
return -policies::raise_overflow_error<float>( return -policies::raise_overflow_error<float>(
"log1p<%1%>(%1%)", nullptr, pol); "log1p<%1%>(%1%)", nullptr, pol);
return ::log1pf(x); return std::log1p(x);
} }
#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
template <class Policy> template <class Policy>
@@ -319,20 +318,7 @@ BOOST_MATH_GPU_ENABLED inline long double log1p(long double x, const Policy& pol
if(x == -1) if(x == -1)
return -policies::raise_overflow_error<long double>( return -policies::raise_overflow_error<long double>(
"log1p<%1%>(%1%)", nullptr, pol); "log1p<%1%>(%1%)", nullptr, pol);
return ::log1pl(x); return std::log1p(x);
}
#endif
#else
template <class Policy>
inline float log1p(float x, const Policy& pol)
{
if(x < -1)
return policies::raise_domain_error<float>(
"log1p<%1%>(%1%)", "log1p(x) requires x > -1, but got x = %1%.", x, pol);
if(x == -1)
return -policies::raise_overflow_error<float>(
"log1p<%1%>(%1%)", nullptr, pol);
return ::log1p(x);
} }
#endif #endif
template <class Policy> template <class Policy>
@@ -344,56 +330,8 @@ BOOST_MATH_GPU_ENABLED inline double log1p(double x, const Policy& pol)
if(x == -1) if(x == -1)
return -policies::raise_overflow_error<double>( return -policies::raise_overflow_error<double>(
"log1p<%1%>(%1%)", nullptr, pol); "log1p<%1%>(%1%)", nullptr, pol);
return ::log1p(x); return std::log1p(x);
} }
#elif defined(_MSC_VER) && (BOOST_MSVC >= 1400)
//
// You should only enable this branch if you are absolutely sure
// that your compilers optimizer won't mess this code up!!
// Currently tested with VC8 and Intel 9.1.
//
template <class Policy>
inline double log1p(double x, const Policy& pol)
{
if(x < -1)
return policies::raise_domain_error<double>(
"log1p<%1%>(%1%)", "log1p(x) requires x > -1, but got x = %1%.", x, pol);
if(x == -1)
return -policies::raise_overflow_error<double>(
"log1p<%1%>(%1%)", nullptr, pol);
double u = 1+x;
if(u == 1.0)
return x;
else
return ::log(u)*(x/(u-1.0));
}
template <class Policy>
inline float log1p(float x, const Policy& pol)
{
return static_cast<float>(boost::math::log1p(static_cast<double>(x), pol));
}
#ifndef _WIN32_WCE
//
// For some reason this fails to compile under WinCE...
// Needs more investigation.
//
template <class Policy>
inline long double log1p(long double x, const Policy& pol)
{
if(x < -1)
return policies::raise_domain_error<long double>(
"log1p<%1%>(%1%)", "log1p(x) requires x > -1, but got x = %1%.", x, pol);
if(x == -1)
return -policies::raise_overflow_error<long double>(
"log1p<%1%>(%1%)", nullptr, pol);
long double u = 1+x;
if(u == 1.0)
return x;
else
return ::logl(u)*(x/(u-1.0));
}
#endif
#endif
template <class T> template <class T>
BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<T>::type log1p(T x) BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<T>::type log1p(T x)
@@ -441,6 +379,37 @@ BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<T>::type log1pmx(T x)
return log1pmx(x, policies::policy<>()); return log1pmx(x, policies::policy<>());
} }
//
// Specific width floating point types:
//
#ifdef __STDCPP_FLOAT32_T__
template <class Policy>
BOOST_MATH_GPU_ENABLED inline std::float32_t log1p(std::float32_t x, const Policy& pol)
{
return boost::math::log1p(static_cast<float>(x), pol);
}
#endif
#ifdef __STDCPP_FLOAT64_T__
template <class Policy>
BOOST_MATH_GPU_ENABLED inline std::float64_t log1p(std::float64_t x, const Policy& pol)
{
return boost::math::log1p(static_cast<double>(x), pol);
}
#endif
#ifdef __STDCPP_FLOAT128_T__
template <class Policy>
BOOST_MATH_GPU_ENABLED inline std::float128_t log1p(std::float128_t x, const Policy& pol)
{
if constexpr (std::numeric_limits<long double>::digits == std::numeric_limits<std::float128_t>::digits)
{
return boost::math::log1p(static_cast<long double>(x), pol);
}
else
{
return boost::math::detail::log1p_imp(x, pol, boost::math::integral_constant<int, 113>());
}
}
#endif
} // namespace math } // namespace math
} // namespace boost } // namespace boost

4
include_private/boost/math/tools/test.hpp
View File

@@ -167,7 +167,7 @@ test_result<typename calculate_result_type<A>::value_type> test(const A& a, F1 t
print_row(row, std::cerr); print_row(row, std::cerr);
BOOST_ERROR("Unexpected non-finite result"); BOOST_ERROR("Unexpected non-finite result");
} }
if(err > 0.5) if(err > 0.5f)
{ {
std::cerr << "CAUTION: Gross error found at entry " << i << ".\n"; std::cerr << "CAUTION: Gross error found at entry " << i << ".\n";
std::cerr << "Found: " << point << " Expected " << expected << " Error: " << err << std::endl; std::cerr << "Found: " << point << " Expected " << expected << " Error: " << err << std::endl;
@@ -240,7 +240,7 @@ test_result<Real> test_hetero(const A& a, F1 test_func, F2 expect_func)
print_row(row, std::cerr); print_row(row, std::cerr);
BOOST_ERROR("Unexpected non-finite result"); BOOST_ERROR("Unexpected non-finite result");
} }
if(err > 0.5) if(err > 0.5f)
{ {
std::cerr << "CAUTION: Gross error found at entry " << i << ".\n"; std::cerr << "CAUTION: Gross error found at entry " << i << ".\n";
std::cerr << "Found: " << point << " Expected " << expected << " Error: " << err << std::endl; std::cerr << "Found: " << point << " Expected " << expected << " Error: " << err << std::endl;

2
test/Jamfile.v2
View File

@@ -178,6 +178,8 @@ test-suite special_fun :
[ run ccmath_fma_test.cpp /boost/test//boost_unit_test_framework : : : [ requires cxx17_if_constexpr ] $(float128_type) ] [ run ccmath_fma_test.cpp /boost/test//boost_unit_test_framework : : : [ requires cxx17_if_constexpr ] $(float128_type) ]
[ run ccmath_signbit_test.cpp /boost/test//boost_unit_test_framework : : : [ requires cxx17_if_constexpr ] $(float128_type) ] [ run ccmath_signbit_test.cpp /boost/test//boost_unit_test_framework : : : [ requires cxx17_if_constexpr ] $(float128_type) ]
[ run log1p_expm1_test.cpp test_instances//test_instances pch_light /boost/test//boost_unit_test_framework ] [ run log1p_expm1_test.cpp test_instances//test_instances pch_light /boost/test//boost_unit_test_framework ]
[ run log1p_expm1_extra_test.cpp pch_light /boost/test//boost_unit_test_framework ]
[ run log1p_expm1_stdfloat_test.cpp pch_light /boost/test//boost_unit_test_framework ]
[ run powm1_sqrtp1m1_test.cpp test_instances//test_instances pch_light /boost/test//boost_unit_test_framework ] [ run powm1_sqrtp1m1_test.cpp test_instances//test_instances pch_light /boost/test//boost_unit_test_framework ]
[ run git_issue_705.cpp /boost/test//boost_unit_test_framework ] [ run git_issue_705.cpp /boost/test//boost_unit_test_framework ]
[ run git_issue_810.cpp /boost/test//boost_unit_test_framework ] [ run git_issue_810.cpp /boost/test//boost_unit_test_framework ]

73
test/log1p_expm1_extra_test.cpp Normal file
View File

@@ -0,0 +1,73 @@
// Copyright John Maddock 2005.
// Copyright Paul A. Bristow 2010
// Use, modification and distribution are subject to the
// Boost Software License, Version 1.0. (See accompanying file
// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
#include <pch_light.hpp>
#define SC_(x) static_cast<T>(BOOST_STRINGIZE(x))
#define BOOST_MATH_OVERFLOW_ERROR_POLICY ignore_error
#define BOOST_TEST_MAIN
#include <boost/test/unit_test.hpp>
#include <boost/test/tools/floating_point_comparison.hpp>
#include <boost/math/special_functions/log1p.hpp>
#include <boost/math/special_functions/expm1.hpp>
#include "log1p_expm1_test.hpp"
#include <boost/multiprecision/cpp_bin_float.hpp>
//
// DESCRIPTION:
// ~~~~~~~~~~~~
//
// This file tests the functions log1p and expm1. The accuracy tests
// use values generated with NTL::RR at 1000-bit precision
// and our generic versions of these functions.
//
// Note that this tests the generic versions of our functions using
// software emulated types. This ensures these still get tested
// even though we mostly defer to std::expm1 and std::log1p at
// these precisions nowadays
//
void expected_results()
{
//
// Define the max and mean errors expected for
// various compilers and platforms.
//
//
// Catch all cases come last:
//
add_expected_result(
".*", // compiler
".*", // stdlib
".*", // platform
".*", // test type(s)
".*", // test data group
".*", // test function
4, // Max Peek error
3); // Max mean error
//
// Finish off by printing out the compiler/stdlib/platform names,
// we do this to make it easier to mark up expected error rates.
//
std::cout << "Tests run with " << BOOST_COMPILER << ", "
<< BOOST_STDLIB << ", " << BOOST_PLATFORM << std::endl;
}
BOOST_AUTO_TEST_CASE( test_main )
{
expected_results();
BOOST_MATH_CONTROL_FP;
test(boost::multiprecision::cpp_bin_float_single(0), "cpp_bin_float_single");
test(boost::multiprecision::cpp_bin_float_double(0), "cpp_bin_float_double");
test(boost::multiprecision::cpp_bin_float_double_extended(0), "cpp_bin_float_double_extended");
test(boost::multiprecision::cpp_bin_float_quad(0), "cpp_bin_float_quad");
}

87
test/log1p_expm1_stdfloat_test.cpp Normal file
View File

@@ -0,0 +1,87 @@
// Copyright John Maddock 2005.
// Copyright Paul A. Bristow 2010
// Use, modification and distribution are subject to the
// Boost Software License, Version 1.0. (See accompanying file
// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
#include <pch_light.hpp>
#if defined __has_include
# if __cplusplus > 202002L || _MSVC_LANG > 202002L
# if __has_include (<stdfloat>)
# include <stdfloat>
# endif
# endif
#endif
#define SC_(x) static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 128, x))
#define BOOST_MATH_OVERFLOW_ERROR_POLICY ignore_error
#define BOOST_TEST_MAIN
#include <boost/test/unit_test.hpp>
#include <boost/test/tools/floating_point_comparison.hpp>
#include <boost/math/special_functions/log1p.hpp>
#include <boost/math/special_functions/expm1.hpp>
#include <boost/math/tools/big_constant.hpp>
#include "log1p_expm1_test.hpp"
//
// DESCRIPTION:
// ~~~~~~~~~~~~
//
// This file tests the functions log1p and expm1. The accuracy tests
// use values generated with NTL::RR at 1000-bit precision
// and our generic versions of these functions.
//
// Note that this tests the generic versions of our functions using
// software emulated types. This ensures these still get tested
// even though we mostly defer to std::expm1 and std::log1p at
// these precisions nowadays
//
void expected_results()
{
//
// Define the max and mean errors expected for
// various compilers and platforms.
//
//
// Catch all cases come last:
//
add_expected_result(
".*", // compiler
".*", // stdlib
".*", // platform
".*", // test type(s)
".*", // test data group
".*", // test function
4, // Max Peek error
3); // Max mean error
//
// Finish off by printing out the compiler/stdlib/platform names,
// we do this to make it easier to mark up expected error rates.
//
std::cout << "Tests run with " << BOOST_COMPILER << ", "
<< BOOST_STDLIB << ", " << BOOST_PLATFORM << std::endl;
}
BOOST_AUTO_TEST_CASE( test_main )
{
expected_results();
BOOST_MATH_CONTROL_FP;
#ifdef __STDCPP_FLOAT32_T__
test(std::float32_t(0), "std::float32_t");
#endif
#ifdef __STDCPP_FLOAT64_T__
test(std::float64_t(0), "std::float64_t");
#endif
#ifdef __STDCPP_FLOAT128_T__
// As of gcc-13 this does not work as float128_t has no std::math function overloads
//test(std::float128_t(0), "std::float128_t");
#endif
}