diff --git a/include/boost/math/distributions/non_central_t.hpp b/include/boost/math/distributions/non_central_t.hpp index cc33128c77..779652f95b 100644 --- a/include/boost/math/distributions/non_central_t.hpp +++ b/include/boost/math/distributions/non_central_t.hpp @@ -74,7 +74,7 @@ namespace boost T term = beta * pois; sum += term; // Don't terminate on first term in case we "fixed" k above: - if((fabs(last_term) > fabs(term)) && fabs(term/sum) < errtol) + if(((fabs(last_term) > fabs(term)) && fabs(term/sum) < errtol) || (v == 2 && i == 0)) break; last_term = term; pois *= (i + 0.5f) / d2; @@ -183,7 +183,8 @@ namespace boost term += beta * pois; pois *= (j + 0.5f) / d2; beta -= xterm; - xterm *= (j) / (x * (v / 2 + j - 1)); + if(!(v == 2 && j == 0)) + xterm *= (j) / (x * (v / 2 + j - 1)); } sum += term; diff --git a/test/Jamfile.v2 b/test/Jamfile.v2 index e2c9c09de1..f9ed071285 100644 --- a/test/Jamfile.v2 +++ b/test/Jamfile.v2 @@ -892,6 +892,7 @@ test-suite distribution_tests : [ compile test_dist_deduction_guides.cpp : [ requires cpp_deduction_guides cpp_variadic_templates ] ] [ run git_issue_800.cpp ../../test/build//boost_unit_test_framework ] [ run git_issue_845.cpp ../../test/build//boost_unit_test_framework ] + [ run scipy_issue_14901.cpp ../../test/build//boost_unit_test_framework ] [ run scipy_issue_17146.cpp ../../test/build//boost_unit_test_framework ] [ run scipy_issue_17388.cpp ../../test/build//boost_unit_test_framework ] ; diff --git a/test/scipy_issue_14901.cpp b/test/scipy_issue_14901.cpp new file mode 100644 index 0000000000..9423f78644 --- /dev/null +++ b/test/scipy_issue_14901.cpp @@ -0,0 +1,54 @@ +// Copyright Matt Borland, 2022 +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. (See accompanying file +// LICENSE or copy at http://www.boost.org/LICENSE_1_0.txt) +// +// See: https://github.com/scipy/scipy/issues/14901 + +#include +#include +#include +#include "math_unit_test.hpp" + +#pragma STDC FENV_ACCESS ON + +int main() +{ + auto nct1 = boost::math::non_central_t(2, 2); + const auto cdf1 = boost::math::cdf(nct1, 0.05); + + CHECK_ULP_CLOSE(cdf1, 0.02528206132724582, 20); + + if (std::fetestexcept(FE_INVALID) || std::fetestexcept(FE_DIVBYZERO)) + { + return 1; + } + + auto nct2 = boost::math::non_central_t(1, 3); + const auto cdf2 = boost::math::cdf(nct2, 0.05); + + CHECK_ULP_CLOSE(cdf2, 0.00154456589169420, 20); + + if (std::fetestexcept(FE_INVALID) || std::fetestexcept(FE_DIVBYZERO)) + { + return 2; + } + + const auto pdf1 = boost::math::pdf(nct1, 0.05); + CHECK_ULP_CLOSE(pdf1, 0.05352074159921797, 20); + + if (std::fetestexcept(FE_INVALID) || std::fetestexcept(FE_DIVBYZERO)) + { + return 3; + } + + const auto quantile1 = boost::math::quantile(nct1, 0.5); + CHECK_ULP_CLOSE(quantile1, 2.33039025947198741, 20); + + if (std::fetestexcept(FE_INVALID) || std::fetestexcept(FE_DIVBYZERO)) + { + return 4; + } + + return boost::math::test::report_errors(); +}