boost_mirror/math
1
0
Fork 0
mirror of https://github.com/boostorg/math.git synced 2025-05-11 21:33:52 +00:00
Code Issues Packages Projects Releases Wiki Activity

Improve erf/expm1/expint coverage.

Browse Source 
In the expm1 case, tighten up error handling and testing.
This commit is contained in:
jzmaddock 2024-03-06 18:47:48 +00:00
parent 7fa77fcac6
commit 9ca40b9f01
6 changed files with 129 additions and 135 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
include/boost/math/special_functions/erf.hpp
View File

@ -690,7 +690,7 @@ T erf_imp(T z, bool invert, const Policy& pol, const std::integral_constant<int,
BOOST_MATH_INSTRUMENT_CODE("113-bit precision erf_imp called");
if ((boost::math::isnan)(z))
return policies::raise_domain_error("boost::math::erf<%1%>(%1%)", "Expected a finite argument but got %1%", z, pol);
return policies::raise_domain_error("boost::math::erf<%1%>(%1%)", "Expected a finite argument but got %1%", z, pol); // LCOV_EXCL_LINE checked as covered.
if(z < 0)
{

110
include/boost/math/special_functions/expint.hpp
View File

@ -58,6 +58,7 @@ T expint_1_rational(const T& z, const std::integral_constant<int, 53>&)
// Maximum Deviation Found: 2.006e-18
// Expected Error Term: 2.006e-18
// Max error found at double precision: 2.760e-17
// LCOV_EXCL_START
static const T Y = 0.66373538970947265625F;
static const T P[6] = {
BOOST_MATH_BIG_CONSTANT(T, 53, 0.0865197248079397976498),
@ -75,6 +76,7 @@ T expint_1_rational(const T& z, const std::integral_constant<int, 53>&)
BOOST_MATH_BIG_CONSTANT(T, 53, 0.000131049900798434683324),
BOOST_MATH_BIG_CONSTANT(T, 53, -0.528611029520217142048e-6)
};
// LCOV_EXCL_STOP
result = tools::evaluate_polynomial(P, z)
/ tools::evaluate_polynomial(Q, z);
result += z - log(z) - Y;
@ -83,6 +85,7 @@ T expint_1_rational(const T& z, const std::integral_constant<int, 53>&)
{
// Maximum Deviation Found (interpolated): 1.444e-17
// Max error found at double precision: 3.119e-17
// LCOV_EXCL_START
static const T P[11] = {
BOOST_MATH_BIG_CONSTANT(T, 53, -0.121013190657725568138e-18),
BOOST_MATH_BIG_CONSTANT(T, 53, -0.999999999999998811143),
@ -110,6 +113,7 @@ T expint_1_rational(const T& z, const std::integral_constant<int, 53>&)
BOOST_MATH_BIG_CONSTANT(T, 53, 1229.20784182403048905),
BOOST_MATH_BIG_CONSTANT(T, 53, -0.776491285282330997549)
};
// LCOV_EXCL_STOP
T recip = 1 / z;
result = 1 + tools::evaluate_polynomial(P, recip)
/ tools::evaluate_polynomial(Q, recip);
@ -132,7 +136,7 @@ T expint_1_rational(const T& z, const std::integral_constant<int, 64>&)
// Maximum Deviation Found: 3.807e-20
// Expected Error Term: 3.807e-20
// Max error found at long double precision: 6.249e-20
// LCOV_EXCL_START
static const T Y = 0.66373538970947265625F;
static const T P[6] = {
BOOST_MATH_BIG_CONSTANT(T, 64, 0.0865197248079397956816),
@ -151,6 +155,7 @@ T expint_1_rational(const T& z, const std::integral_constant<int, 64>&)
BOOST_MATH_BIG_CONSTANT(T, 64, -0.202872781770207871975e-5),
BOOST_MATH_BIG_CONSTANT(T, 64, 0.52779248094603709945e-7)
};
// LCOV_EXCL_STOP
result = tools::evaluate_polynomial(P, z)
/ tools::evaluate_polynomial(Q, z);
result += z - log(z) - Y;
@ -159,6 +164,7 @@ T expint_1_rational(const T& z, const std::integral_constant<int, 64>&)
{
// Maximum Deviation Found (interpolated): 2.220e-20
// Max error found at long double precision: 1.346e-19
// LCOV_EXCL_START
static const T P[14] = {
BOOST_MATH_BIG_CONSTANT(T, 64, -0.534401189080684443046e-23),
BOOST_MATH_BIG_CONSTANT(T, 64, -0.999999999999999999905),
@ -191,6 +197,7 @@ T expint_1_rational(const T& z, const std::integral_constant<int, 64>&)
BOOST_MATH_BIG_CONSTANT(T, 64, 73930.2995984054930821),
BOOST_MATH_BIG_CONSTANT(T, 64, 2063.86994219629165937)
};
// LCOV_EXCL_STOP
T recip = 1 / z;
result = 1 + tools::evaluate_polynomial(P, recip)
/ tools::evaluate_polynomial(Q, recip);
@ -213,7 +220,7 @@ T expint_1_rational(const T& z, const std::integral_constant<int, 113>&)
// Maximum Deviation Found: 2.477e-35
// Expected Error Term: 2.477e-35
// Max error found at long double precision: 6.810e-35
// LCOV_EXCL_START
static const T Y = 0.66373538970947265625F;
static const T P[10] = {
BOOST_MATH_BIG_CONSTANT(T, 113, 0.0865197248079397956434879099175975937),
@ -239,6 +246,7 @@ T expint_1_rational(const T& z, const std::integral_constant<int, 113>&)
BOOST_MATH_BIG_CONSTANT(T, 113, 0.369373328141051577845488477377890236e-9),
BOOST_MATH_BIG_CONSTANT(T, 113, -0.274149801370933606409282434677600112e-12)
};
// LCOV_EXCL_STOP
result = tools::evaluate_polynomial(P, z)
/ tools::evaluate_polynomial(Q, z);
result += z - log(z) - Y;
@ -247,7 +255,7 @@ T expint_1_rational(const T& z, const std::integral_constant<int, 113>&)
{
// Max error in interpolated form: 5.614e-35
// Max error found at long double precision: 7.979e-35
// LCOV_EXCL_START
static const T Y = 0.70190334320068359375F;
static const T P[16] = {
@ -286,6 +294,7 @@ T expint_1_rational(const T& z, const std::integral_constant<int, 113>&)
BOOST_MATH_BIG_CONSTANT(T, 113, 169.845369689596739824177412096477219),
BOOST_MATH_BIG_CONSTANT(T, 113, 2.17607292280092201170768401876895354)
};
// LCOV_EXCL_STOP
T recip = 1 / z;
result = Y + tools::evaluate_polynomial(P, recip)
/ tools::evaluate_polynomial(Q, recip);
@ -295,7 +304,7 @@ T expint_1_rational(const T& z, const std::integral_constant<int, 113>&)
{
// Max error in interpolated form: 4.413e-35
// Max error found at long double precision: 8.928e-35
// LCOV_EXCL_START
static const T P[19] = {
BOOST_MATH_BIG_CONSTANT(T, 113, -0.559148411832951463689610809550083986e-40),
BOOST_MATH_BIG_CONSTANT(T, 113, -0.999999999999999999999999999999999997),
@ -339,6 +348,7 @@ T expint_1_rational(const T& z, const std::integral_constant<int, 113>&)
BOOST_MATH_BIG_CONSTANT(T, 113, 70242279152.8241187845178443118302693),
BOOST_MATH_BIG_CONSTANT(T, 113, -37633302.9409263839042721539363416685)
};
// LCOV_EXCL_STOP
T recip = 1 / z;
result = 1 + tools::evaluate_polynomial(P, recip)
/ tools::evaluate_polynomial(Q, recip);
@ -1486,94 +1496,6 @@ T expint_i_imp(T z, const Policy& pol, const std::integral_constant<int, 113>& t
return result;
}
template <class T, class Policy, class tag>
struct expint_i_initializer
{
struct init
{
init()
{
do_init(tag());
}
static void do_init(const std::integral_constant<int, 0>&){}
static void do_init(const std::integral_constant<int, 53>&)
{
boost::math::expint(T(5), Policy());
boost::math::expint(T(7), Policy());
boost::math::expint(T(18), Policy());
boost::math::expint(T(38), Policy());
boost::math::expint(T(45), Policy());
}
static void do_init(const std::integral_constant<int, 64>&)
{
boost::math::expint(T(5), Policy());
boost::math::expint(T(7), Policy());
boost::math::expint(T(18), Policy());
boost::math::expint(T(38), Policy());
boost::math::expint(T(45), Policy());
}
static void do_init(const std::integral_constant<int, 113>&)
{
boost::math::expint(T(5), Policy());
boost::math::expint(T(7), Policy());
boost::math::expint(T(17), Policy());
boost::math::expint(T(25), Policy());
boost::math::expint(T(40), Policy());
boost::math::expint(T(50), Policy());
boost::math::expint(T(80), Policy());
boost::math::expint(T(200), Policy());
boost::math::expint(T(220), Policy());
}
void force_instantiate()const{}
};
static const init initializer;
static void force_instantiate()
{
initializer.force_instantiate();
}
};
template <class T, class Policy, class tag>
const typename expint_i_initializer<T, Policy, tag>::init expint_i_initializer<T, Policy, tag>::initializer;
template <class T, class Policy, class tag>
struct expint_1_initializer
{
struct init
{
init()
{
do_init(tag());
}
static void do_init(const std::integral_constant<int, 0>&){}
static void do_init(const std::integral_constant<int, 53>&)
{
boost::math::expint(1, T(0.5), Policy());
boost::math::expint(1, T(2), Policy());
}
static void do_init(const std::integral_constant<int, 64>&)
{
boost::math::expint(1, T(0.5), Policy());
boost::math::expint(1, T(2), Policy());
}
static void do_init(const std::integral_constant<int, 113>&)
{
boost::math::expint(1, T(0.5), Policy());
boost::math::expint(1, T(2), Policy());
boost::math::expint(1, T(6), Policy());
}
void force_instantiate()const{}
};
static const init initializer;
static void force_instantiate()
{
initializer.force_instantiate();
}
};
template <class T, class Policy, class tag>
const typename expint_1_initializer<T, Policy, tag>::init expint_1_initializer<T, Policy, tag>::initializer;
template <class T, class Policy>
inline typename tools::promote_args<T>::type
expint_forwarder(T z, const Policy& /*pol*/, std::true_type const&)
@ -1594,8 +1516,6 @@ inline typename tools::promote_args<T>::type
precision_type::value <= 113 ? 113 : 0
> tag_type;
expint_i_initializer<value_type, forwarding_policy, tag_type>::force_instantiate();
return policies::checked_narrowing_cast<result_type, forwarding_policy>(detail::expint_i_imp(
static_cast<value_type>(z),
forwarding_policy(),
@ -1631,8 +1551,6 @@ inline typename tools::promote_args<T>::type
precision_type::value <= 113 ? 113 : 0
> tag_type;
detail::expint_1_initializer<value_type, forwarding_policy, tag_type>::force_instantiate();
return policies::checked_narrowing_cast<result_type, forwarding_policy>(detail::expint_imp(
n,
static_cast<value_type>(z),

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

@ -69,37 +69,6 @@ namespace detail
expm1_series& operator=(const expm1_series&) = delete;
};
template <class T, class Policy, class tag>
struct expm1_initializer
{
struct init
{
init()
{
do_init(tag());
}
template <int N>
static void do_init(const std::integral_constant<int, N>&){}
static void do_init(const std::integral_constant<int, 64>&)
{
expm1(T(0.5));
}
static void do_init(const std::integral_constant<int, 113>&)
{
expm1(T(0.5));
}
void force_instantiate()const{}
};
static const init initializer;
static void force_instantiate()
{
initializer.force_instantiate();
}
};
template <class T, class Policy, class tag>
const typename expm1_initializer<T, Policy, tag>::init expm1_initializer<T, Policy, tag>::initializer;
//
// Algorithm expm1 is part of C99, but is not yet provided by many compilers.
//
@ -142,6 +111,10 @@ T expm1_imp(T x, const std::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>())
@ -169,6 +142,10 @@ T expm1_imp(T x, const std::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>())
@ -212,6 +189,10 @@ T expm1_imp(T x, const std::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>())
@ -278,8 +259,6 @@ inline typename tools::promote_args<T>::type expm1(T x, const Policy& /* pol */)
precision_type::value <= 113 ? 113 : 0
> tag_type;
detail::expm1_initializer<value_type, forwarding_policy, tag_type>::force_instantiate();
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%)");
@ -294,14 +273,85 @@ inline typename tools::promote_args<T>::type expm1(T x, const Policy& /* pol */)
#if defined(BOOST_HAS_EXPM1) && !(defined(__osf__) && defined(__DECCXX_VER))
# ifdef BOOST_MATH_USE_C99
inline float expm1(float x, const policies::policy<>&){ return ::expm1f(x); }
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 ::expm1f(x);
}
# ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
inline long double expm1(long double x, const policies::policy<>&){ return ::expm1l(x); }
template <class Policy>
inline long double 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<long double>())
return policies::raise_overflow_error<long double>("boost::math::expm1<%1%>(%1%)", nullptr, Policy());
}
return ::expm1l(x);
}
# endif
# else
inline float expm1(float x, const policies::policy<>&){ return static_cast<float>(::expm1(x)); }
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
inline double expm1(double x, const policies::policy<>&){ return ::expm1(x); }
template <class Policy>
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
template <class T>

4
test/log1p_expm1_test.hpp
View File

@ -89,5 +89,9 @@ void test(T, const char* type_name)
#endif
#endif
}
if (std::numeric_limits<T>::has_quiet_NaN)
{
BOOST_MATH_CHECK_THROW(boost::math::expm1(std::numeric_limits<T>::quiet_NaN()), std::domain_error);
}
}

2
test/test_expint.cpp
View File

@ -79,7 +79,7 @@ void expected_results()
"float|double|long double", // test type(s)
".*Ei.*", // test data group
".*", 6, 3); // test function
if(std::numeric_limits<long double>::digits > 100)
BOOST_IF_CONSTEXPR (std::numeric_limits<long double>::digits > 100)
{
add_expected_result(
".*", // compiler

22
test/test_expint.hpp
View File

@ -190,5 +190,27 @@ void test_spots(T, const char* t)
BOOST_CHECK_CLOSE(::boost::math::expint(static_cast<T>(-0.5)), static_cast<T>(-0.559773594776160811746795939315085235226846890316353515248293L), tolerance);
BOOST_CHECK_CLOSE(::boost::math::expint(static_cast<T>(-1)), static_cast<T>(-0.219383934395520273677163775460121649031047293406908207577979L), tolerance);
BOOST_CHECK_CLOSE(::boost::math::expint(static_cast<T>(-50.5)), static_cast<T>(-2.27237132932219350440719707268817831250090574830769670186618e-24L), tolerance);
//
// Extra coverage cases:
//
BOOST_CHECK_EQUAL(boost::math::expint(-boost::math::tools::max_value<T>()), T(0));
BOOST_CHECK_THROW(boost::math::expint(1, T(-1)), std::domain_error);
BOOST_CHECK_THROW(boost::math::expint(2, T(-1)), std::domain_error);
BOOST_CHECK_EQUAL(boost::math::expint(2, T(0)), T(1));
BOOST_CHECK_EQUAL(boost::math::expint(3, T(0)), T(0.5));
BOOST_IF_CONSTEXPR(std::numeric_limits<T>::has_infinity)
{
BOOST_CHECK_EQUAL(boost::math::expint(1, T(0)), std::numeric_limits<T>::infinity());
BOOST_CHECK_EQUAL(boost::math::expint(T(0)), -std::numeric_limits<T>::infinity());
BOOST_CHECK_EQUAL(boost::math::expint(boost::math::tools::log_max_value<T>() * 2), std::numeric_limits<T>::infinity());
BOOST_CHECK_EQUAL(boost::math::expint(boost::math::tools::log_max_value<T>() + T(38)), std::numeric_limits<T>::infinity());
}
else
{
BOOST_CHECK_GE(boost::math::expint(1, T(0)), boost::math::tools::max_value<T>());
BOOST_CHECK_LE(boost::math::expint(T(0)), -boost::math::tools::max_value<T>());
BOOST_CHECK_GE(boost::math::expint(boost::math::tools::log_max_value<T>() * 2), boost::math::tools::max_value<T>());
BOOST_CHECK_GE(boost::math::expint(boost::math::tools::log_max_value<T>() + T(38)), boost::math::tools::max_value<T>());
}
}