From 7503be4cfaf38533fff0c6206b0b1461503743a0 Mon Sep 17 00:00:00 2001
From: Kilian Kilger <kkilger@gmail.com>
Date: Tue, 1 Apr 2025 21:18:02 +0200
Subject: [PATCH 1/3] Fix lower incomplete gamma functions with x = 0

In this case, the errno error handling did not work correctly,
as internal functions where accidently setting it, although no
overflow happens.

Fixes #1249.
---
 .../boost/math/special_functions/gamma.hpp    |  14 +-
 test/Jamfile.v2                               |   1 +
 test/git_issue_1249.cpp                       | 145 ++++++++++++++++++
 3 files changed, 155 insertions(+), 5 deletions(-)
 create mode 100644 test/git_issue_1249.cpp

diff --git a/include/boost/math/special_functions/gamma.hpp b/include/boost/math/special_functions/gamma.hpp
index 809d610c18..6fb916708b 100644
--- a/include/boost/math/special_functions/gamma.hpp
+++ b/include/boost/math/special_functions/gamma.hpp
@@ -908,7 +908,7 @@ BOOST_MATH_GPU_ENABLED T full_igamma_prefix(T a, T z, const Policy& pol)
 {
    BOOST_MATH_STD_USING
 
-   if (z > tools::max_value<T>())
+   if (z > tools::max_value<T>() || (a > 0 && z == 0))
       return 0;
 
    T alz = a * log(z);
@@ -962,7 +962,7 @@ template <class T, class Policy, class Lanczos>
 BOOST_MATH_GPU_ENABLED T regularised_gamma_prefix(T a, T z, const Policy& pol, const Lanczos& l)
 {
    BOOST_MATH_STD_USING
-   if (z >= tools::max_value<T>())
+   if (z >= tools::max_value<T>() || (a > 0 && z == 0))
       return 0;
    T agh = a + static_cast<T>(Lanczos::g()) - T(0.5);
    T prefix;
@@ -1297,7 +1297,11 @@ BOOST_MATH_GPU_ENABLED T gamma_incomplete_imp_final(T a, T x, bool normalised, b
 
    int eval_method;
 
-   if(is_int && (x > 0.6))
+   if (x == 0)
+   {
+      eval_method = 2;
+   }
+   else if(is_int && (x > 0.6))
    {
       // calculate Q via finite sum:
       invert = !invert;
@@ -1627,7 +1631,7 @@ BOOST_MATH_GPU_ENABLED T gamma_incomplete_imp_final(T a, T x, bool normalised, b
       #endif
       result = gam - result;
    }
-   if(p_derivative)
+   if(p_derivative && x > 0)
    {
       //
       // Need to convert prefix term to derivative:
@@ -1660,7 +1664,7 @@ BOOST_MATH_GPU_ENABLED T gamma_incomplete_imp(T a, T x, bool normalised, bool in
 
    T result = 0; // Just to avoid warning C4701: potentially uninitialized local variable 'result' used
 
-   if(a >= max_factorial<T>::value && !normalised)
+   if(x > 0 && a >= max_factorial<T>::value && !normalised)
    {
       //
       // When we're computing the non-normalized incomplete gamma
diff --git a/test/Jamfile.v2 b/test/Jamfile.v2
index 4adb29d160..8a4ffb3725 100644
--- a/test/Jamfile.v2
+++ b/test/Jamfile.v2
@@ -194,6 +194,7 @@ test-suite special_fun :
    [ run git_issue_1139.cpp ]
    [ run git_issue_1175.cpp ]
    [ run git_issue_1194.cpp ]
+   [ run git_issue_1249.cpp /boost/test//boost_unit_test_framework  ]
    [ run special_functions_test.cpp /boost/test//boost_unit_test_framework  ]
    [ run test_airy.cpp test_instances//test_instances pch_light /boost/test//boost_unit_test_framework  ]
    [ run test_bessel_j.cpp test_instances//test_instances pch_light /boost/test//boost_unit_test_framework  ]
diff --git a/test/git_issue_1249.cpp b/test/git_issue_1249.cpp
new file mode 100644
index 0000000000..e9af0b1701
--- /dev/null
+++ b/test/git_issue_1249.cpp
@@ -0,0 +1,145 @@
+//  (C) Copyright Kilian Kilger 2025.
+//  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_TEST_MAIN
+
+#include <boost/test/unit_test.hpp>
+#include <boost/test/tools/floating_point_comparison.hpp>
+#include <boost/test/unit_test.hpp>
+#include <boost/test/results_collector.hpp>
+#include <boost/math/special_functions/gamma.hpp>
+#include <boost/multiprecision/cpp_bin_float.hpp>
+
+using namespace std;
+using namespace boost::math;
+using namespace boost::math::policies;
+using namespace boost::multiprecision;
+
+typedef policy<
+   policies::domain_error<errno_on_error>,
+   policies::pole_error<errno_on_error>,
+   policies::overflow_error<errno_on_error>,
+   policies::evaluation_error<errno_on_error>
+> c_policy;
+
+template<typename T>
+struct test_lower
+{
+   T operator()(T a, T x) const
+   {
+      return tgamma_lower(a, x, c_policy());
+   }
+
+   T expected(T a) const
+   {
+      return T(0.0);
+   }
+};
+
+template<typename T>
+struct test_upper
+{
+   T operator()(T a, T x) const
+   {
+      return tgamma(a, x, c_policy());
+   }
+   T expected(T a) const
+   {
+      return tgamma(a, c_policy());
+   }
+};
+
+template<typename T>
+struct test_gamma_p
+{
+   T operator()(T a, T x) const
+   {
+      return gamma_p(a, x, c_policy());
+   }
+   T expected(T) const
+   {
+      return T(0.0);
+   }
+};
+
+template<typename T>
+struct test_gamma_q
+{
+   T operator()(T a, T x) const
+   {
+      return gamma_q(a, x, c_policy());
+   }
+   T expected(T) const
+   {
+      return T(1.0);
+   }
+};
+
+template<typename T, template<typename> class Fun>
+void test_impl(T a)
+{
+   Fun<T> fn;
+   errno = 0;
+   T x = T(0.0);
+   T result = fn(a, x);
+   int saveErrno = errno;
+
+   errno = 0;
+
+   T expected = fn.expected(a);
+
+   BOOST_CHECK(errno == saveErrno);
+   BOOST_CHECK_EQUAL(result, expected);
+}
+
+template<template<typename> class Fun, typename T>
+void test_type_dispatch(T val)
+{
+   if (val <= (std::numeric_limits<float>::max)())
+      test_impl<float, Fun>(static_cast<float>(val));
+   if (val <= (std::numeric_limits<double>::max)())
+      test_impl<double, Fun>(static_cast<double>(val));
+#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
+   test_impl<long double, Fun>(static_cast<long double>(val));
+#endif
+   test_impl<cpp_bin_float_50, Fun>(static_cast<cpp_bin_float_50>(val));
+}
+
+template<template<typename> class Fun>
+void test_impl()
+{
+   test_type_dispatch<Fun, double>(1.0);
+   test_type_dispatch<Fun, double>(0.1);
+   test_type_dispatch<Fun, double>(0.5);
+   test_type_dispatch<Fun, double>(0.6);
+   test_type_dispatch<Fun, double>(1.3);
+   test_type_dispatch<Fun, double>(1.5);
+   test_type_dispatch<Fun, double>(2);
+   test_type_dispatch<Fun, double>(100);
+   test_type_dispatch<Fun, double>((std::numeric_limits<float>::max)());
+   test_type_dispatch<Fun, double>((std::numeric_limits<double>::max)());
+#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS   
+   test_type_dispatch<Fun, long double>((std::numeric_limits<long double>::max)());
+#endif
+}
+
+void test_derivative()
+{
+   double derivative = 0;
+   double result = boost::math::detail::gamma_incomplete_imp(1.0, 0.0, 
+     true, false, c_policy(), &derivative);
+   BOOST_CHECK(errno == 0);
+   BOOST_CHECK_EQUAL(derivative, 0);
+   BOOST_CHECK_EQUAL(result, 0);
+}
+
+BOOST_AUTO_TEST_CASE( test_main )
+{
+   test_impl<test_lower>();
+   test_impl<test_upper>();
+   test_impl<test_gamma_p>();
+   test_impl<test_gamma_q>();
+   test_derivative();
+}

From 410f74821b2eedc5884fe82a578fee55626edf59 Mon Sep 17 00:00:00 2001
From: Kilian Kilger <kkilger@gmail.com>
Date: Fri, 18 Apr 2025 12:04:17 +0200
Subject: [PATCH 2/3] Implement suggestions by jzmaddock

---
 .../boost/math/special_functions/gamma.hpp    | 24 +++++++++----------
 test/git_issue_1249.cpp                       |  6 ++---
 2 files changed, 15 insertions(+), 15 deletions(-)

diff --git a/include/boost/math/special_functions/gamma.hpp b/include/boost/math/special_functions/gamma.hpp
index 6fb916708b..8e8ebd33ad 100644
--- a/include/boost/math/special_functions/gamma.hpp
+++ b/include/boost/math/special_functions/gamma.hpp
@@ -1481,14 +1481,14 @@ BOOST_MATH_GPU_ENABLED T gamma_incomplete_imp_final(T a, T x, bool normalised, b
                #ifdef BOOST_MATH_HAS_NVRTC
                if (boost::math::is_same_v<T, float>)
                {
-                  init_value = (normalised ? 1 : ::tgammaf(a));
+                  init_value = (normalised ? T(1) : ::tgammaf(a));
                }
                else
                {
-                  init_value = (normalised ? 1 : ::tgamma(a));
+                  init_value = (normalised ? T(1) : ::tgamma(a));
                }
                #else
-               init_value = (normalised ? 1 : boost::math::tgamma(a, pol));
+               init_value = (normalised ? T(1) : boost::math::tgamma(a, pol));
                #endif
 
                if(normalised || (result >= 1) || (tools::max_value<T>() * result > init_value))
@@ -1620,29 +1620,29 @@ BOOST_MATH_GPU_ENABLED T gamma_incomplete_imp_final(T a, T x, bool normalised, b
       T gam;
       if (boost::math::is_same_v<T, float>)
       {
-         gam = normalised ? 1 : ::tgammaf(a);
+         gam = normalised ? T(1) : ::tgammaf(a);
       }
       else
       {
-         gam = normalised ? 1 : ::tgamma(a);
+         gam = normalised ? T(1) : ::tgamma(a);
       }
       #else
-      T gam = normalised ? 1 : boost::math::tgamma(a, pol);
+      T gam = normalised ? T(1) : boost::math::tgamma(a, pol);
       #endif
       result = gam - result;
    }
-   if(p_derivative && x > 0)
+   if(p_derivative)
    {
       //
       // Need to convert prefix term to derivative:
       //
-      if((x < 1) && (tools::max_value<T>() * x < *p_derivative))
+      if(x == 0 || ((x < 1) && (tools::max_value<T>() * x < *p_derivative)))
       {
          // overflow, just return an arbitrarily large value:
          *p_derivative = tools::max_value<T>() / 2;
       }
-
-      *p_derivative /= x;
+      else
+         *p_derivative /= x;
    }
 
    return result;
@@ -2110,8 +2110,8 @@ BOOST_MATH_GPU_ENABLED T gamma_p_derivative_imp(T a, T x, const Policy& pol)
    //
    if(x == 0)
    {
-      return (a > 1) ? 0 :
-         (a == 1) ? 1 : policies::raise_overflow_error<T>("boost::math::gamma_p_derivative<%1%>(%1%, %1%)", nullptr, pol);
+      return (a > 1) ? T(0) :
+         (a == 1) ? T(1) : policies::raise_overflow_error<T>("boost::math::gamma_p_derivative<%1%>(%1%, %1%)", nullptr, pol);
    }
    //
    // Normal case:
diff --git a/test/git_issue_1249.cpp b/test/git_issue_1249.cpp
index e9af0b1701..5e0fa90d2e 100644
--- a/test/git_issue_1249.cpp
+++ b/test/git_issue_1249.cpp
@@ -127,11 +127,11 @@ void test_impl()
 
 void test_derivative()
 {
+   using namespace boost::math::detail;
    double derivative = 0;
-   double result = boost::math::detail::gamma_incomplete_imp(1.0, 0.0, 
-     true, false, c_policy(), &derivative);
+   double result = gamma_incomplete_imp(1.0, 0.0, true, false, c_policy(), &derivative);
    BOOST_CHECK(errno == 0);
-   BOOST_CHECK_EQUAL(derivative, 0);
+   BOOST_CHECK_EQUAL(derivative, tools::max_value<double>() / 2);
    BOOST_CHECK_EQUAL(result, 0);
 }
 

From 0d6ff719b1bcf6efc45c8b5f4de6f52c3d1a354a Mon Sep 17 00:00:00 2001
From: Kilian Kilger <kkilger@gmail.com>
Date: Tue, 29 Apr 2025 09:57:58 +0200
Subject: [PATCH 3/3] Incomplete Gamma: Restore previous overflow check

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

diff --git a/include/boost/math/special_functions/gamma.hpp b/include/boost/math/special_functions/gamma.hpp
index 7856fd2360..c49e2b3d46 100644
--- a/include/boost/math/special_functions/gamma.hpp
+++ b/include/boost/math/special_functions/gamma.hpp
@@ -1657,22 +1657,13 @@ BOOST_MATH_GPU_ENABLED T gamma_incomplete_imp_final(T a, T x, bool normalised, b
    }
    if(p_derivative)
    {
-      /*
-      * We can never reach this:
-      if((x < 1) && (tools::max_value<T>() * x < *p_derivative))
+      if((x == 0) || ((x < 1) && (tools::max_value<T>() * x < *p_derivative)))
       {
          // overflow, just return an arbitrarily large value:
          *p_derivative = tools::max_value<T>() / 2;
       }
-      */
-      BOOST_MATH_ASSERT((x >= 1) || (tools::max_value<T>() * x >= *p_derivative));
-      //
-      // Need to convert prefix term to derivative:
-      //
-      if (x > 0)
+      else
          *p_derivative /= x;
-      else 
-         *p_derivative = tools::max_value<T>() / 2;
    }
 
    return result;