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

Merge pull request #1315 from boostorg/edges_c99_funcs

Browse Source 
Repair neg-zero and tiny-arg erf and erfc
This commit is contained in:
Christopher Kormanyos
2025-08-29 14:05:57 +02:00
committed by GitHub
parent 0992ba6792 c1d2becb56
commit ab265d474f
2 changed files with 36 additions and 10 deletions
Download Patch File Download Diff File Expand all files Collapse all files

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

@@ -153,7 +153,19 @@ T erf_imp(T z, bool invert, const Policy& pol, const Tag& t)
if ((boost::math::isnan)(z))
return policies::raise_domain_error("boost::math::erf<%1%>(%1%)", "Expected a finite argument but got %1%", z, pol);
if(z < 0)
if (fabs(z) < tools::root_epsilon<T>())
{
// Series[Erf[x], {x, 0, 4}]
// Series[Erfc[x], {x, 0, 4}]
const T term2 { z * 2 / constants::root_pi<T>() };
return invert ? 1 - term2 : term2;
}
const bool signbit_result = ((boost::math::signbit)(z) != 0);
if (signbit_result)
{
if(!invert)
return -erf_imp(T(-z), invert, pol, t);
@@ -172,32 +184,32 @@ T erf_imp(T z, bool invert, const Policy& pol, const Tag& t)
}
else
{
T x = z * z;
const T z_sq { z * z };
if(z < 1.3f)
{
// Compute P:
// This is actually good for z p to 2 or so, but the cutoff given seems
// to be the best compromise. Performance wise, this is way quicker than anything else...
// to be the best compromise. Regarding performance, this is way quicker than anything else...
result = erf_series_near_zero_sum(z, pol);
}
else if(x > 1 / tools::epsilon<T>())
else if(z_sq > 1 / tools::epsilon<T>())
{
// http://functions.wolfram.com/06.27.06.0006.02
invert = !invert;
result = exp(-x) / (constants::root_pi<T>() * z);
result = exp(-z_sq) / (constants::root_pi<T>() * z);
}
else
{
// Compute Q:
invert = !invert;
result = z * exp(-x);
result = z * exp(-z_sq);
result /= boost::math::constants::root_pi<T>();
result *= upper_gamma_fraction(T(0.5f), x, policies::get_epsilon<T, Policy>());
result *= upper_gamma_fraction(T(0.5f), z_sq, policies::get_epsilon<T, Policy>());
}
}
if(invert)
result = 1 - result;
return result;
return ((!invert) ? result : 1 - result);
}
// LCOV_EXCL_STOP

14
test/test_erf.hpp
View File

@@ -174,6 +174,20 @@ void test_erf(T, const char* name)
do_test_erfc_inv<T>(erfc_inv_big_data, name, "Inverse Erfc Function: extreme values");
}
BOOST_CHECK_EQUAL(boost::math::erf(T(0)), T(0));
BOOST_CHECK_EQUAL(boost::math::erfc(T(0)), T(1));
BOOST_CHECK(boost::math::erf(boost::math::tools::root_epsilon<T>() / 8) > T(0));
BOOST_CHECK(boost::math::erf(boost::math::tools::epsilon<T>() / 32) >= T(0));
BOOST_CHECK(boost::math::erfc(boost::math::tools::epsilon<T>() / 32) <= T(1));
const bool has_negative_zero { ((boost::math::signbit)(-T(0)) != 0) };
if(has_negative_zero)
{
BOOST_CHECK_EQUAL(boost::math::erf(-T(0)), -T(0));
}
BOOST_IF_CONSTEXPR(std::numeric_limits<T>::has_quiet_NaN)
{
BOOST_CHECK_THROW(boost::math::erf(std::numeric_limits<T>::quiet_NaN()), std::domain_error);