diff --git a/include/boost/math/distributions/non_central_beta.hpp b/include/boost/math/distributions/non_central_beta.hpp
index 238297aba0..9d1eedecca 100644
--- a/include/boost/math/distributions/non_central_beta.hpp
+++ b/include/boost/math/distributions/non_central_beta.hpp
@@ -17,6 +17,7 @@
 #include <boost/math/distributions/detail/generic_mode.hpp>
 #include <boost/math/distributions/detail/common_error_handling.hpp> // error checks
 #include <boost/math/special_functions/fpclassify.hpp> // isnan.
+#include <boost/math/special_functions/trunc.hpp>
 #include <boost/math/tools/roots.hpp> // for root finding.
 #include <boost/math/tools/series.hpp>
 
@@ -48,7 +49,7 @@ namespace boost
             // k to zero, when l2 is small, as forward iteration
             // is unstable:
             //
-            int k = itrunc(l2);
+            long long k = lltrunc(l2);
             if(k == 0)
                k = 1;
                // Starting Poisson weight:
@@ -75,7 +76,7 @@ namespace boost
             //
             T last_term = 0;
             std::uintmax_t count = k;
-            for(int i = k; i >= 0; --i)
+            for(auto i = k; i >= 0; --i)
             {
                T term = beta * pois;
                sum += term;
@@ -94,7 +95,7 @@ namespace boost
                
                last_term = term;
             }
-            for(int i = k + 1; ; ++i)
+            for(auto i = k + 1; ; ++i)
             {
                poisf *= l2 / i;
                xtermf *= (x * (a + b + i - 2)) / (a + i - 1);
@@ -131,7 +132,7 @@ namespace boost
             // k is the starting point for iteration, and is the
             // maximum of the poisson weighting term:
             //
-            int k = itrunc(l2);
+            long long k = lltrunc(l2);
             T pois;
             if(k <= 30)
             {
@@ -174,7 +175,7 @@ namespace boost
             //
             T last_term = 0;
             std::uintmax_t count = 0;
-            for(int i = k + 1; ; ++i)
+            for(auto i = k + 1; ; ++i)
             {
                poisf *= l2 / i;
                xtermf *= (x * (a + b + i - 2)) / (a + i - 1);
@@ -195,7 +196,7 @@ namespace boost
                }
                last_term = term;
             }
-            for(int i = k; i >= 0; --i)
+            for(auto i = k; i >= 0; --i)
             {
                T term = beta * pois;
                sum += term;
@@ -535,7 +536,7 @@ namespace boost
             // k is the starting point for iteration, and is the
             // maximum of the poisson weighting term:
             //
-            int k = itrunc(l2);
+            long long k = lltrunc(l2);
             // Starting Poisson weight:
             T pois = gamma_p_derivative(T(k+1), l2, pol);
             // Starting beta term:
@@ -550,7 +551,7 @@ namespace boost
             // Stable backwards recursion first:
             //
             std::uintmax_t count = k;
-            for(int i = k; i >= 0; --i)
+            for(auto i = k; i >= 0; --i)
             {
                T term = beta * pois;
                sum += term;
@@ -566,7 +567,7 @@ namespace boost
                   beta *= (a + i - 1) / (x * (a + i + b - 1));
                }
             }
-            for(int i = k + 1; ; ++i)
+            for(auto i = k + 1; ; ++i)
             {
                poisf *= l2 / i;
                betaf *= x * (a + b + i - 1) / (a + i - 1);
diff --git a/include/boost/math/distributions/non_central_t.hpp b/include/boost/math/distributions/non_central_t.hpp
index 779652f95b..9468a56c5e 100644
--- a/include/boost/math/distributions/non_central_t.hpp
+++ b/include/boost/math/distributions/non_central_t.hpp
@@ -15,6 +15,7 @@
 #include <boost/math/distributions/normal.hpp> // for normal CDF and quantile
 #include <boost/math/distributions/students_t.hpp>
 #include <boost/math/distributions/detail/generic_quantile.hpp> // quantile
+#include <boost/math/special_functions/trunc.hpp>
 
 namespace boost
 {
@@ -43,7 +44,7 @@ namespace boost
             // cancellation errors later (test case is v = 1621286869049072.3
             // delta = 0.16212868690490723, x = 0.86987415482475994).
             //
-            int k = itrunc(d2);
+            long long k = lltrunc(d2);
             T pois;
             if(k == 0) k = 1;
             // Starting Poisson weight:
@@ -69,7 +70,7 @@ namespace boost
             //
             std::uintmax_t count = 0;
             T last_term = 0;
-            for(int i = k; i >= 0; --i)
+            for(auto i = k; i >= 0; --i)
             {
                T term = beta * pois;
                sum += term;
@@ -83,7 +84,7 @@ namespace boost
                ++count;
             }
             last_term = 0;
-            for(int i = k + 1; ; ++i)
+            for(auto i = k + 1; ; ++i)
             {
                poisf *= d2 / (i + 0.5f);
                xtermf *= (x * (v / 2 + i - 1)) / (i);
@@ -121,11 +122,11 @@ namespace boost
             // (test case is v = 561908036470413.25, delta = 0.056190803647041321,
             // x = 1.6155232703966216):
             //
-            int k = itrunc(d2);
+            long long k = lltrunc(d2);
             if(k == 0) k = 1;
             // Starting Poisson weight:
             T pois;
-            if((k < static_cast<int>(max_factorial<T>::value)) && (d2 < tools::log_max_value<T>()) && (log(d2) * k < tools::log_max_value<T>()))
+            if((k < static_cast<long long>(max_factorial<T>::value)) && (d2 < tools::log_max_value<T>()) && (log(d2) * k < tools::log_max_value<T>()))
             {
                //
                // For small k we can optimise this calculation by using
@@ -171,7 +172,7 @@ namespace boost
             //
             std::uintmax_t count = 0;
             T last_term = 0;
-            for(int i = k + 1, j = k; ; ++i, --j)
+            for(auto i = k + 1, j = k; ; ++i, --j)
             {
                poisf *= d2 / (i + 0.5f);
                xtermf *= (x * (v / 2 + i - 1)) / (i);
@@ -387,7 +388,7 @@ namespace boost
             // k is the starting point for iteration, and is the
             // maximum of the poisson weighting term:
             //
-            int k = itrunc(d2);
+            long long k = lltrunc(d2);
             T pois, xterm;
             if(k == 0)
                k = 1;
@@ -409,7 +410,7 @@ namespace boost
             // direction for recursion:
             //
             std::uintmax_t count = 0;
-            for(int i = k; i >= 0; --i)
+            for(auto i = k; i >= 0; --i)
             {
                T term = xterm * pois;
                sum += term;
@@ -425,7 +426,7 @@ namespace boost
                      "Series did not converge, closest value was %1%", sum, pol);
                }
             }
-            for(int i = k + 1; ; ++i)
+            for(auto i = k + 1; ; ++i)
             {
                poisf *= d2 / (i + 0.5f);
                xtermf *= (x * (n / 2 + i)) / (i);
diff --git a/test/Jamfile.v2 b/test/Jamfile.v2
index b35f652dbc..8aa2bdcff0 100644
--- a/test/Jamfile.v2
+++ b/test/Jamfile.v2
@@ -899,6 +899,7 @@ test-suite distribution_tests :
    [ run scipy_issue_17146.cpp ../../test/build//boost_unit_test_framework  ]
    [ run scipy_issue_17388.cpp ../../test/build//boost_unit_test_framework  ]
    [ run scipy_issue_17916.cpp ../../test/build//boost_unit_test_framework  ]
+   [ run scipy_issue_17916_nct.cpp ../../test/build//boost_unit_test_framework  ]
 ;
 
 test-suite mp :
diff --git a/test/scipy_issue_17916_nct.cpp b/test/scipy_issue_17916_nct.cpp
new file mode 100644
index 0000000000..24f87cfdd1
--- /dev/null
+++ b/test/scipy_issue_17916_nct.cpp
@@ -0,0 +1,31 @@
+// Copyright Matt Borland, 2023
+// 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/17916#issuecomment-1437826800
+
+#include <boost/math/distributions/non_central_t.hpp>
+#include "math_unit_test.hpp"
+
+int main(void)
+{
+    auto dist = boost::math::non_central_t(2.0, 482023264.0);
+    double test_pdf;
+    double test_cdf;
+
+    try
+    {
+        test_pdf = boost::math::pdf(dist, 2.0);
+        test_cdf = boost::math::cdf(dist, 2.0);
+    }
+    catch (...)
+    {
+        return 1;
+    }
+
+    CHECK_ULP_CLOSE(test_pdf, 0.0, 1);
+    CHECK_ULP_CLOSE(test_cdf, 0.0, 1);
+
+    return boost::math::test::report_errors();
+}