From dac677cb6f334a6d5fb39f15538e0aa6f32ed76b Mon Sep 17 00:00:00 2001
From: jzmaddock <john@johnmaddock.co.uk>
Date: Mon, 25 Aug 2025 12:06:57 +0100
Subject: [PATCH] Correct garbage result in non-central beta and T. Fixes
 github.com/boostorg/math/issues/1308.

---
 .../math/distributions/non_central_beta.hpp   |  4 +++
 .../math/distributions/non_central_t.hpp      |  4 +++
 test/Jamfile.v2                               |  1 +
 test/git_issue_1308.cpp                       | 36 +++++++++++++++++++
 4 files changed, 45 insertions(+)
 create mode 100644 test/git_issue_1308.cpp

diff --git a/include/boost/math/distributions/non_central_beta.hpp b/include/boost/math/distributions/non_central_beta.hpp
index f311f7773..ea0c5a168 100644
--- a/include/boost/math/distributions/non_central_beta.hpp
+++ b/include/boost/math/distributions/non_central_beta.hpp
@@ -112,11 +112,15 @@ namespace boost
                last_term = term;
             }
             last_term = 0;
+            T betaf_lim = betaf * tools::epsilon<T>() * 4;
             for(auto i = k + 1; ; ++i)
             {
                poisf *= l2 / i;
                xtermf *= (x * (a + b + i - 2)) / (a + i - 1);
                betaf -= xtermf;
+               
+               if (betaf < betaf_lim)
+                  break; // nothing but garbage bits in betaf!!
 
                T term = poisf * betaf;
                sum += term;
diff --git a/include/boost/math/distributions/non_central_t.hpp b/include/boost/math/distributions/non_central_t.hpp
index a6c988d75..85581a4c8 100644
--- a/include/boost/math/distributions/non_central_t.hpp
+++ b/include/boost/math/distributions/non_central_t.hpp
@@ -100,11 +100,15 @@ namespace boost
                ++count;
             }
             last_term = 0;
+            T betaf_lim = betaf * tools::epsilon<T>() * 4;
             for(auto i = k + 1; ; ++i)
             {
                poisf *= d2 / (i + 0.5f);
                xtermf *= (x * (v / 2 + i - 1)) / (i);
                betaf -= xtermf;
+               if (betaf < betaf_lim)
+                  break; // Nothing but garbage left in betaf now!!
+
                T term = poisf * betaf;
                sum += term;
                if((fabs(last_term) >= fabs(term)) && (fabs(term/sum) < errtol))
diff --git a/test/Jamfile.v2 b/test/Jamfile.v2
index 0c43f2084..10f962305 100644
--- a/test/Jamfile.v2
+++ b/test/Jamfile.v2
@@ -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_1255.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 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  ]
diff --git a/test/git_issue_1308.cpp b/test/git_issue_1308.cpp
new file mode 100644
index 000000000..66cc35a7f
--- /dev/null
+++ b/test/git_issue_1308.cpp
@@ -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);
+}
+