diff --git a/include/boost/math/distributions/students_t.hpp b/include/boost/math/distributions/students_t.hpp index d768175b6c..02309bfb26 100644 --- a/include/boost/math/distributions/students_t.hpp +++ b/include/boost/math/distributions/students_t.hpp @@ -466,7 +466,7 @@ BOOST_MATH_GPU_ENABLED RealType calculate_hill_df_guess( RealType za = quantile(complement(n, alpha)); RealType zb = quantile(complement(n, beta)); - RealType v_norm = pow((za + zb) * sd / difference_from_mean, 2) - 1; + RealType v_norm = static_cast(pow((za + zb) * sd / difference_from_mean, 2) - 1); RealType hint = v_norm + (za * za - za * zb + zb * zb + 1) / 2; return (hint <= 0) ? RealType(1) : hint; diff --git a/include/boost/math/special_functions/detail/t_distribution_inv.hpp b/include/boost/math/special_functions/detail/t_distribution_inv.hpp index 4b55b4fe21..281f4e1816 100644 --- a/include/boost/math/special_functions/detail/t_distribution_inv.hpp +++ b/include/boost/math/special_functions/detail/t_distribution_inv.hpp @@ -509,7 +509,14 @@ BOOST_MATH_GPU_ENABLED T fast_students_t_quantile_imp(T df, T p, const Policy& p // Get cdf from incomplete beta result: T p0 = f0 / 2 - p; // Get pdf from derivative: - T p1 = f1 * sqrt(y * xb * xb * xb / df); + T rxb = sqrt(xb); + T p1 = f1 * sqrt(y / df); + p1 *= rxb; + p1 *= rxb; + p1 *= rxb; + // If p1 has underflowed, all subsequent calculations will fail: + if (p1 < tools::min_value()) + return t; // // Second derivative divided by p1: // diff --git a/test/test_students_t.cpp b/test/test_students_t.cpp index dfc0921bb2..fe5aaa423b 100644 --- a/test/test_students_t.cpp +++ b/test/test_students_t.cpp @@ -422,6 +422,24 @@ void test_spots(RealType) BOOST_CHECK_THROW(boost::math::quantile(students_t_distribution((std::numeric_limits::min)() / 2), static_cast(0.0025f)), std::overflow_error); } BOOST_CHECK_CLOSE(boost::math::cdf(students_t_distribution(static_cast(1)), ldexp(RealType(1), -30)), static_cast(0.5000000002964491827262478614651948102240210317156775562512071908L), tolerance); + // https://github.com/boostorg/math/issues/1436 + // We have quite limited precision in this area: + if (boost::math::tools::digits() > boost::math::tools::digits()) + { + using std::ldexp; + BOOST_CHECK_CLOSE_FRACTION(boost::math::quantile(students_t_distribution(static_cast(3)), ldexp(RealType(1), -800)), static_cast(-1.94452005735447553162080266906e+80L), 1e-8); + BOOST_CHECK_CLOSE_FRACTION(boost::math::quantile(students_t_distribution(static_cast(3)), ldexp(RealType(1), -1000)), static_cast(-2.27760708321880957407749124767e+100L), 1e-8); + BOOST_CHECK_CLOSE_FRACTION(boost::math::quantile(students_t_distribution(static_cast(5)), ldexp(RealType(1), -1000)), static_cast(-2.52030967154944272146770707689e+60L), 1e-8); + BOOST_CHECK_CLOSE_FRACTION(boost::math::quantile(students_t_distribution(static_cast(7)), ldexp(RealType(1), -1000)), static_cast(-2.02884419021075115948274137329e+43L), 1e-8); + BOOST_CHECK_CLOSE_FRACTION(boost::math::quantile(students_t_distribution(static_cast(10)), ldexp(RealType(1), -1000)), static_cast(-3.25092256692541451558538754799e+30L), 1e-8); + BOOST_CHECK_CLOSE_FRACTION(boost::math::quantile(students_t_distribution(static_cast(30)), ldexp(RealType(1), -1000)), static_cast(-54306462721.8246714347211540407L), 1e-8); + using nopromote = boost::math::policies::policy>; + BOOST_CHECK_CLOSE_FRACTION(boost::math::quantile(students_t_distribution(static_cast(30)), ldexp(RealType(1), -1000)), static_cast(-54306462721.8246714347211540407L), 1e-8); + if (boost::math::tools::digits() > boost::math::tools::digits()) + { + BOOST_CHECK_CLOSE_FRACTION(boost::math::quantile(students_t_distribution(static_cast(30)), ldexp(RealType(1), -16000)), static_cast(-1.77766265021682828433490865393e+161L), 1e-8); + } + } } // Student's t pdf tests.