From eaf876c81ec37e8881f3ea354b93c30357ea2a7c Mon Sep 17 00:00:00 2001
From: jzmaddock <john@johnmaddock.co.uk>
Date: Fri, 13 Oct 2023 18:45:03 +0100
Subject: [PATCH] Correct non-central-t series convergence bug. Fixes
 https://github.com/boostorg/math/issues/1035. See also
 https://github.com/scipy/scipy/issues/19348. Accuracy in left tail is still
 poor, and the reflection formula appears to be to blame as it's use causes
 the series to cancel out the first term, but it appears we have no real
 choice in the matter here.  At least we do now get a few digits correct.

---
 include/boost/math/distributions/non_central_t.hpp | 13 ++++++++++++-
 test/test_nc_t.hpp                                 |  9 +++++++++
 2 files changed, 21 insertions(+), 1 deletion(-)

diff --git a/include/boost/math/distributions/non_central_t.hpp b/include/boost/math/distributions/non_central_t.hpp
index cb9af7e2e5..48483f0cae 100644
--- a/include/boost/math/distributions/non_central_t.hpp
+++ b/include/boost/math/distributions/non_central_t.hpp
@@ -396,10 +396,12 @@ namespace boost
             pois = gamma_p_derivative(T(k+1), d2, pol)
                * tgamma_delta_ratio(T(k + 1), T(0.5f))
                * delta / constants::root_two<T>();
+            BOOST_MATH_INSTRUMENT_VARIABLE(pois);
             // Starting beta term:
             xterm = x < y
                ? ibeta_derivative(T(k + 1), n / 2, x, pol)
                : ibeta_derivative(n / 2, T(k + 1), y, pol);
+            BOOST_MATH_INSTRUMENT_VARIABLE(xterm);
             T poisf(pois), xtermf(xterm);
             T sum = init_val;
             if((pois == 0) || (xterm == 0))
@@ -410,12 +412,16 @@ namespace boost
             // direction for recursion:
             //
             std::uintmax_t count = 0;
+            T old_ratio = 1;  // arbitrary "large" value
             for(auto i = k; i >= 0; --i)
             {
                T term = xterm * pois;
                sum += term;
-               if(((fabs(term/sum) < errtol) && (i != k)) || (term == 0))
+               BOOST_MATH_INSTRUMENT_VARIABLE(sum);
+               T ratio = fabs(term / sum);
+               if(((ratio < errtol) && (i != k) && (ratio < old_ratio)) || (term == 0))
                   break;
+               old_ratio = ratio;
                pois *= (i + 0.5f) / d2;
                xterm *= (i) / (x * (n / 2 + i));
                ++count;
@@ -426,6 +432,7 @@ namespace boost
                      "Series did not converge, closest value was %1%", sum, pol);
                }
             }
+            BOOST_MATH_INSTRUMENT_VARIABLE(sum);
             for(auto i = k + 1; ; ++i)
             {
                poisf *= d2 / (i + 0.5f);
@@ -442,6 +449,7 @@ namespace boost
                      "Series did not converge, closest value was %1%", sum, pol);
                }
             }
+            BOOST_MATH_INSTRUMENT_VARIABLE(sum);
             return sum;
          }
 
@@ -505,9 +513,12 @@ namespace boost
             // Calculate pdf:
             //
             T dt = n * t / (n * n + 2 * n * t * t + t * t * t * t);
+            BOOST_MATH_INSTRUMENT_VARIABLE(dt);
             T result = non_central_beta_pdf(a, b, d2, x, y, pol);
+            BOOST_MATH_INSTRUMENT_VARIABLE(result);
             T tol = tools::epsilon<T>() * result * 500;
             result = non_central_t2_pdf(n, delta, x, y, pol, result);
+            BOOST_MATH_INSTRUMENT_VARIABLE(result);
             if(result <= tol)
                result = 0;
             result *= dt;
diff --git a/test/test_nc_t.hpp b/test/test_nc_t.hpp
index 8fb6accdf3..b04b10ba6f 100644
--- a/test/test_nc_t.hpp
+++ b/test/test_nc_t.hpp
@@ -320,6 +320,15 @@ void test_spots(RealType)
       BOOST_CHECK(boost::math::isnan(pdf(d2, 0.5)));
       BOOST_CHECK(boost::math::isnan(cdf(d2, 0.5)));
    }
+
+   // Bug cases, 
+   // https://github.com/scipy/scipy/issues/19348
+   //
+   {
+      distro1 d(8.0f, 8.5f);
+      BOOST_CHECK_CLOSE(pdf(d, -1), static_cast<RealType>(6.1747948083757028903541988987716621647020752431287e-20), 2e-5);  // Can we do better on accuracy here?
+   }
+
 } // template <class RealType>void test_spots(RealType)
 
 template <class T>