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

Correct garbage result in non-central beta and T.

Browse Source 
Fixes github.com/boostorg/math/issues/1308.
This commit is contained in:
jzmaddock
2025-08-25 12:06:57 +01:00
parent c2284b4887
commit dac677cb6f
4 changed files with 45 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

4
include/boost/math/distributions/non_central_beta.hpp
View File

@@ -112,11 +112,15 @@ namespace boost
last_term = term; last_term = term;
} }
last_term = 0; last_term = 0;
T betaf_lim = betaf * tools::epsilon<T>() * 4;
for(auto i = k + 1; ; ++i) for(auto i = k + 1; ; ++i)
{ {
poisf *= l2 / i; poisf *= l2 / i;
xtermf *= (x * (a + b + i - 2)) / (a + i - 1); xtermf *= (x * (a + b + i - 2)) / (a + i - 1);
betaf -= xtermf; betaf -= xtermf;
if (betaf < betaf_lim)
break; // nothing but garbage bits in betaf!!
T term = poisf * betaf; T term = poisf * betaf;
sum += term; sum += term;

4
include/boost/math/distributions/non_central_t.hpp
View File

@@ -100,11 +100,15 @@ namespace boost
++count; ++count;
} }
last_term = 0; last_term = 0;
T betaf_lim = betaf * tools::epsilon<T>() * 4;
for(auto i = k + 1; ; ++i) for(auto i = k + 1; ; ++i)
{ {
poisf *= d2 / (i + 0.5f); poisf *= d2 / (i + 0.5f);
xtermf *= (x * (v / 2 + i - 1)) / (i); xtermf *= (x * (v / 2 + i - 1)) / (i);
betaf -= xtermf; betaf -= xtermf;
if (betaf < betaf_lim)
break; // Nothing but garbage left in betaf now!!
T term = poisf * betaf; T term = poisf * betaf;
sum += term; sum += term;
if((fabs(last_term) >= fabs(term)) && (fabs(term/sum) < errtol)) if((fabs(last_term) >= fabs(term)) && (fabs(term/sum) < errtol))

1
test/Jamfile.v2
View File

@@ -201,6 +201,7 @@ test-suite special_fun :
[ run git_issue_1249.cpp /boost/test//boost_unit_test_framework : : : [ check-target-builds ../config//is_ci_sanitizer_run "Santizer Build" : : <build>no ] ] [ run git_issue_1249.cpp /boost/test//boost_unit_test_framework : : : [ check-target-builds ../config//is_ci_sanitizer_run "Santizer Build" : : <build>no ] ]
[ run git_issue_1255.cpp ] [ run git_issue_1255.cpp ]
[ run git_issue_1247.cpp ] [ run git_issue_1247.cpp ]
[ run git_issue_1308.cpp /boost/test//boost_unit_test_framework ]
[ run special_functions_test.cpp /boost/test//boost_unit_test_framework ] [ run special_functions_test.cpp /boost/test//boost_unit_test_framework ]
[ run test_airy.cpp test_instances//test_instances pch_light /boost/test//boost_unit_test_framework ] [ run test_airy.cpp test_instances//test_instances pch_light /boost/test//boost_unit_test_framework ]
[ run test_bessel_j.cpp test_instances//test_instances pch_light /boost/test//boost_unit_test_framework ] [ run test_bessel_j.cpp test_instances//test_instances pch_light /boost/test//boost_unit_test_framework ]

36
test/git_issue_1308.cpp Normal file
View File

@@ -0,0 +1,36 @@
// Copyright John Maddock 2025.
// 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)
#define BOOST_TEST_MAIN
#include <boost/math/distributions.hpp>
#include <boost/multiprecision/cpp_bin_float.hpp>
#include <iostream>
#include <boost/test/tools/floating_point_comparison.hpp>
#include <boost/test/unit_test.hpp>
typedef boost::math::policies::policy<
boost::math::policies::promote_double<false>
> SciPyPolicy;
using mp_t = boost::multiprecision::cpp_bin_float_quad;
BOOST_AUTO_TEST_CASE(test_main)
{
double v = 980.0;
double l = 38.0;
double x = 1.5;
// With policy
double y_policy = boost::math::cdf(
boost::math::non_central_t_distribution<double, SciPyPolicy>(v, l), x);
double reference = 1.1824454111413493e-291;
std::cout << std::setprecision(17);
std::cout << "With SciPyPolicy: cdf(noncentral_t; x=" << x << ", v=" << v << ", l=" << l << ") = " << y_policy << std::endl;
double tolerance = 1e-8; // precision is limited due to cancellation
BOOST_CHECK_CLOSE_FRACTION(y_policy, reference, tolerance);
}