From cf12642591b97256a074914cd4782875316b5025 Mon Sep 17 00:00:00 2001 From: ckormanyos Date: Fri, 29 Aug 2025 09:52:01 +0200 Subject: [PATCH 1/3] Repair neg-zero and tiny-arg erf and erfc --- include/boost/math/special_functions/erf.hpp | 81 ++++++++++++-------- 1 file changed, 50 insertions(+), 31 deletions(-) diff --git a/include/boost/math/special_functions/erf.hpp b/include/boost/math/special_functions/erf.hpp index c1d9061db..fbc80f53b 100644 --- a/include/boost/math/special_functions/erf.hpp +++ b/include/boost/math/special_functions/erf.hpp @@ -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()) + { + // Series[Erf[x], {x, 0, 4}] + // Series[Erfc[x], {x, 0, 4}] + + const T term2 { 2 * z / constants::root_pi() }; + + 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); @@ -161,43 +173,50 @@ T erf_imp(T z, bool invert, const Policy& pol, const Tag& t) return 1 + erf_imp(T(-z), false, pol, t); } - T result; + if (z > 0) + { + T result; - if(!invert && (z > detail::erf_asymptotic_limit())) - { - detail::erf_asympt_series_t s(z); - std::uintmax_t max_iter = policies::get_max_series_iterations(); - result = boost::math::tools::sum_series(s, policies::get_epsilon(), max_iter, 1); - policies::check_series_iterations("boost::math::erf<%1%>(%1%, %1%)", max_iter, pol); - } - else - { - T x = z * z; - if(z < 1.3f) + if(!invert && (z > detail::erf_asymptotic_limit())) { - // 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... - result = erf_series_near_zero_sum(z, pol); - } - else if(x > 1 / tools::epsilon()) - { - // http://functions.wolfram.com/06.27.06.0006.02 - invert = !invert; - result = exp(-x) / (constants::root_pi() * z); + detail::erf_asympt_series_t s(z); + std::uintmax_t max_iter = policies::get_max_series_iterations(); + result = boost::math::tools::sum_series(s, policies::get_epsilon(), max_iter, 1); + policies::check_series_iterations("boost::math::erf<%1%>(%1%, %1%)", max_iter, pol); } else { - // Compute Q: - invert = !invert; - result = z * exp(-x); - result /= boost::math::constants::root_pi(); - result *= upper_gamma_fraction(T(0.5f), x, policies::get_epsilon()); + T x = 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... + result = erf_series_near_zero_sum(z, pol); + } + else if(x > 1 / tools::epsilon()) + { + // http://functions.wolfram.com/06.27.06.0006.02 + invert = !invert; + result = exp(-x) / (constants::root_pi() * z); + } + else + { + // Compute Q: + invert = !invert; + result = z * exp(-x); + result /= boost::math::constants::root_pi(); + result *= upper_gamma_fraction(T(0.5f), x, policies::get_epsilon()); + } } + + return ((!invert) ? result : 1 - result); + } + else + { + // Handle z == 0. + return ((!invert) ? T(0) : T(1)); } - if(invert) - result = 1 - result; - return result; } // LCOV_EXCL_STOP From 739f483394ec5dff5003fc3a2736a503a3835e16 Mon Sep 17 00:00:00 2001 From: ckormanyos Date: Fri, 29 Aug 2025 11:54:51 +0200 Subject: [PATCH 2/3] Further simplify erf and add tests --- include/boost/math/special_functions/erf.hpp | 75 +++++++++----------- test/test_erf.hpp | 14 ++++ 2 files changed, 48 insertions(+), 41 deletions(-) diff --git a/include/boost/math/special_functions/erf.hpp b/include/boost/math/special_functions/erf.hpp index fbc80f53b..dffb0c59c 100644 --- a/include/boost/math/special_functions/erf.hpp +++ b/include/boost/math/special_functions/erf.hpp @@ -153,12 +153,12 @@ 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(fabs(z) < tools::root_epsilon()) + if (fabs(z) < tools::root_epsilon()) { // Series[Erf[x], {x, 0, 4}] // Series[Erfc[x], {x, 0, 4}] - const T term2 { 2 * z / constants::root_pi() }; + const T term2 { z * 2 / constants::root_pi() }; return invert ? 1 - term2 : term2; } @@ -173,50 +173,43 @@ T erf_imp(T z, bool invert, const Policy& pol, const Tag& t) return 1 + erf_imp(T(-z), false, pol, t); } - if (z > 0) + T result; + + if(!invert && (z > detail::erf_asymptotic_limit())) { - T result; - - if(!invert && (z > detail::erf_asymptotic_limit())) - { - detail::erf_asympt_series_t s(z); - std::uintmax_t max_iter = policies::get_max_series_iterations(); - result = boost::math::tools::sum_series(s, policies::get_epsilon(), max_iter, 1); - policies::check_series_iterations("boost::math::erf<%1%>(%1%, %1%)", max_iter, pol); - } - else - { - T x = 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... - result = erf_series_near_zero_sum(z, pol); - } - else if(x > 1 / tools::epsilon()) - { - // http://functions.wolfram.com/06.27.06.0006.02 - invert = !invert; - result = exp(-x) / (constants::root_pi() * z); - } - else - { - // Compute Q: - invert = !invert; - result = z * exp(-x); - result /= boost::math::constants::root_pi(); - result *= upper_gamma_fraction(T(0.5f), x, policies::get_epsilon()); - } - } - - return ((!invert) ? result : 1 - result); + detail::erf_asympt_series_t s(z); + std::uintmax_t max_iter = policies::get_max_series_iterations(); + result = boost::math::tools::sum_series(s, policies::get_epsilon(), max_iter, 1); + policies::check_series_iterations("boost::math::erf<%1%>(%1%, %1%)", max_iter, pol); } else { - // Handle z == 0. - return ((!invert) ? T(0) : T(1)); + 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. Regarding performance, this is way quicker than anything else... + result = erf_series_near_zero_sum(z, pol); + } + else if(z_sq > 1 / tools::epsilon()) + { + // http://functions.wolfram.com/06.27.06.0006.02 + invert = !invert; + result = exp(-z_sq) / (constants::root_pi() * z); + } + else + { + // Compute Q: + invert = !invert; + result = z * exp(-z_sq); + result /= boost::math::constants::root_pi(); + result *= upper_gamma_fraction(T(0.5f), z_sq, policies::get_epsilon()); + } } + + return ((!invert) ? result : 1 - result); } // LCOV_EXCL_STOP diff --git a/test/test_erf.hpp b/test/test_erf.hpp index b70c73953..f242250ad 100644 --- a/test/test_erf.hpp +++ b/test/test_erf.hpp @@ -174,6 +174,20 @@ void test_erf(T, const char* name) do_test_erfc_inv(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() / 8) > T(0)); + BOOST_CHECK(boost::math::erf(boost::math::tools::epsilon() / 32) >= T(0)); + BOOST_CHECK(boost::math::erfc(boost::math::tools::epsilon() / 32) <= T(1)); + + const bool has_negative_zero { (boost::math::signbit)(-T(0)) }; + + if(has_negative_zero) + { + BOOST_CHECK_EQUAL(boost::math::erf(-T(0)), -T(0)); + } + BOOST_IF_CONSTEXPR(std::numeric_limits::has_quiet_NaN) { BOOST_CHECK_THROW(boost::math::erf(std::numeric_limits::quiet_NaN()), std::domain_error); From c1d2becb56cd40d9f1e3c1f6771ab93e7b1f2265 Mon Sep 17 00:00:00 2001 From: ckormanyos Date: Fri, 29 Aug 2025 12:31:18 +0200 Subject: [PATCH 3/3] Repair erroneous implicit cast --- test/test_erf.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_erf.hpp b/test/test_erf.hpp index f242250ad..2147938aa 100644 --- a/test/test_erf.hpp +++ b/test/test_erf.hpp @@ -181,7 +181,7 @@ void test_erf(T, const char* name) BOOST_CHECK(boost::math::erf(boost::math::tools::epsilon() / 32) >= T(0)); BOOST_CHECK(boost::math::erfc(boost::math::tools::epsilon() / 32) <= T(1)); - const bool has_negative_zero { (boost::math::signbit)(-T(0)) }; + const bool has_negative_zero { ((boost::math::signbit)(-T(0)) != 0) }; if(has_negative_zero) {