diff --git a/include/boost/math/distributions/students_t.hpp b/include/boost/math/distributions/students_t.hpp index bc1c7139eb..69e5dde43c 100644 --- a/include/boost/math/distributions/students_t.hpp +++ b/include/boost/math/distributions/students_t.hpp @@ -129,6 +129,11 @@ BOOST_MATH_GPU_ENABLED inline RealType pdf(const students_t_distribution sqrt(tools::max_value())) + { + // Exact limit: the pdf underflows to zero long before x * x can overflow. + return 0; + } RealType basem1 = x * x / df; if(basem1 < 0.125) { @@ -146,6 +151,8 @@ BOOST_MATH_GPU_ENABLED inline RealType pdf(const students_t_distribution BOOST_MATH_GPU_ENABLED inline RealType cdf(const students_t_distribution& dist, const RealType& x) { + BOOST_MATH_STD_USING // for ADL of std functions + RealType error_result; // degrees_of_freedom > 0 or infinity check: RealType df = dist.degrees_of_freedom(); @@ -201,6 +208,11 @@ BOOST_MATH_GPU_ENABLED inline RealType cdf(const students_t_distribution sqrt(tools::max_value())) + { + // Exact tail values: avoids a spurious intermediate overflow in x * x below. + return (x < 0) ? static_cast(0) : static_cast(1); + } RealType x2 = x * x; RealType probability; if(df > 2 * x2) diff --git a/test/Jamfile.v2 b/test/Jamfile.v2 index c5bde01635..5b67e474df 100644 --- a/test/Jamfile.v2 +++ b/test/Jamfile.v2 @@ -895,6 +895,7 @@ test-suite distribution_tests : : test_poisson_real_concept ] [ run test_rayleigh.cpp /boost/test//boost_unit_test_framework ] [ run test_students_t.cpp /boost/test//boost_unit_test_framework ] + [ run test_students_t_overflow.cpp /boost/test//boost_unit_test_framework ] [ run test_skew_normal.cpp /boost/test//boost_unit_test_framework ] [ run test_triangular.cpp pch /boost/test//boost_unit_test_framework ] [ run test_uniform.cpp pch /boost/test//boost_unit_test_framework ] diff --git a/test/test_students_t_overflow.cpp b/test/test_students_t_overflow.cpp new file mode 100644 index 0000000000..d4020d3bf4 --- /dev/null +++ b/test/test_students_t_overflow.cpp @@ -0,0 +1,133 @@ +// Copyright Anton Leontev 2026. +// 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) +// +// Regression test for spurious intermediate overflow of x*x in the +// Student's t pdf() and cdf(): for |x| > sqrt(max_value) the functions +// must return the exact tail values (0 for the pdf, 0/1 for the cdf and +// its complement) without raising FE_OVERFLOW and without invoking any +// error handler. + +#define BOOST_TEST_MAIN +#include +#include +#include + +#include +#include +#include + +#ifndef BOOST_MATH_STANDALONE +#include // for real_concept +#include // for cpp_bin_float_50 +#endif + +#include +#include +#include + +#pragma STDC FENV_ACCESS ON + +template +RealType overflowing_deviate() +{ + BOOST_MATH_STD_USING // ADL sqrt for built-in and UDT types + using boost::math::tools::max_value; + return 2 * sqrt(max_value()); +} + +// For |x| > sqrt(max_value) the mathematically exact results have already +// saturated: pdf == 0, cdf == 0 or 1. These must be returned directly. +// A throwing policy is used so that any error-handler invocation would +// fail the test with an unexpected exception. +template +void test_tail_values(RealType) +{ + using namespace boost::math; + using namespace boost::math::policies; + + typedef policy > strict_policy; + typedef students_t_distribution dist_t; + + dist_t dist(static_cast(5)); + const RealType big = overflowing_deviate(); + + BOOST_CHECK((boost::math::isfinite)(big)); + + BOOST_CHECK_EQUAL(pdf(dist, big), static_cast(0)); + BOOST_CHECK_EQUAL(pdf(dist, -big), static_cast(0)); + + BOOST_CHECK_EQUAL(cdf(dist, big), static_cast(1)); + BOOST_CHECK_EQUAL(cdf(dist, -big), static_cast(0)); + + BOOST_CHECK_EQUAL(cdf(complement(dist, big)), static_cast(0)); + BOOST_CHECK_EQUAL(cdf(complement(dist, -big)), static_cast(1)); +} + +// Built-in types: the calls must not leave FE_OVERFLOW set -- avoiding the +// spurious intermediate overflow is the whole point of the guard. +template +void test_no_fp_exceptions(RealType) +{ + using namespace boost::math; + + students_t_distribution dist(static_cast(5)); + const RealType big = overflowing_deviate(); + + std::feclearexcept(FE_ALL_EXCEPT); + RealType p = pdf(dist, big); + RealType c1 = cdf(dist, big); + RealType c2 = cdf(dist, -big); + BOOST_CHECK_EQUAL(std::fetestexcept(FE_OVERFLOW), 0); + BOOST_CHECK_EQUAL(p, static_cast(0)); + BOOST_CHECK_EQUAL(c1, static_cast(1)); + BOOST_CHECK_EQUAL(c2, static_cast(0)); +} + +// Ordinary arguments must be unaffected by the guard. +template +void test_no_regression(RealType) +{ + using namespace boost::math; + + students_t_distribution dist(static_cast(10)); + const RealType tol = boost::math::tools::epsilon() * 100; + + BOOST_CHECK_CLOSE_FRACTION(cdf(dist, static_cast(0)), + static_cast(0.5), tol); + + RealType p = pdf(dist, static_cast(1.5)); + BOOST_CHECK((boost::math::isfinite)(p)); + BOOST_CHECK(p > 0); + + RealType big_but_safe = static_cast(1e6); + BOOST_CHECK_CLOSE_FRACTION(cdf(dist, big_but_safe), + static_cast(1), tol); +} + +BOOST_AUTO_TEST_CASE(students_t_overflow_test) +{ + test_tail_values(0.0F); + test_tail_values(0.0); + test_tail_values(0.0L); + + test_no_fp_exceptions(0.0F); + test_no_fp_exceptions(0.0); + test_no_fp_exceptions(0.0L); + + test_no_regression(0.0F); + test_no_regression(0.0); + test_no_regression(0.0L); + + // real_concept and multiprecision compute distribution constants via + // lexical_cast, which standalone mode disables, so both are skipped there. +#ifndef BOOST_MATH_STANDALONE + test_tail_values(boost::math::concepts::real_concept(0)); + test_no_regression(boost::math::concepts::real_concept(0)); + + test_tail_values(boost::multiprecision::cpp_bin_float_50(0)); + test_no_regression(boost::multiprecision::cpp_bin_float_50(0)); +#endif +}