Fix for Non-Central T Distribution FE_INVALID

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;

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  ]
 ;

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 <cfenv>
+#include <iostream>
+#include <boost/math/distributions/non_central_t.hpp>
+#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();
+}