From fb85bdb5104842afc314bff1e3266b94faaca8d1 Mon Sep 17 00:00:00 2001
From: Matt Borland <matt@mattborland.com>
Date: Fri, 3 Feb 2023 09:10:47 -0800
Subject: [PATCH] Fix for scipy issue 17916 (#939)

---
 .../distributions/non_central_chi_squared.hpp |  8 ++---
 test/Jamfile.v2                               |  1 +
 test/scipy_issue_17916.cpp                    | 31 +++++++++++++++++++
 3 files changed, 36 insertions(+), 4 deletions(-)
 create mode 100644 test/scipy_issue_17916.cpp

diff --git a/include/boost/math/distributions/non_central_chi_squared.hpp b/include/boost/math/distributions/non_central_chi_squared.hpp
index b848851ce..5a84c734f 100644
--- a/include/boost/math/distributions/non_central_chi_squared.hpp
+++ b/include/boost/math/distributions/non_central_chi_squared.hpp
@@ -13,7 +13,7 @@
 #include <boost/math/distributions/fwd.hpp>
 #include <boost/math/special_functions/gamma.hpp> // for incomplete gamma. gamma_q
 #include <boost/math/special_functions/bessel.hpp> // for cyl_bessel_i
-#include <boost/math/special_functions/round.hpp> // for iround
+#include <boost/math/special_functions/round.hpp> // for llround
 #include <boost/math/distributions/complement.hpp> // complements
 #include <boost/math/distributions/chi_squared.hpp> // central distribution
 #include <boost/math/distributions/detail/common_error_handling.hpp> // error checks
@@ -71,7 +71,7 @@ namespace boost
             // k is chosen as the peek of the Poisson weights, which
             // will occur *before* the largest term.
             //
-            int k = iround(lambda, pol);
+            long long k = llround(lambda, pol);
             // Forwards and backwards Poisson weights:
             T poisf = boost::math::gamma_p_derivative(static_cast<T>(1 + k), lambda, pol);
             T poisb = poisf * k / lambda;
@@ -215,7 +215,7 @@ namespace boost
             // function, which ocurrs *after* the largest term in the
             // sum.
             //
-            int k = iround(del, pol);
+            long long k = llround(del, pol);
             T a = n / 2 + k;
             // Central chi squared term for forward iteration:
             T gamkf = boost::math::gamma_p(a, x, pol);
@@ -294,7 +294,7 @@ namespace boost
             T n2 = n / 2;
             T l2 = lambda / 2;
             T sum = 0;
-            int k = itrunc(l2);
+            long long k = lltrunc(l2);
             T pois = gamma_p_derivative(static_cast<T>(k + 1), l2, pol) * gamma_p_derivative(static_cast<T>(n2 + k), x2);
             if(pois == 0)
                return 0;
diff --git a/test/Jamfile.v2 b/test/Jamfile.v2
index afce7272b..4afae3bf7 100644
--- a/test/Jamfile.v2
+++ b/test/Jamfile.v2
@@ -896,6 +896,7 @@ test-suite distribution_tests :
    [ 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  ]
+   [ run scipy_issue_17916.cpp ../../test/build//boost_unit_test_framework  ]
 ;
 
 test-suite mp :
diff --git a/test/scipy_issue_17916.cpp b/test/scipy_issue_17916.cpp
new file mode 100644
index 000000000..d3b769fb3
--- /dev/null
+++ b/test/scipy_issue_17916.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
+
+#include <boost/math/distributions/non_central_chi_squared.hpp>
+#include "math_unit_test.hpp"
+
+int main(void)
+{
+    auto dist = boost::math::non_central_chi_squared(2.0, 4820232647677555.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();
+}