From d4afedc638cb34690365cbbb8853545c4ad1e77f Mon Sep 17 00:00:00 2001
From: Jacob Hass <jacobhass8@gmail.com>
Date: Sun, 28 Dec 2025 14:06:20 -0800
Subject: [PATCH 01/24] boost::math::detail log incomplete gamma function
 implemented

---
 .../boost/math/special_functions/gamma.hpp    | 45 +++++++++++++++++++
 1 file changed, 45 insertions(+)

diff --git a/include/boost/math/special_functions/gamma.hpp b/include/boost/math/special_functions/gamma.hpp
index 47b4f3d68..037c7cba0 100644
--- a/include/boost/math/special_functions/gamma.hpp
+++ b/include/boost/math/special_functions/gamma.hpp
@@ -1772,6 +1772,51 @@ BOOST_MATH_GPU_ENABLED T gamma_incomplete_imp(T a, T x, bool normalised, bool in
    return gamma_incomplete_imp_final(T(a), T(x), normalised, invert, pol, p_derivative);
 }
 
+// Calculate log of incomplete gamma function
+template <class T, class Policy>
+T lgamma_incomplete_imp_final(T a, T x, const Policy& pol)
+{
+   using namespace boost::math;  // temporary until we're in the right namespace
+
+   BOOST_MATH_STD_USING_CORE
+
+   if (((x > 1000) && ((a < x) || (fabs(a - 50) / x < 1))) || ((x > tools::log_max_value<T>() - 10) && (x > a)))
+   {
+      //
+      // Take the logarithmic version of the asymtotic expansion:
+      //
+      return log(detail::incomplete_tgamma_large_x(a, x, pol)) + a * log(x) - x - lgamma(a, pol) - log(x);
+   }
+   //
+   // Can't do better than taking the log of Q, but...
+   //
+   // Figure out whether we need P or Q, since if we calculate Q and it's too close to unity
+   // we will loose precision in the result, selection logic here is extracted from gamma_incomplete_imp_final:
+   //
+   bool need_p = false;
+   if ((x < 0.5) && (T(-0.4) / log(x) < a))
+      need_p = true;
+   else if ((x < 1.1) && (x >= 0.5) && (x * 0.75f < a))
+      need_p = true;
+   else if ((x < a) && (x >= 1.1))
+      need_p = true;
+
+   if (need_p)
+      return log1p(-gamma_p(a, x, pol), pol);
+   return log(gamma_q(a, x, pol));
+}
+
+template <class T, class Policy> 
+T lgamma_incomplete_imp(T a, T x, const Policy& pol){
+   constexpr auto function = "boost::math::lgamma_p<%1%>(%1%, %1%)";
+   if(a <= 0)
+      return policies::raise_domain_error<T>(function, "Argument a to the incomplete gamma function must be greater than zero (got a=%1%).", a, pol);
+   if(x < 0)
+      return policies::raise_domain_error<T>(function, "Argument x to the incomplete gamma function must be >= 0 (got x=%1%).", x, pol);
+
+   // If input is valid proceed as normal
+   return lgamma_incomplete_imp_final(T(a), T(x), pol)
+}
 //
 // Ratios of two gamma functions:
 //

From 49950ef147437e660637be0e4b28f224f4c8a4c9 Mon Sep 17 00:00:00 2001
From: Jacob Hass <jacobhass8@gmail.com>
Date: Sun, 28 Dec 2025 14:47:50 -0800
Subject: [PATCH 02/24] Implemented log incomplete gamma with/without policy
 [skip ci]

---
 .../boost/math/special_functions/gamma.hpp    | 25 ++++++++++++++++++-
 1 file changed, 24 insertions(+), 1 deletion(-)

diff --git a/include/boost/math/special_functions/gamma.hpp b/include/boost/math/special_functions/gamma.hpp
index 037c7cba0..a8347494b 100644
--- a/include/boost/math/special_functions/gamma.hpp
+++ b/include/boost/math/special_functions/gamma.hpp
@@ -1808,7 +1808,7 @@ T lgamma_incomplete_imp_final(T a, T x, const Policy& pol)
 
 template <class T, class Policy> 
 T lgamma_incomplete_imp(T a, T x, const Policy& pol){
-   constexpr auto function = "boost::math::lgamma_p<%1%>(%1%, %1%)";
+   constexpr auto function = "boost::math::lgamma_q<%1%>(%1%, %1%)";
    if(a <= 0)
       return policies::raise_domain_error<T>(function, "Argument a to the incomplete gamma function must be greater than zero (got a=%1%).", a, pol);
    if(x < 0)
@@ -2435,6 +2435,29 @@ BOOST_MATH_GPU_ENABLED inline tools::promote_args_t<T1, T2>
 {
    return gamma_q(a, z, policies::policy<>());
 }
+
+template <class T1, class T2, class Policy>
+inline tools::promote_args_t<T1, T2> lgamma_q(T1 a, T2 z, const Policy& /* pol */)
+{
+   typedef tools::promote_args_t<T1, T2> result_type;
+   typedef typename policies::evaluation<result_type, Policy>::type value_type;
+   typedef typename policies::normalise<
+      Policy,
+      policies::promote_float<false>,
+      policies::promote_double<false>,
+      policies::discrete_quantile<>,
+      policies::assert_undefined<> >::type forwarding_policy;
+
+   return policies::checked_narrowing_cast<result_type, forwarding_policy>(
+      detail::lgamma_incomplete_imp(static_cast<value_type>(a),
+      static_cast<value_type>(z), forwarding_policy()), "lgamma_q<%1%>(%1%, %1%)");
+}
+
+template <class T1, class T2, class Policy>
+inline tools::promote_args_t<T1, T2> lgamma_q(T1 a, T2 z)
+{
+   return lgamma_q(a, z, policies::policy<>());
+}
 //
 // Regularised lower incomplete gamma:
 //

From 4cdc487bf7999681fa0db7943e92e990c7da15af Mon Sep 17 00:00:00 2001
From: Jacob Hass <jacobhass8@gmail.com>
Date: Sun, 28 Dec 2025 16:57:55 -0800
Subject: [PATCH 03/24] Small bug fixes [skip ci]

---
 include/boost/math/special_functions/gamma.hpp | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/include/boost/math/special_functions/gamma.hpp b/include/boost/math/special_functions/gamma.hpp
index a8347494b..cb2511840 100644
--- a/include/boost/math/special_functions/gamma.hpp
+++ b/include/boost/math/special_functions/gamma.hpp
@@ -1815,7 +1815,7 @@ T lgamma_incomplete_imp(T a, T x, const Policy& pol){
       return policies::raise_domain_error<T>(function, "Argument x to the incomplete gamma function must be >= 0 (got x=%1%).", x, pol);
 
    // If input is valid proceed as normal
-   return lgamma_incomplete_imp_final(T(a), T(x), pol)
+   return lgamma_incomplete_imp_final(T(a), T(x), pol);
 }
 //
 // Ratios of two gamma functions:
@@ -2453,7 +2453,7 @@ inline tools::promote_args_t<T1, T2> lgamma_q(T1 a, T2 z, const Policy& /* pol *
       static_cast<value_type>(z), forwarding_policy()), "lgamma_q<%1%>(%1%, %1%)");
 }
 
-template <class T1, class T2, class Policy>
+template <class T1, class T2>
 inline tools::promote_args_t<T1, T2> lgamma_q(T1 a, T2 z)
 {
    return lgamma_q(a, z, policies::policy<>());

From 27300a6b60b9000e50d94ad198eaf70a620143bd Mon Sep 17 00:00:00 2001
From: Jacob Hass <jacobhass8@gmail.com>
Date: Wed, 31 Dec 2025 11:50:07 -0800
Subject: [PATCH 04/24] Merged imp_final and imp function [skip ci]

---
 .../boost/math/special_functions/gamma.hpp    | 20 ++++++++-----------
 1 file changed, 8 insertions(+), 12 deletions(-)

diff --git a/include/boost/math/special_functions/gamma.hpp b/include/boost/math/special_functions/gamma.hpp
index cb2511840..dede0dbf4 100644
--- a/include/boost/math/special_functions/gamma.hpp
+++ b/include/boost/math/special_functions/gamma.hpp
@@ -1774,12 +1774,19 @@ BOOST_MATH_GPU_ENABLED T gamma_incomplete_imp(T a, T x, bool normalised, bool in
 
 // Calculate log of incomplete gamma function
 template <class T, class Policy>
-T lgamma_incomplete_imp_final(T a, T x, const Policy& pol)
+T lgamma_incomplete_imp(T a, T x, const Policy& pol)
 {
    using namespace boost::math;  // temporary until we're in the right namespace
 
    BOOST_MATH_STD_USING_CORE
 
+   // Check for invalid inputs (a < 0 or x < 0)
+   constexpr auto function = "boost::math::lgamma_q<%1%>(%1%, %1%)";
+   if(a <= 0)
+      return policies::raise_domain_error<T>(function, "Argument a to the incomplete gamma function must be greater than zero (got a=%1%).", a, pol);
+   if(x < 0)
+      return policies::raise_domain_error<T>(function, "Argument x to the incomplete gamma function must be >= 0 (got x=%1%).", x, pol);
+
    if (((x > 1000) && ((a < x) || (fabs(a - 50) / x < 1))) || ((x > tools::log_max_value<T>() - 10) && (x > a)))
    {
       //
@@ -1806,17 +1813,6 @@ T lgamma_incomplete_imp_final(T a, T x, const Policy& pol)
    return log(gamma_q(a, x, pol));
 }
 
-template <class T, class Policy> 
-T lgamma_incomplete_imp(T a, T x, const Policy& pol){
-   constexpr auto function = "boost::math::lgamma_q<%1%>(%1%, %1%)";
-   if(a <= 0)
-      return policies::raise_domain_error<T>(function, "Argument a to the incomplete gamma function must be greater than zero (got a=%1%).", a, pol);
-   if(x < 0)
-      return policies::raise_domain_error<T>(function, "Argument x to the incomplete gamma function must be >= 0 (got x=%1%).", x, pol);
-
-   // If input is valid proceed as normal
-   return lgamma_incomplete_imp_final(T(a), T(x), pol);
-}
 //
 // Ratios of two gamma functions:
 //

From 595eb6850be82302414a1bbec0f776f825f4e08c Mon Sep 17 00:00:00 2001
From: Jacob Hass <jacobhass8@gmail.com>
Date: Wed, 31 Dec 2025 14:31:06 -0800
Subject: [PATCH 05/24] Added 3 test cases for lgamma_q [skip ci]

---
 test/test_igamma.hpp | 7 +++++++
 1 file changed, 7 insertions(+)

diff --git a/test/test_igamma.hpp b/test/test_igamma.hpp
index 3459f71d9..1d3be9fd4 100644
--- a/test/test_igamma.hpp
+++ b/test/test_igamma.hpp
@@ -255,6 +255,13 @@ void test_spots(T)
    //
    BOOST_CHECK_EQUAL(::boost::math::gamma_q(static_cast<T>(1770), static_cast<T>(1e-12)), 1);
    BOOST_CHECK_EQUAL(::boost::math::gamma_p(static_cast<T>(1770), static_cast<T>(1e-12)), 0);
+
+   //
+   // Check that lgamma_q returns correct values
+   //
+   BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast<T>(5), static_cast<T>(100)), static_cast<T>(log(1.6139305336977304790405739225035685228527400976549e-37L)), tolerance);
+   BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast<T>(22.5), static_cast<T>(2000)), static_cast<T>(-1883.4897732037716195918619632721L), tolerance * 10); // calculated via mpmath
+   BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast<T>(501.2), static_cast<T>(2000)), static_cast<T>(-810.31461624182202285737730562687L), tolerance * 10); // calculated via mpmath
    //
    // Coverage:
    //

From 8ba98b79f67c233931a4e4b28c7ff0d34ab4b1da Mon Sep 17 00:00:00 2001
From: jzmaddock <john@johnmaddock.co.uk>
Date: Thu, 1 Jan 2026 10:48:31 +0000
Subject: [PATCH 06/24] Use test values from wolframapha

---
 test/test_igamma.hpp | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/test/test_igamma.hpp b/test/test_igamma.hpp
index 1d3be9fd4..7de179c13 100644
--- a/test/test_igamma.hpp
+++ b/test/test_igamma.hpp
@@ -257,11 +257,11 @@ void test_spots(T)
    BOOST_CHECK_EQUAL(::boost::math::gamma_p(static_cast<T>(1770), static_cast<T>(1e-12)), 0);
 
    //
-   // Check that lgamma_q returns correct values
+   // Check that lgamma_q returns correct values with spot values calculated via wolframalpha  log(Q[a, x])
    //
-   BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast<T>(5), static_cast<T>(100)), static_cast<T>(log(1.6139305336977304790405739225035685228527400976549e-37L)), tolerance);
-   BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast<T>(22.5), static_cast<T>(2000)), static_cast<T>(-1883.4897732037716195918619632721L), tolerance * 10); // calculated via mpmath
-   BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast<T>(501.2), static_cast<T>(2000)), static_cast<T>(-810.31461624182202285737730562687L), tolerance * 10); // calculated via mpmath
+   BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast<T>(5), static_cast<T>(100)), static_cast<T>(-84.71697591169848944613823640968965801339401810393519310714864307L), tolerance);
+   BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast<T>(22.5), static_cast<T>(2000)), static_cast<T>(-1883.489773203771543025750308264545743305089849873060383828767138L), tolerance * 10);
+   BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast<T>(501.2), static_cast<T>(2000)), static_cast<T>(-810.2453406781655559126505101822969531699112391075198076300675402L), tolerance * 10);
    //
    // Coverage:
    //

From 82e6d3e3c7e5e02e80471a5a3aee2d420b1f2acc Mon Sep 17 00:00:00 2001
From: jzmaddock <john@johnmaddock.co.uk>
Date: Thu, 1 Jan 2026 11:04:47 +0000
Subject: [PATCH 07/24] Correct test input

---
 test/test_igamma.hpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/test/test_igamma.hpp b/test/test_igamma.hpp
index 7de179c13..546170ee4 100644
--- a/test/test_igamma.hpp
+++ b/test/test_igamma.hpp
@@ -261,7 +261,7 @@ void test_spots(T)
    //
    BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast<T>(5), static_cast<T>(100)), static_cast<T>(-84.71697591169848944613823640968965801339401810393519310714864307L), tolerance);
    BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast<T>(22.5), static_cast<T>(2000)), static_cast<T>(-1883.489773203771543025750308264545743305089849873060383828767138L), tolerance * 10);
-   BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast<T>(501.2), static_cast<T>(2000)), static_cast<T>(-810.2453406781655559126505101822969531699112391075198076300675402L), tolerance * 10);
+   BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast<T>(501.25), static_cast<T>(2000)), static_cast<T>(-810.2453406781655559126505101822969531699112391075198076300675402L), tolerance * 10);
    //
    // Coverage:
    //

From 40a0c541909e617727a6b9768b4df4a8a5a6fdc5 Mon Sep 17 00:00:00 2001
From: jzmaddock <john@johnmaddock.co.uk>
Date: Thu, 1 Jan 2026 14:18:26 +0000
Subject: [PATCH 08/24] More tests for code coverage.

---
 test/test_igamma.hpp | 17 +++++++++++++++--
 1 file changed, 15 insertions(+), 2 deletions(-)

diff --git a/test/test_igamma.hpp b/test/test_igamma.hpp
index 546170ee4..caee33e6a 100644
--- a/test/test_igamma.hpp
+++ b/test/test_igamma.hpp
@@ -260,8 +260,11 @@ void test_spots(T)
    // Check that lgamma_q returns correct values with spot values calculated via wolframalpha  log(Q[a, x])
    //
    BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast<T>(5), static_cast<T>(100)), static_cast<T>(-84.71697591169848944613823640968965801339401810393519310714864307L), tolerance);
-   BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast<T>(22.5), static_cast<T>(2000)), static_cast<T>(-1883.489773203771543025750308264545743305089849873060383828767138L), tolerance * 10);
-   BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast<T>(501.25), static_cast<T>(2000)), static_cast<T>(-810.2453406781655559126505101822969531699112391075198076300675402L), tolerance * 10);
+   BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast<T>(22.5), static_cast<T>(2000)), static_cast<T>(-1883.489773203771543025750308264545743305089849873060383828767138L), tolerance);
+   BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast<T>(501.25), static_cast<T>(2000)), static_cast<T>(-810.2453406781655559126505101822969531699112391075198076300675402L), tolerance);
+   BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast<T>(20), static_cast<T>(0.25)), static_cast<T>(-2.946458104491857816330873290969917497748067639461638294404e-31L), tolerance);
+   BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast<T>(40), static_cast<T>(0.75)), static_cast<T>(-5.930604927955460343652485525435087275997461623988991819824e-54L), tolerance);
+   BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast<T>(50), static_cast<T>(2)), static_cast<T>(-5.214301903317168085381693412994550732094621576607843973832e-51L), tolerance);
    //
    // Coverage:
    //
@@ -272,6 +275,11 @@ void test_spots(T)
    BOOST_CHECK_THROW(boost::math::gamma_q(static_cast<T>(1), static_cast<T>(-2)), std::domain_error);
    BOOST_CHECK_THROW(boost::math::gamma_p(static_cast<T>(0), static_cast<T>(2)), std::domain_error);
    BOOST_CHECK_THROW(boost::math::gamma_q(static_cast<T>(0), static_cast<T>(2)), std::domain_error);
+
+   BOOST_CHECK_THROW(boost::math::lgamma_q(static_cast<T>(-1), static_cast<T>(2)), std::domain_error);
+   BOOST_CHECK_THROW(boost::math::lgamma_q(static_cast<T>(1), static_cast<T>(-2)), std::domain_error);
+   BOOST_CHECK_THROW(boost::math::lgamma_q(static_cast<T>(0), static_cast<T>(2)), std::domain_error);
+
    BOOST_CHECK_THROW(boost::math::gamma_p_derivative(static_cast<T>(-1), static_cast<T>(2)), std::domain_error);
    BOOST_CHECK_THROW(boost::math::gamma_p_derivative(static_cast<T>(1), static_cast<T>(-2)), std::domain_error);
    BOOST_CHECK_THROW(boost::math::gamma_p_derivative(static_cast<T>(0), static_cast<T>(2)), std::domain_error);
@@ -282,6 +290,11 @@ void test_spots(T)
    BOOST_CHECK((boost::math::isnan)(boost::math::gamma_q(static_cast<T>(1), static_cast<T>(-2))));
    BOOST_CHECK((boost::math::isnan)(boost::math::gamma_p(static_cast<T>(0), static_cast<T>(2))));
    BOOST_CHECK((boost::math::isnan)(boost::math::gamma_q(static_cast<T>(0), static_cast<T>(2))));
+
+   BOOST_CHECK((boost::math::isnan)(boost::math::lgamma_q(static_cast<T>(-1), static_cast<T>(2))));
+   BOOST_CHECK((boost::math::isnan)(boost::math::lgamma_q(static_cast<T>(1), static_cast<T>(-2))));
+   BOOST_CHECK((boost::math::isnan)(boost::math::lgamma_q(static_cast<T>(0), static_cast<T>(2))));
+
    BOOST_CHECK((boost::math::isnan)(boost::math::gamma_p_derivative(static_cast<T>(-1), static_cast<T>(2))));
    BOOST_CHECK((boost::math::isnan)(boost::math::gamma_p_derivative(static_cast<T>(1), static_cast<T>(-2))));
    BOOST_CHECK((boost::math::isnan)(boost::math::gamma_p_derivative(static_cast<T>(0), static_cast<T>(2))));

From 562f4cc793e6a00ec3f2e3a2165eabde0f1674b9 Mon Sep 17 00:00:00 2001
From: Jacob Hass <jacobhass8@gmail.com>
Date: Thu, 1 Jan 2026 10:44:48 -0800
Subject: [PATCH 09/24] Added documentation for

---
 doc/sf/igamma.qbk | 25 +++++++++++++++++++++++++
 1 file changed, 25 insertions(+)

diff --git a/doc/sf/igamma.qbk b/doc/sf/igamma.qbk
index 4675928e6..3895d64c3 100644
--- a/doc/sf/igamma.qbk
+++ b/doc/sf/igamma.qbk
@@ -20,6 +20,12 @@
    template <class T1, class T2, class ``__Policy``>
    BOOST_MATH_GPU_ENABLED ``__sf_result`` gamma_q(T1 a, T2 z, const ``__Policy``&);
    
+   template <class T1, class T2>
+   BOOST_MATH_GPU_ENABLED ``__sf_result`` lgamma_q(T1 a, T2 z);
+   
+   template <class T1, class T2, class ``__Policy``>
+   BOOST_MATH_GPU_ENABLED ``__sf_result`` lgamma_q(T1 a, T2 z, const ``__Policy``&);
+
    template <class T1, class T2>
    BOOST_MATH_GPU_ENABLED ``__sf_result`` tgamma_lower(T1 a, T2 z);
    
@@ -80,6 +86,15 @@ This function changes rapidly from 1 to 0 around the point z == a:
 
 [graph gamma_q]
 
+   template <class T1, class T2>
+   BOOST_MATH_GPU_ENABLED ``__sf_result`` lgamma_q(T1 a, T2 z);
+   
+   template <class T1, class T2, class ``__Policy``>
+   BOOST_MATH_GPU_ENABLED ``__sf_result`` lgamma_q(T1 a, T2 z, const ``__Policy``&);
+
+Returns the natural log of the normalized upper incomplete gamma function 
+of a and z.  
+
    template <class T1, class T2>
    BOOST_MATH_GPU_ENABLED ``__sf_result`` tgamma_lower(T1 a, T2 z);
 
@@ -263,6 +278,16 @@ large a and x the errors will still get you eventually, although this does
 delay the inevitable much longer than other methods.  Use of /log(1+x)-x/ here
 is inspired by Temme (see references below).
 
+The natural log of the normalized upper incomplete gamma function is computed 
+as expected except when the normalized upper incomplete gamma function 
+begins to underflow. This approximately occurs at 
+
+   ((x > 1000) && ((a < x) || (fabs(a - 50) / x < 1))) || ((x > log_max_value<T>() - 10) && (x > a))
+
+in which case an expansion, for large x, of the (non-normalised) upper 
+incomplete gamma function is used. The return is then normalised by subtracting
+the log of the gamma function and adding /a log(x)-x-log(x)/.  
+
 [h4 References]
 
 * N. M. Temme, A Set of Algorithms for the Incomplete Gamma Functions,

From 1d28f1f8a77f1aba3dd8cdffd547ad6d9a959388 Mon Sep 17 00:00:00 2001
From: Jacob Hass <jacobhass8@gmail.com>
Date: Thu, 1 Jan 2026 11:24:38 -0800
Subject: [PATCH 10/24] Added test to be right below when approximation is
 taken

---
 test/test_igamma.hpp | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/test/test_igamma.hpp b/test/test_igamma.hpp
index caee33e6a..93ffce2cb 100644
--- a/test/test_igamma.hpp
+++ b/test/test_igamma.hpp
@@ -18,7 +18,7 @@
 #include <boost/type_traits/is_floating_point.hpp>
 #include <boost/array.hpp>
 #include "functor.hpp"
-
+#include <iostream>
 #include "handle_test_result.hpp"
 #include "table_type.hpp"
 
@@ -265,6 +265,7 @@ void test_spots(T)
    BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast<T>(20), static_cast<T>(0.25)), static_cast<T>(-2.946458104491857816330873290969917497748067639461638294404e-31L), tolerance);
    BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast<T>(40), static_cast<T>(0.75)), static_cast<T>(-5.930604927955460343652485525435087275997461623988991819824e-54L), tolerance);
    BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast<T>(50), static_cast<T>(2)), static_cast<T>(-5.214301903317168085381693412994550732094621576607843973832e-51L), tolerance);
+   BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast<T>(20), static_cast<T>(999)), static_cast<T>(-907.0923609280941476865481041233627829242287760441299044050137514L), tolerance);
    //
    // Coverage:
    //

From f400fcff5572f06246bed003dea686a072b83547 Mon Sep 17 00:00:00 2001
From: jzmaddock <john@johnmaddock.co.uk>
Date: Sat, 3 Jan 2026 14:12:31 +0000
Subject: [PATCH 11/24] Simplify lgamma_q internal selection logic. Add more
 tests where the logic changes method for large a,x.

---
 include/boost/math/special_functions/gamma.hpp |  2 +-
 test/test_igamma.hpp                           | 17 ++++++++++++++++-
 2 files changed, 17 insertions(+), 2 deletions(-)

diff --git a/include/boost/math/special_functions/gamma.hpp b/include/boost/math/special_functions/gamma.hpp
index dede0dbf4..3e2642abe 100644
--- a/include/boost/math/special_functions/gamma.hpp
+++ b/include/boost/math/special_functions/gamma.hpp
@@ -1787,7 +1787,7 @@ T lgamma_incomplete_imp(T a, T x, const Policy& pol)
    if(x < 0)
       return policies::raise_domain_error<T>(function, "Argument x to the incomplete gamma function must be >= 0 (got x=%1%).", x, pol);
 
-   if (((x > 1000) && ((a < x) || (fabs(a - 50) / x < 1))) || ((x > tools::log_max_value<T>() - 10) && (x > a)))
+   if (((x > 1000) || (x > tools::log_max_value<T>() - 10)) && (a + 50 < x))
    {
       //
       // Take the logarithmic version of the asymtotic expansion:
diff --git a/test/test_igamma.hpp b/test/test_igamma.hpp
index 93ffce2cb..c72cd962c 100644
--- a/test/test_igamma.hpp
+++ b/test/test_igamma.hpp
@@ -265,7 +265,22 @@ void test_spots(T)
    BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast<T>(20), static_cast<T>(0.25)), static_cast<T>(-2.946458104491857816330873290969917497748067639461638294404e-31L), tolerance);
    BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast<T>(40), static_cast<T>(0.75)), static_cast<T>(-5.930604927955460343652485525435087275997461623988991819824e-54L), tolerance);
    BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast<T>(50), static_cast<T>(2)), static_cast<T>(-5.214301903317168085381693412994550732094621576607843973832e-51L), tolerance);
-   BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast<T>(20), static_cast<T>(999)), static_cast<T>(-907.0923609280941476865481041233627829242287760441299044050137514L), tolerance);
+   // Pairs of tests that bisect the crossover condition in our code at double and then quad precision:
+   BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast<T>(10), static_cast<T>(698.75)), static_cast<T>(-652.5952453102824132865663191324423994628428404928732148525545721L), tolerance);
+   BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast<T>(10), static_cast<T>(699.25)), static_cast<T>(-653.0888168445921483147208556398158677077537551419551652934287016L), tolerance);
+   BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast<T>(10), static_cast<T>(999.75)), static_cast<T>(-950.3752463850600415679327136010171192193400042422096029239012176L), tolerance);
+   BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast<T>(10), static_cast<T>(1000.25)), static_cast<T>(-950.8707509166152482936275883547176592196662090187561681198668099L), tolerance);
+   BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast<T>(50), static_cast<T>(698.75)), static_cast<T>(-522.3277960730837166223131189587863209730608668858212533099139269L), tolerance);
+   BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast<T>(50), static_cast<T>(699.25)), static_cast<T>(-522.7927997457481265511084805522296021540768033975669071565674196L), tolerance);
+   BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast<T>(50), static_cast<T>(999.75)), static_cast<T>(-805.7977867938474339107474131612354353193501692041340771552419902L), tolerance);
+   BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast<T>(50), static_cast<T>(1000.25)), static_cast<T>(-806.2733124989172792095030711884568388681331032891850662521501582L), tolerance);
+   BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast<T>(800), static_cast<T>(999.75)), static_cast<T>(-24.33274293617739453303937714319703386675839030466670622049929011L), tolerance);
+   BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast<T>(800), static_cast<T>(1000.25)), static_cast<T>(-24.43514173634027477093666725985191846106997808357863808910970142L), tolerance);
+   // Once we get large a,x then error start to accumulate no matter what we do:
+   BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast<T>(1200), static_cast<T>(1249.75)), static_cast<T>(-2.565496161584661216769813239648606145255794643472303927896044375L), tolerance * (std::is_floating_point<T>::value ? 1 : 2));
+   BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast<T>(1200), static_cast<T>(1250.25)), static_cast<T>(-2.591934862117586205519309712218581885256650074210410262843591453L), tolerance * (std::is_same<T, float>::value ? 1 : 50));
+   BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast<T>(2200), static_cast<T>(2249.75)), static_cast<T>(-1.933779894897391651410597618307863427927461116308937004149240320L), tolerance * (std::is_floating_point<T>::value ? 1 : 10));
+   BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast<T>(2200), static_cast<T>(2250.25)), static_cast<T>(-1.950346484067948344620463026377077515919992808640737320057812268L), tolerance * (std::is_same<T, float>::value ? 1 : (std::is_floating_point<T>::value ? 50 : 100)));
    //
    // Coverage:
    //

From 6be6ac15f3f1e8aa70a0f277d954285aed180c4f Mon Sep 17 00:00:00 2001
From: jzmaddock <john@johnmaddock.co.uk>
Date: Sat, 3 Jan 2026 16:44:42 +0000
Subject: [PATCH 12/24] Add diagnostic info to tests. So we can find the cause
 of failures easier.

---
 test/test_igamma.cpp | 8 ++++----
 test/test_igamma.hpp | 3 ++-
 2 files changed, 6 insertions(+), 5 deletions(-)

diff --git a/test/test_igamma.cpp b/test/test_igamma.cpp
index ddc37f075..b113eb064 100644
--- a/test/test_igamma.cpp
+++ b/test/test_igamma.cpp
@@ -394,13 +394,13 @@ BOOST_AUTO_TEST_CASE( test_main )
    BOOST_MATH_CONTROL_FP;
 
 #ifndef BOOST_MATH_BUGGY_LARGE_FLOAT_CONSTANTS
-   test_spots(0.0F);
+   test_spots(0.0F, "float");
 #endif
-   test_spots(0.0);
+   test_spots(0.0, "double");
 #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
-   test_spots(0.0L);
+   test_spots(0.0L, "long double");
 #ifndef BOOST_MATH_NO_REAL_CONCEPT_TESTS
-   test_spots(boost::math::concepts::real_concept(0.1));
+   test_spots(boost::math::concepts::real_concept(0.1), "real_concept");
 #endif
 #endif
 
diff --git a/test/test_igamma.hpp b/test/test_igamma.hpp
index c72cd962c..a746d0e59 100644
--- a/test/test_igamma.hpp
+++ b/test/test_igamma.hpp
@@ -141,8 +141,9 @@ void test_gamma(T, const char* name)
 }
 
 template <class T>
-void test_spots(T)
+void test_spots(T, const char* name = nullptr)
 {
+   std::cout << "Testing spot values with type " << name << std::endl;
    //
    // basic sanity checks, tolerance is 10 epsilon expressed as a percentage:
    //

From 46f8b604cff6a33ff128f9c198db81aa8c21f240 Mon Sep 17 00:00:00 2001
From: jzmaddock <john@johnmaddock.co.uk>
Date: Sat, 3 Jan 2026 18:50:37 +0000
Subject: [PATCH 13/24] Adjust expected error rates.

---
 test/test_igamma.hpp | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/test/test_igamma.hpp b/test/test_igamma.hpp
index a746d0e59..3ae30c9a7 100644
--- a/test/test_igamma.hpp
+++ b/test/test_igamma.hpp
@@ -276,12 +276,12 @@ void test_spots(T, const char* name = nullptr)
    BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast<T>(50), static_cast<T>(999.75)), static_cast<T>(-805.7977867938474339107474131612354353193501692041340771552419902L), tolerance);
    BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast<T>(50), static_cast<T>(1000.25)), static_cast<T>(-806.2733124989172792095030711884568388681331032891850662521501582L), tolerance);
    BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast<T>(800), static_cast<T>(999.75)), static_cast<T>(-24.33274293617739453303937714319703386675839030466670622049929011L), tolerance);
-   BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast<T>(800), static_cast<T>(1000.25)), static_cast<T>(-24.43514173634027477093666725985191846106997808357863808910970142L), tolerance);
+   BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast<T>(800), static_cast<T>(1000.25)), static_cast<T>(-24.43514173634027477093666725985191846106997808357863808910970142L), tolerance * (boost::math::tools::digits<T>() > 54 ? 20 : 1));
    // Once we get large a,x then error start to accumulate no matter what we do:
    BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast<T>(1200), static_cast<T>(1249.75)), static_cast<T>(-2.565496161584661216769813239648606145255794643472303927896044375L), tolerance * (std::is_floating_point<T>::value ? 1 : 2));
    BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast<T>(1200), static_cast<T>(1250.25)), static_cast<T>(-2.591934862117586205519309712218581885256650074210410262843591453L), tolerance * (std::is_same<T, float>::value ? 1 : 50));
    BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast<T>(2200), static_cast<T>(2249.75)), static_cast<T>(-1.933779894897391651410597618307863427927461116308937004149240320L), tolerance * (std::is_floating_point<T>::value ? 1 : 10));
-   BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast<T>(2200), static_cast<T>(2250.25)), static_cast<T>(-1.950346484067948344620463026377077515919992808640737320057812268L), tolerance * (std::is_same<T, float>::value ? 1 : (std::is_floating_point<T>::value ? 50 : 100)));
+   BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast<T>(2200), static_cast<T>(2250.25)), static_cast<T>(-1.950346484067948344620463026377077515919992808640737320057812268L), tolerance * (std::is_same<T, float>::value ? 1 : (std::is_floating_point<T>::value ? 100 : 200)));
    //
    // Coverage:
    //

From f5f403a4a6a3e552a82bcad7d883c3355d9d171c Mon Sep 17 00:00:00 2001
From: jzmaddock <john@johnmaddock.co.uk>
Date: Sun, 4 Jan 2026 11:47:48 +0000
Subject: [PATCH 14/24] Tweak expected error rates.

---
 test/test_igamma.hpp | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/test/test_igamma.hpp b/test/test_igamma.hpp
index 3ae30c9a7..9f86fdb3f 100644
--- a/test/test_igamma.hpp
+++ b/test/test_igamma.hpp
@@ -275,10 +275,10 @@ void test_spots(T, const char* name = nullptr)
    BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast<T>(50), static_cast<T>(699.25)), static_cast<T>(-522.7927997457481265511084805522296021540768033975669071565674196L), tolerance);
    BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast<T>(50), static_cast<T>(999.75)), static_cast<T>(-805.7977867938474339107474131612354353193501692041340771552419902L), tolerance);
    BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast<T>(50), static_cast<T>(1000.25)), static_cast<T>(-806.2733124989172792095030711884568388681331032891850662521501582L), tolerance);
-   BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast<T>(800), static_cast<T>(999.75)), static_cast<T>(-24.33274293617739453303937714319703386675839030466670622049929011L), tolerance);
+   BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast<T>(800), static_cast<T>(999.75)), static_cast<T>(-24.33274293617739453303937714319703386675839030466670622049929011L), tolerance * 2);
    BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast<T>(800), static_cast<T>(1000.25)), static_cast<T>(-24.43514173634027477093666725985191846106997808357863808910970142L), tolerance * (boost::math::tools::digits<T>() > 54 ? 20 : 1));
    // Once we get large a,x then error start to accumulate no matter what we do:
-   BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast<T>(1200), static_cast<T>(1249.75)), static_cast<T>(-2.565496161584661216769813239648606145255794643472303927896044375L), tolerance * (std::is_floating_point<T>::value ? 1 : 2));
+   BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast<T>(1200), static_cast<T>(1249.75)), static_cast<T>(-2.565496161584661216769813239648606145255794643472303927896044375L), tolerance * (std::is_floating_point<T>::value ? 1 : 4));
    BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast<T>(1200), static_cast<T>(1250.25)), static_cast<T>(-2.591934862117586205519309712218581885256650074210410262843591453L), tolerance * (std::is_same<T, float>::value ? 1 : 50));
    BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast<T>(2200), static_cast<T>(2249.75)), static_cast<T>(-1.933779894897391651410597618307863427927461116308937004149240320L), tolerance * (std::is_floating_point<T>::value ? 1 : 10));
    BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast<T>(2200), static_cast<T>(2250.25)), static_cast<T>(-1.950346484067948344620463026377077515919992808640737320057812268L), tolerance * (std::is_same<T, float>::value ? 1 : (std::is_floating_point<T>::value ? 100 : 200)));

From 3363b21473d1e5b30efe87016e569745245598c2 Mon Sep 17 00:00:00 2001
From: Jacob Hass <jacobhass8@gmail.com>
Date: Mon, 5 Jan 2026 10:22:02 -0800
Subject: [PATCH 15/24] Added tests from scipy issue #8424

---
 test/test_igamma.hpp | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/test/test_igamma.hpp b/test/test_igamma.hpp
index 9f86fdb3f..d98f99c81 100644
--- a/test/test_igamma.hpp
+++ b/test/test_igamma.hpp
@@ -256,7 +256,6 @@ void test_spots(T, const char* name = nullptr)
    //
    BOOST_CHECK_EQUAL(::boost::math::gamma_q(static_cast<T>(1770), static_cast<T>(1e-12)), 1);
    BOOST_CHECK_EQUAL(::boost::math::gamma_p(static_cast<T>(1770), static_cast<T>(1e-12)), 0);
-
    //
    // Check that lgamma_q returns correct values with spot values calculated via wolframalpha  log(Q[a, x])
    //
@@ -266,6 +265,8 @@ void test_spots(T, const char* name = nullptr)
    BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast<T>(20), static_cast<T>(0.25)), static_cast<T>(-2.946458104491857816330873290969917497748067639461638294404e-31L), tolerance);
    BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast<T>(40), static_cast<T>(0.75)), static_cast<T>(-5.930604927955460343652485525435087275997461623988991819824e-54L), tolerance);
    BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast<T>(50), static_cast<T>(2)), static_cast<T>(-5.214301903317168085381693412994550732094621576607843973832e-51L), tolerance);
+   BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast<T>(500), static_cast<T>(10)), static_cast<T>(-3.79666711621207197039397438773960431648625558027046365463e-639L), tolerance * 3);
+   BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast<T>(5), static_cast<T>(1000)), static_cast<T>(-975.5430287171020511929200293377669175923128826278957569928895945L), tolerance);
    // Pairs of tests that bisect the crossover condition in our code at double and then quad precision:
    BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast<T>(10), static_cast<T>(698.75)), static_cast<T>(-652.5952453102824132865663191324423994628428404928732148525545721L), tolerance);
    BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast<T>(10), static_cast<T>(699.25)), static_cast<T>(-653.0888168445921483147208556398158677077537551419551652934287016L), tolerance);

From d8d0c0c5cac1254882821d401d46236ed542ab6d Mon Sep 17 00:00:00 2001
From: Matt Borland <matt@mattborland.com>
Date: Mon, 5 Jan 2026 13:31:08 -0500
Subject: [PATCH 16/24] Fix typo

---
 include/boost/math/special_functions/gamma.hpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/include/boost/math/special_functions/gamma.hpp b/include/boost/math/special_functions/gamma.hpp
index 3e2642abe..51bb1441c 100644
--- a/include/boost/math/special_functions/gamma.hpp
+++ b/include/boost/math/special_functions/gamma.hpp
@@ -1798,7 +1798,7 @@ T lgamma_incomplete_imp(T a, T x, const Policy& pol)
    // Can't do better than taking the log of Q, but...
    //
    // Figure out whether we need P or Q, since if we calculate Q and it's too close to unity
-   // we will loose precision in the result, selection logic here is extracted from gamma_incomplete_imp_final:
+   // we will lose precision in the result, selection logic here is extracted from gamma_incomplete_imp_final:
    //
    bool need_p = false;
    if ((x < 0.5) && (T(-0.4) / log(x) < a))

From f979925eb5b67740a2649c67421289acac947bdc Mon Sep 17 00:00:00 2001
From: Matt Borland <matt@mattborland.com>
Date: Mon, 5 Jan 2026 13:53:26 -0500
Subject: [PATCH 17/24] Add CUDA markers and testing

---
 .../boost/math/special_functions/gamma.hpp    |   6 +-
 test/cuda_jamfile                             |   2 +
 test/test_lgamma_q_double.cu                  | 102 ++++++++++++++++++
 test/test_lgamma_q_float.cu                   | 102 ++++++++++++++++++
 4 files changed, 209 insertions(+), 3 deletions(-)
 create mode 100644 test/test_lgamma_q_double.cu
 create mode 100644 test/test_lgamma_q_float.cu

diff --git a/include/boost/math/special_functions/gamma.hpp b/include/boost/math/special_functions/gamma.hpp
index 51bb1441c..458dc40a3 100644
--- a/include/boost/math/special_functions/gamma.hpp
+++ b/include/boost/math/special_functions/gamma.hpp
@@ -1774,7 +1774,7 @@ BOOST_MATH_GPU_ENABLED T gamma_incomplete_imp(T a, T x, bool normalised, bool in
 
 // Calculate log of incomplete gamma function
 template <class T, class Policy>
-T lgamma_incomplete_imp(T a, T x, const Policy& pol)
+BOOST_MATH_GPU_ENABLED T lgamma_incomplete_imp(T a, T x, const Policy& pol)
 {
    using namespace boost::math;  // temporary until we're in the right namespace
 
@@ -2433,7 +2433,7 @@ BOOST_MATH_GPU_ENABLED inline tools::promote_args_t<T1, T2>
 }
 
 template <class T1, class T2, class Policy>
-inline tools::promote_args_t<T1, T2> lgamma_q(T1 a, T2 z, const Policy& /* pol */)
+BOOST_MATH_GPU_ENABLED inline tools::promote_args_t<T1, T2> lgamma_q(T1 a, T2 z, const Policy& /* pol */)
 {
    typedef tools::promote_args_t<T1, T2> result_type;
    typedef typename policies::evaluation<result_type, Policy>::type value_type;
@@ -2450,7 +2450,7 @@ inline tools::promote_args_t<T1, T2> lgamma_q(T1 a, T2 z, const Policy& /* pol *
 }
 
 template <class T1, class T2>
-inline tools::promote_args_t<T1, T2> lgamma_q(T1 a, T2 z)
+BOOST_MATH_GPU_ENABLED inline tools::promote_args_t<T1, T2> lgamma_q(T1 a, T2 z)
 {
    return lgamma_q(a, z, policies::policy<>());
 }
diff --git a/test/cuda_jamfile b/test/cuda_jamfile
index 02dcea838..4082057fa 100644
--- a/test/cuda_jamfile
+++ b/test/cuda_jamfile
@@ -369,6 +369,8 @@ run test_gamma_p_derivative_double.cu ;
 run test_gamma_p_derivative_float.cu ;
 run test_gamma_p_inv_double.cu ;
 run test_gamma_p_inv_float.cu ;
+run test_lgamma_q_double.cu ;
+run test_lgamma_q_float.cu ;
 
 run test_log1p_double.cu ;
 run test_log1p_float.cu ;
diff --git a/test/test_lgamma_q_double.cu b/test/test_lgamma_q_double.cu
new file mode 100644
index 000000000..4326587ab
--- /dev/null
+++ b/test/test_lgamma_q_double.cu
@@ -0,0 +1,102 @@
+
+//  Copyright John Maddock 2016.
+//  Copyright Matt Borland 2024.
+//  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_MATH_PROMOTE_DOUBLE_POLICY false
+
+#include <iostream>
+#include <iomanip>
+#include <vector>
+#include <boost/math/special_functions.hpp>
+#include "cuda_managed_ptr.hpp"
+#include "stopwatch.hpp"
+
+// For the CUDA runtime routines (prefixed with "cuda_")
+#include <cuda_runtime.h>
+
+typedef double float_type;
+
+/**
+ * CUDA Kernel Device code
+ *
+ */
+__global__ void cuda_test(const float_type *in, float_type *out, int numElements)
+{
+    using std::cos;
+    int i = blockDim.x * blockIdx.x + threadIdx.x;
+
+    if (i < numElements)
+    {
+        out[i] = boost::math::lgamma_q(in[i], in[i]);
+    }
+}
+
+/**
+ * Host main routine
+ */
+int main(void)
+{
+    // Error code to check return values for CUDA calls
+    cudaError_t err = cudaSuccess;
+
+    // Print the vector length to be used, and compute its size
+    int numElements = 50000;
+    std::cout << "[Vector operation on " << numElements << " elements]" << std::endl;
+
+    // Allocate the managed input vector A
+    cuda_managed_ptr<float_type> input_vector(numElements);
+
+    // Allocate the managed output vector C
+    cuda_managed_ptr<float_type> output_vector(numElements);
+
+    // Initialize the input vectors
+    for (int i = 0; i < numElements; ++i)
+    {
+        input_vector[i] = rand()/(float_type)RAND_MAX;
+    }
+
+    // Launch the Vector Add CUDA Kernel
+    int threadsPerBlock = 1024;
+    int blocksPerGrid =(numElements + threadsPerBlock - 1) / threadsPerBlock;
+    std::cout << "CUDA kernel launch with " << blocksPerGrid << " blocks of " << threadsPerBlock << " threads" << std::endl;
+
+    watch w;
+
+    cuda_test<<<blocksPerGrid, threadsPerBlock>>>(input_vector.get(), output_vector.get(), numElements);
+    cudaDeviceSynchronize();
+
+    std::cout << "CUDA kernal done in: " << w.elapsed() << "s" << std::endl;
+
+    err = cudaGetLastError();
+
+    if (err != cudaSuccess)
+    {
+        std::cerr << "Failed to launch vectorAdd kernel (error code " << cudaGetErrorString(err) << ")!" << std::endl;
+        return EXIT_FAILURE;
+    }
+
+    // Verify that the result vector is correct
+    std::vector<float_type> results;
+    results.reserve(numElements);
+    w.reset();
+    for(int i = 0; i < numElements; ++i)
+       results.push_back(boost::math::lgamma_q(input_vector[i], input_vector[i]));
+    double t = w.elapsed();
+    // check the results
+    for(int i = 0; i < numElements; ++i)
+    {
+        if (boost::math::epsilon_difference(output_vector[i], results[i]) > 10)
+        {
+            std::cerr << "Result verification failed at element " << i << "!" << std::endl;
+            return EXIT_FAILURE;
+        }
+    }
+
+    std::cout << "Test PASSED, normal calculation time: " << t << "s" << std::endl;
+    std::cout << "Done\n";
+
+    return 0;
+}
diff --git a/test/test_lgamma_q_float.cu b/test/test_lgamma_q_float.cu
new file mode 100644
index 000000000..a38d8e13f
--- /dev/null
+++ b/test/test_lgamma_q_float.cu
@@ -0,0 +1,102 @@
+
+//  Copyright John Maddock 2016.
+//  Copyright Matt Borland 2024.
+//  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_MATH_PROMOTE_DOUBLE_POLICY false
+
+#include <iostream>
+#include <iomanip>
+#include <vector>
+#include <boost/math/special_functions.hpp>
+#include "cuda_managed_ptr.hpp"
+#include "stopwatch.hpp"
+
+// For the CUDA runtime routines (prefixed with "cuda_")
+#include <cuda_runtime.h>
+
+typedef float float_type;
+
+/**
+ * CUDA Kernel Device code
+ *
+ */
+__global__ void cuda_test(const float_type *in, float_type *out, int numElements)
+{
+    using std::cos;
+    int i = blockDim.x * blockIdx.x + threadIdx.x;
+
+    if (i < numElements)
+    {
+        out[i] = boost::math::lgamma_q(in[i], in[i]);
+    }
+}
+
+/**
+ * Host main routine
+ */
+int main(void)
+{
+    // Error code to check return values for CUDA calls
+    cudaError_t err = cudaSuccess;
+
+    // Print the vector length to be used, and compute its size
+    int numElements = 50000;
+    std::cout << "[Vector operation on " << numElements << " elements]" << std::endl;
+
+    // Allocate the managed input vector A
+    cuda_managed_ptr<float_type> input_vector(numElements);
+
+    // Allocate the managed output vector C
+    cuda_managed_ptr<float_type> output_vector(numElements);
+
+    // Initialize the input vectors
+    for (int i = 0; i < numElements; ++i)
+    {
+        input_vector[i] = rand()/(float_type)RAND_MAX;
+    }
+
+    // Launch the Vector Add CUDA Kernel
+    int threadsPerBlock = 1024;
+    int blocksPerGrid =(numElements + threadsPerBlock - 1) / threadsPerBlock;
+    std::cout << "CUDA kernel launch with " << blocksPerGrid << " blocks of " << threadsPerBlock << " threads" << std::endl;
+
+    watch w;
+
+    cuda_test<<<blocksPerGrid, threadsPerBlock>>>(input_vector.get(), output_vector.get(), numElements);
+    cudaDeviceSynchronize();
+
+    std::cout << "CUDA kernal done in: " << w.elapsed() << "s" << std::endl;
+
+    err = cudaGetLastError();
+
+    if (err != cudaSuccess)
+    {
+        std::cerr << "Failed to launch vectorAdd kernel (error code " << cudaGetErrorString(err) << ")!" << std::endl;
+        return EXIT_FAILURE;
+    }
+
+    // Verify that the result vector is correct
+    std::vector<float_type> results;
+    results.reserve(numElements);
+    w.reset();
+    for(int i = 0; i < numElements; ++i)
+       results.push_back(boost::math::lgamma_q(input_vector[i], input_vector[i]));
+    double t = w.elapsed();
+    // check the results
+    for(int i = 0; i < numElements; ++i)
+    {
+        if (boost::math::epsilon_difference(output_vector[i], results[i]) > 10)
+        {
+            std::cerr << "Result verification failed at element " << i << "!" << std::endl;
+            return EXIT_FAILURE;
+        }
+    }
+
+    std::cout << "Test PASSED, normal calculation time: " << t << "s" << std::endl;
+    std::cout << "Done\n";
+
+    return 0;
+}

From 582d769a7add27a53d4c024748b602f4b3dec9c1 Mon Sep 17 00:00:00 2001
From: Jacob Hass <jacobhass8@gmail.com>
Date: Tue, 6 Jan 2026 14:51:41 -0800
Subject: [PATCH 18/24] Relaxed precision for 128bit long doubles

---
 test/test_igamma.hpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/test/test_igamma.hpp b/test/test_igamma.hpp
index d98f99c81..febe16f21 100644
--- a/test/test_igamma.hpp
+++ b/test/test_igamma.hpp
@@ -280,7 +280,7 @@ void test_spots(T, const char* name = nullptr)
    BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast<T>(800), static_cast<T>(1000.25)), static_cast<T>(-24.43514173634027477093666725985191846106997808357863808910970142L), tolerance * (boost::math::tools::digits<T>() > 54 ? 20 : 1));
    // Once we get large a,x then error start to accumulate no matter what we do:
    BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast<T>(1200), static_cast<T>(1249.75)), static_cast<T>(-2.565496161584661216769813239648606145255794643472303927896044375L), tolerance * (std::is_floating_point<T>::value ? 1 : 4));
-   BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast<T>(1200), static_cast<T>(1250.25)), static_cast<T>(-2.591934862117586205519309712218581885256650074210410262843591453L), tolerance * (std::is_same<T, float>::value ? 1 : 50));
+   BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast<T>(1200), static_cast<T>(1250.25)), static_cast<T>(-2.591934862117586205519309712218581885256650074210410262843591453L), tolerance * ((std::numeric_limits<T>::max_digits10 >= 36) ? 50 : (std::is_same<T, float>::value ? 1 : 50)));
    BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast<T>(2200), static_cast<T>(2249.75)), static_cast<T>(-1.933779894897391651410597618307863427927461116308937004149240320L), tolerance * (std::is_floating_point<T>::value ? 1 : 10));
    BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast<T>(2200), static_cast<T>(2250.25)), static_cast<T>(-1.950346484067948344620463026377077515919992808640737320057812268L), tolerance * (std::is_same<T, float>::value ? 1 : (std::is_floating_point<T>::value ? 100 : 200)));
    //

From da1a7e40315fd56251d9dff47e47869fca7c4ffa Mon Sep 17 00:00:00 2001
From: Jacob Hass <72838561+JacobHass8@users.noreply.github.com>
Date: Tue, 6 Jan 2026 19:47:56 -0800
Subject: [PATCH 19/24] Increased tolerance more for lgamma_q test case

---
 test/test_igamma.hpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/test/test_igamma.hpp b/test/test_igamma.hpp
index febe16f21..c2b6820bc 100644
--- a/test/test_igamma.hpp
+++ b/test/test_igamma.hpp
@@ -280,7 +280,7 @@ void test_spots(T, const char* name = nullptr)
    BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast<T>(800), static_cast<T>(1000.25)), static_cast<T>(-24.43514173634027477093666725985191846106997808357863808910970142L), tolerance * (boost::math::tools::digits<T>() > 54 ? 20 : 1));
    // Once we get large a,x then error start to accumulate no matter what we do:
    BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast<T>(1200), static_cast<T>(1249.75)), static_cast<T>(-2.565496161584661216769813239648606145255794643472303927896044375L), tolerance * (std::is_floating_point<T>::value ? 1 : 4));
-   BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast<T>(1200), static_cast<T>(1250.25)), static_cast<T>(-2.591934862117586205519309712218581885256650074210410262843591453L), tolerance * ((std::numeric_limits<T>::max_digits10 >= 36) ? 50 : (std::is_same<T, float>::value ? 1 : 50)));
+   BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast<T>(1200), static_cast<T>(1250.25)), static_cast<T>(-2.591934862117586205519309712218581885256650074210410262843591453L), tolerance * ((std::numeric_limits<T>::max_digits10 >= 36) ? 500 : (std::is_same<T, float>::value ? 1 : 50)));
    BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast<T>(2200), static_cast<T>(2249.75)), static_cast<T>(-1.933779894897391651410597618307863427927461116308937004149240320L), tolerance * (std::is_floating_point<T>::value ? 1 : 10));
    BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast<T>(2200), static_cast<T>(2250.25)), static_cast<T>(-1.950346484067948344620463026377077515919992808640737320057812268L), tolerance * (std::is_same<T, float>::value ? 1 : (std::is_floating_point<T>::value ? 100 : 200)));
    //

From a3bae106849a9d03cd5741b0724ebf1ddae44d40 Mon Sep 17 00:00:00 2001
From: jzmaddock <john@johnmaddock.co.uk>
Date: Wed, 7 Jan 2026 12:39:34 +0000
Subject: [PATCH 20/24] lgamma_q: Tweak expected error rates. Hook up concept
 checks. Hook up include test. Add forward declarations.

---
 include/boost/math/special_functions/math_fwd.hpp | 9 +++++++++
 test/compile_test/instantiate.hpp                 | 8 ++++++++
 test/compile_test/sf_gamma_incl_test.cpp          | 6 ++++++
 test/test_igamma.hpp                              | 2 +-
 4 files changed, 24 insertions(+), 1 deletion(-)

diff --git a/include/boost/math/special_functions/math_fwd.hpp b/include/boost/math/special_functions/math_fwd.hpp
index 0e0815790..ca16fc8d9 100644
--- a/include/boost/math/special_functions/math_fwd.hpp
+++ b/include/boost/math/special_functions/math_fwd.hpp
@@ -561,6 +561,12 @@ namespace boost
    template <class RT1, class RT2, class Policy>
    BOOST_MATH_GPU_ENABLED tools::promote_args_t<RT1, RT2> gamma_q(RT1 a, RT2 z, const Policy&);
 
+   template <class RT1, class RT2>
+   BOOST_MATH_GPU_ENABLED tools::promote_args_t<RT1, RT2> lgamma_q(RT1 a, RT2 z);
+
+   template <class RT1, class RT2, class Policy>
+   BOOST_MATH_GPU_ENABLED tools::promote_args_t<RT1, RT2> lgamma_q(RT1 a, RT2 z, const Policy&);
+
    template <class RT1, class RT2>
    BOOST_MATH_GPU_ENABLED tools::promote_args_t<RT1, RT2> gamma_p(RT1 a, RT2 z);
 
@@ -1516,6 +1522,9 @@ namespace boost
 \
    template <class RT1, class RT2>\
    BOOST_MATH_GPU_ENABLED inline boost::math::tools::promote_args_t<RT1, RT2> gamma_q(RT1 a, RT2 z){ return boost::math::gamma_q(a, z, Policy()); }\
+\
+   template <class RT1, class RT2>\
+   BOOST_MATH_GPU_ENABLED inline boost::math::tools::promote_args_t<RT1, RT2> lgamma_q(RT1 a, RT2 z){ return boost::math::lgamma_q(a, z, Policy()); }\
 \
    template <class RT1, class RT2>\
    BOOST_MATH_GPU_ENABLED inline boost::math::tools::promote_args_t<RT1, RT2> gamma_p(RT1 a, RT2 z){ return boost::math::gamma_p(a, z, Policy()); }\
diff --git a/test/compile_test/instantiate.hpp b/test/compile_test/instantiate.hpp
index 26fa53d7b..1d186aff7 100644
--- a/test/compile_test/instantiate.hpp
+++ b/test/compile_test/instantiate.hpp
@@ -263,6 +263,7 @@ void instantiate(RealType)
    boost::math::tgamma_lower(v1, v2);
    boost::math::gamma_p(v1, v2);
    boost::math::gamma_q(v1, v2);
+   boost::math::lgamma_q(v1, v2);
    boost::math::gamma_p_inv(v1, v2);
    boost::math::gamma_q_inv(v1, v2);
    boost::math::gamma_p_inva(v1, v2);
@@ -542,6 +543,7 @@ void instantiate(RealType)
    boost::math::tgamma_lower(v1 * 1, v2 - 0);
    boost::math::gamma_p(v1 * 1, v2 + 0);
    boost::math::gamma_q(v1 * 1, v2 + 0);
+   boost::math::lgamma_q(v1 * 1, v2 + 0);
    boost::math::gamma_p_inv(v1 * 1, v2 + 0);
    boost::math::gamma_q_inv(v1 * 1, v2 + 0);
    boost::math::gamma_p_inva(v1 * 1, v2 + 0);
@@ -793,6 +795,7 @@ void instantiate(RealType)
    boost::math::tgamma_lower(v1, v2, pol);
    boost::math::gamma_p(v1, v2, pol);
    boost::math::gamma_q(v1, v2, pol);
+   boost::math::lgamma_q(v1, v2, pol);
    boost::math::gamma_p_inv(v1, v2, pol);
    boost::math::gamma_q_inv(v1, v2, pol);
    boost::math::gamma_p_inva(v1, v2, pol);
@@ -1070,6 +1073,7 @@ void instantiate(RealType)
    test::tgamma_lower(v1, v2);
    test::gamma_p(v1, v2);
    test::gamma_q(v1, v2);
+   test::lgamma_q(v1, v2);
    test::gamma_p_inv(v1, v2);
    test::gamma_q_inv(v1, v2);
    test::gamma_p_inva(v1, v2);
@@ -1351,6 +1355,7 @@ void instantiate_mixed(RealType)
    boost::math::gamma_p(i, s);
    boost::math::gamma_p(fr, lr);
    boost::math::gamma_q(i, s);
+   boost::math::lgamma_q(i, s);
    boost::math::gamma_q(fr, lr);
    boost::math::gamma_p_inv(i, fr);
    boost::math::gamma_q_inv(s, fr);
@@ -1566,6 +1571,7 @@ void instantiate_mixed(RealType)
    boost::math::gamma_p(i, s, pol);
    boost::math::gamma_p(fr, lr, pol);
    boost::math::gamma_q(i, s, pol);
+   boost::math::lgamma_q(i, s, pol);
    boost::math::gamma_q(fr, lr, pol);
    boost::math::gamma_p_inv(i, fr, pol);
    boost::math::gamma_q_inv(s, fr, pol);
@@ -1777,7 +1783,9 @@ void instantiate_mixed(RealType)
    test::gamma_p(i, s);
    test::gamma_p(fr, lr);
    test::gamma_q(i, s);
+   test::lgamma_q(i, s);
    test::gamma_q(fr, lr);
+   test::lgamma_q(fr, lr);
    test::gamma_p_inv(i, fr);
    test::gamma_q_inv(s, fr);
    test::gamma_p_inva(i, lr);
diff --git a/test/compile_test/sf_gamma_incl_test.cpp b/test/compile_test/sf_gamma_incl_test.cpp
index f7ad18e9c..74ab85b2d 100644
--- a/test/compile_test/sf_gamma_incl_test.cpp
+++ b/test/compile_test/sf_gamma_incl_test.cpp
@@ -39,6 +39,12 @@ void compile_and_link_test()
    check_result<long double>(boost::math::gamma_q<long double>(l, l));
 #endif
 
+   check_result<float>(boost::math::lgamma_q<float>(f, f));
+   check_result<double>(boost::math::lgamma_q<double>(d, d));
+#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
+   check_result<long double>(boost::math::lgamma_q<long double>(l, l));
+#endif
+
    check_result<float>(boost::math::gamma_p_inv<float>(f, f));
    check_result<double>(boost::math::gamma_p_inv<double>(d, d));
 #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
diff --git a/test/test_igamma.hpp b/test/test_igamma.hpp
index c2b6820bc..85b86f421 100644
--- a/test/test_igamma.hpp
+++ b/test/test_igamma.hpp
@@ -280,7 +280,7 @@ void test_spots(T, const char* name = nullptr)
    BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast<T>(800), static_cast<T>(1000.25)), static_cast<T>(-24.43514173634027477093666725985191846106997808357863808910970142L), tolerance * (boost::math::tools::digits<T>() > 54 ? 20 : 1));
    // Once we get large a,x then error start to accumulate no matter what we do:
    BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast<T>(1200), static_cast<T>(1249.75)), static_cast<T>(-2.565496161584661216769813239648606145255794643472303927896044375L), tolerance * (std::is_floating_point<T>::value ? 1 : 4));
-   BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast<T>(1200), static_cast<T>(1250.25)), static_cast<T>(-2.591934862117586205519309712218581885256650074210410262843591453L), tolerance * ((std::numeric_limits<T>::max_digits10 >= 36) ? 500 : (std::is_same<T, float>::value ? 1 : 50)));
+   BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast<T>(1200), static_cast<T>(1250.25)), static_cast<T>(-2.591934862117586205519309712218581885256650074210410262843591453L), tolerance * ((std::numeric_limits<T>::max_digits10 >= 36) ? 750 : (std::is_same<T, float>::value ? 1 : 50)));
    BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast<T>(2200), static_cast<T>(2249.75)), static_cast<T>(-1.933779894897391651410597618307863427927461116308937004149240320L), tolerance * (std::is_floating_point<T>::value ? 1 : 10));
    BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast<T>(2200), static_cast<T>(2250.25)), static_cast<T>(-1.950346484067948344620463026377077515919992808640737320057812268L), tolerance * (std::is_same<T, float>::value ? 1 : (std::is_floating_point<T>::value ? 100 : 200)));
    //

From 266bf67061d2d4ef1e3dc04b0354a17257f4fa17 Mon Sep 17 00:00:00 2001
From: Jacob Hass <jacobhass8@gmail.com>
Date: Wed, 7 Jan 2026 10:25:44 -0800
Subject: [PATCH 21/24] Relaxed error rates more for type real_concept

---
 test/test_igamma.hpp | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/test/test_igamma.hpp b/test/test_igamma.hpp
index 85b86f421..b03419667 100644
--- a/test/test_igamma.hpp
+++ b/test/test_igamma.hpp
@@ -18,7 +18,6 @@
 #include <boost/type_traits/is_floating_point.hpp>
 #include <boost/array.hpp>
 #include "functor.hpp"
-#include <iostream>
 #include "handle_test_result.hpp"
 #include "table_type.hpp"
 
@@ -280,7 +279,7 @@ void test_spots(T, const char* name = nullptr)
    BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast<T>(800), static_cast<T>(1000.25)), static_cast<T>(-24.43514173634027477093666725985191846106997808357863808910970142L), tolerance * (boost::math::tools::digits<T>() > 54 ? 20 : 1));
    // Once we get large a,x then error start to accumulate no matter what we do:
    BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast<T>(1200), static_cast<T>(1249.75)), static_cast<T>(-2.565496161584661216769813239648606145255794643472303927896044375L), tolerance * (std::is_floating_point<T>::value ? 1 : 4));
-   BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast<T>(1200), static_cast<T>(1250.25)), static_cast<T>(-2.591934862117586205519309712218581885256650074210410262843591453L), tolerance * ((std::numeric_limits<T>::max_digits10 >= 36) ? 750 : (std::is_same<T, float>::value ? 1 : 50)));
+   BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast<T>(1200), static_cast<T>(1250.25)), static_cast<T>(-2.591934862117586205519309712218581885256650074210410262843591453L), tolerance * ((std::numeric_limits<T>::max_digits10 >= 36 || std::is_same<T, boost::math::concepts::real_concept>::value) ? 750 : (std::is_same<T, float>::value ? 1 : 50))); // Test fails on ARM64 and s390x long doubles and real_concept types unless tolerance is adjusted 
    BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast<T>(2200), static_cast<T>(2249.75)), static_cast<T>(-1.933779894897391651410597618307863427927461116308937004149240320L), tolerance * (std::is_floating_point<T>::value ? 1 : 10));
    BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast<T>(2200), static_cast<T>(2250.25)), static_cast<T>(-1.950346484067948344620463026377077515919992808640737320057812268L), tolerance * (std::is_same<T, float>::value ? 1 : (std::is_floating_point<T>::value ? 100 : 200)));
    //

From ce9aba42a8c4b227590f763584d5d10614c29342 Mon Sep 17 00:00:00 2001
From: jzmaddock <john@johnmaddock.co.uk>
Date: Thu, 8 Jan 2026 12:26:13 +0000
Subject: [PATCH 22/24] lgamma_q: windows gcc tolerance changes.

---
 test/test_igamma.hpp | 7 ++++++-
 1 file changed, 6 insertions(+), 1 deletion(-)

diff --git a/test/test_igamma.hpp b/test/test_igamma.hpp
index b03419667..2d13448aa 100644
--- a/test/test_igamma.hpp
+++ b/test/test_igamma.hpp
@@ -263,7 +263,12 @@ void test_spots(T, const char* name = nullptr)
    BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast<T>(501.25), static_cast<T>(2000)), static_cast<T>(-810.2453406781655559126505101822969531699112391075198076300675402L), tolerance);
    BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast<T>(20), static_cast<T>(0.25)), static_cast<T>(-2.946458104491857816330873290969917497748067639461638294404e-31L), tolerance);
    BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast<T>(40), static_cast<T>(0.75)), static_cast<T>(-5.930604927955460343652485525435087275997461623988991819824e-54L), tolerance);
-   BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast<T>(50), static_cast<T>(2)), static_cast<T>(-5.214301903317168085381693412994550732094621576607843973832e-51L), tolerance);
+   BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast<T>(50), static_cast<T>(2)), static_cast<T>(-5.214301903317168085381693412994550732094621576607843973832e-51L), 
+#if defined(__CYGWIN__) || defined(__MINGW32__)
+      tolerance * 2);
+#else
+      tolerance);
+#endif
    BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast<T>(500), static_cast<T>(10)), static_cast<T>(-3.79666711621207197039397438773960431648625558027046365463e-639L), tolerance * 3);
    BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast<T>(5), static_cast<T>(1000)), static_cast<T>(-975.5430287171020511929200293377669175923128826278957569928895945L), tolerance);
    // Pairs of tests that bisect the crossover condition in our code at double and then quad precision:

From c366b2d9050356eab5db2bcfbaca55b4cf41a848 Mon Sep 17 00:00:00 2001
From: jzmaddock <john@johnmaddock.co.uk>
Date: Thu, 8 Jan 2026 15:56:12 +0000
Subject: [PATCH 23/24] lgamma_q: correct PP logic.

---
 test/test_igamma.hpp | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/test/test_igamma.hpp b/test/test_igamma.hpp
index 2d13448aa..d208072d0 100644
--- a/test/test_igamma.hpp
+++ b/test/test_igamma.hpp
@@ -263,12 +263,12 @@ void test_spots(T, const char* name = nullptr)
    BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast<T>(501.25), static_cast<T>(2000)), static_cast<T>(-810.2453406781655559126505101822969531699112391075198076300675402L), tolerance);
    BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast<T>(20), static_cast<T>(0.25)), static_cast<T>(-2.946458104491857816330873290969917497748067639461638294404e-31L), tolerance);
    BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast<T>(40), static_cast<T>(0.75)), static_cast<T>(-5.930604927955460343652485525435087275997461623988991819824e-54L), tolerance);
-   BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast<T>(50), static_cast<T>(2)), static_cast<T>(-5.214301903317168085381693412994550732094621576607843973832e-51L), 
 #if defined(__CYGWIN__) || defined(__MINGW32__)
-      tolerance * 2);
+   T gcc_win_mul = 2;
 #else
-      tolerance);
+   T gcc_win_mul = 1;
 #endif
+   BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast<T>(50), static_cast<T>(2)), static_cast<T>(-5.214301903317168085381693412994550732094621576607843973832e-51L), tolerance * gcc_win_mul);
    BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast<T>(500), static_cast<T>(10)), static_cast<T>(-3.79666711621207197039397438773960431648625558027046365463e-639L), tolerance * 3);
    BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast<T>(5), static_cast<T>(1000)), static_cast<T>(-975.5430287171020511929200293377669175923128826278957569928895945L), tolerance);
    // Pairs of tests that bisect the crossover condition in our code at double and then quad precision:

From 2a6fb336dcc874bb1264b7d055ed043fb889568e Mon Sep 17 00:00:00 2001
From: jzmaddock <john@johnmaddock.co.uk>
Date: Thu, 8 Jan 2026 16:43:19 +0000
Subject: [PATCH 24/24] lgamma_q: one more go at getting the tolerances
 correct.

---
 test/test_igamma.hpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/test/test_igamma.hpp b/test/test_igamma.hpp
index d208072d0..88f0bca17 100644
--- a/test/test_igamma.hpp
+++ b/test/test_igamma.hpp
@@ -269,7 +269,7 @@ void test_spots(T, const char* name = nullptr)
    T gcc_win_mul = 1;
 #endif
    BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast<T>(50), static_cast<T>(2)), static_cast<T>(-5.214301903317168085381693412994550732094621576607843973832e-51L), tolerance * gcc_win_mul);
-   BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast<T>(500), static_cast<T>(10)), static_cast<T>(-3.79666711621207197039397438773960431648625558027046365463e-639L), tolerance * 3);
+   BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast<T>(500), static_cast<T>(10)), static_cast<T>(-3.79666711621207197039397438773960431648625558027046365463e-639L), tolerance * 3 * gcc_win_mul);
    BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast<T>(5), static_cast<T>(1000)), static_cast<T>(-975.5430287171020511929200293377669175923128826278957569928895945L), tolerance);
    // Pairs of tests that bisect the crossover condition in our code at double and then quad precision:
    BOOST_CHECK_CLOSE(::boost::math::lgamma_q(static_cast<T>(10), static_cast<T>(698.75)), static_cast<T>(-652.5952453102824132865663191324423994628428404928732148525545721L), tolerance);