From 20fc2addea3415275ce1b32b4e676c4542a8bc6e Mon Sep 17 00:00:00 2001
From: jzmaddock <john@johnmaddock.co.uk>
Date: Sun, 27 Nov 2022 19:14:14 +0000
Subject: [PATCH] More and better handling of ibeta small args.

---
 include/boost/math/special_functions/beta.hpp |  3 ++
 .../detail/ibeta_inverse.hpp                  | 40 +++++++++++++++++--
 .../boost/math/special_functions/gamma.hpp    |  2 +-
 test/test_ibeta_inv.hpp                       | 11 +++++
 4 files changed, 52 insertions(+), 4 deletions(-)

diff --git a/include/boost/math/special_functions/beta.hpp b/include/boost/math/special_functions/beta.hpp
index 4f0a9b61b..f4019bd3d 100644
--- a/include/boost/math/special_functions/beta.hpp
+++ b/include/boost/math/special_functions/beta.hpp
@@ -570,6 +570,9 @@ T ibeta_series(T a, T b, T x, T s0, const Lanczos&, bool normalised, T* p_deriva
       T cgh = static_cast<T>(c + Lanczos::g() - 0.5f);
       result = Lanczos::lanczos_sum_expG_scaled(c) / (Lanczos::lanczos_sum_expG_scaled(a) * Lanczos::lanczos_sum_expG_scaled(b));
 
+      if (!(boost::math::isfinite)(result))
+         result = 0;
+
       T l1 = log(cgh / bgh) * (b - 0.5f);
       T l2 = log(x * cgh / agh) * a;
       //
diff --git a/include/boost/math/special_functions/detail/ibeta_inverse.hpp b/include/boost/math/special_functions/detail/ibeta_inverse.hpp
index 9b7bb50b8..9f257793d 100644
--- a/include/boost/math/special_functions/detail/ibeta_inverse.hpp
+++ b/include/boost/math/special_functions/detail/ibeta_inverse.hpp
@@ -638,6 +638,12 @@ T ibeta_inv_imp(T a, T b, T p, T q, const Policy& pol, T* py)
 #endif
                {
                   bet = boost::math::beta(a, b, pol);
+
+                  typedef typename Policy::overflow_error_type overflow_type;
+
+                  BOOST_IF_CONSTEXPR(overflow_type::value != boost::math::policies::throw_on_error)
+                     if(bet > tools::max_value<T>())
+                        bet = tools::max_value<T>();
                }
 #ifndef BOOST_NO_EXCEPTIONS
                catch (const std::overflow_error&)
@@ -686,7 +692,7 @@ T ibeta_inv_imp(T a, T b, T p, T q, const Policy& pol, T* py)
          invert = !invert;
          xs = 1 - xs;
       }
-      if (a < tools::min_value<T>())
+      if ((a < tools::min_value<T>()) && (b > tools::min_value<T>()))
       {
          if (py)
          {
@@ -715,6 +721,8 @@ T ibeta_inv_imp(T a, T b, T p, T q, const Policy& pol, T* py)
       if (overflow || !(boost::math::isfinite)(bet))
       {
          xg = exp((boost::math::lgamma(a + 1, pol) + boost::math::lgamma(b, pol) - boost::math::lgamma(a + b, pol) + log(p)) / a);
+         if (xg > 2 / tools::epsilon<T>())
+            xg = 2 / tools::epsilon<T>();
       }
       else
          xg = pow(a * p * bet, 1/a);
@@ -814,7 +822,20 @@ T ibeta_inv_imp(T a, T b, T p, T q, const Policy& pol, T* py)
       }
       if(pow(p, 1/a) < 0.5)
       {
-         x = pow(p * a * boost::math::beta(a, b, pol), 1 / a);
+#ifndef BOOST_NO_EXCEPTIONS
+         try 
+         {
+#endif
+            x = pow(p * a * boost::math::beta(a, b, pol), 1 / a);
+            if ((x > 1) || !(boost::math::isfinite)(x))
+               x = 1;
+#ifndef BOOST_NO_EXCEPTIONS
+         }
+         catch (const std::overflow_error&)
+         {
+            x = 1;
+         }
+#endif
          if(x == 0)
             x = boost::math::tools::min_value<T>();
          y = 1 - x;
@@ -822,7 +843,20 @@ T ibeta_inv_imp(T a, T b, T p, T q, const Policy& pol, T* py)
       else /*if(pow(q, 1/b) < 0.1)*/
       {
          // model a distorted quarter circle:
-         y = pow(1 - pow(p, b * boost::math::beta(a, b, pol)), 1/b);
+#ifndef BOOST_NO_EXCEPTIONS
+         try 
+         {
+#endif
+            y = pow(1 - pow(p, b * boost::math::beta(a, b, pol)), 1/b);
+            if ((y > 1) || !(boost::math::isfinite)(y))
+               y = 1;
+#ifndef BOOST_NO_EXCEPTIONS
+         }
+         catch (const std::overflow_error&)
+         {
+            y = 1;
+         }
+#endif
          if(y == 0)
             y = boost::math::tools::min_value<T>();
          x = 1 - y;
diff --git a/include/boost/math/special_functions/gamma.hpp b/include/boost/math/special_functions/gamma.hpp
index f2c5a83eb..68afe66e3 100644
--- a/include/boost/math/special_functions/gamma.hpp
+++ b/include/boost/math/special_functions/gamma.hpp
@@ -882,7 +882,7 @@ T regularised_gamma_prefix(T a, T z, const Policy& pol, const Lanczos& l)
       //
       // TODO: is this still required?  Lanczos approx should be better now?
       //
-      if(z <= tools::log_min_value<T>())
+      if((z <= tools::log_min_value<T>()) || (a < 1 / tools::max_value<T>()))
       {
          // Oh dear, have to use logs, should be free of cancellation errors though:
          return exp(a * log(z) - z - lgamma_imp(a, pol, l));
diff --git a/test/test_ibeta_inv.hpp b/test/test_ibeta_inv.hpp
index 1e1c841f2..9e303e923 100644
--- a/test/test_ibeta_inv.hpp
+++ b/test/test_ibeta_inv.hpp
@@ -305,5 +305,16 @@ void test_spots(T)
       BOOST_MATH_CHECK_THROW(::boost::math::ibeta_inv(static_cast<T>(2.125), -n, static_cast<T>(0.125)), std::domain_error);
       BOOST_MATH_CHECK_THROW(::boost::math::ibeta_inv(static_cast<T>(2.125), static_cast<T>(1.125), -n), std::domain_error);
    }
+   if (std::numeric_limits<T>::has_denorm)
+   {
+      T m = std::numeric_limits<T>::denorm_min();
+      T small = 2 * (std::numeric_limits<T>::min)();
+      BOOST_CHECK((boost::math::isfinite)(boost::math::ibeta_inv(m, static_cast<T>(2.125), static_cast<T>(0.125))));
+      BOOST_CHECK((boost::math::isfinite)(boost::math::ibeta_inv(m, m, static_cast<T>(0.125))));
+      BOOST_CHECK_LT(boost::math::ibeta_inv(m, static_cast<T>(12.125), static_cast<T>(0.125)), small);
+      BOOST_CHECK((boost::math::isfinite)(boost::math::ibeta_inv(static_cast<T>(2.125), m, static_cast<T>(0.125))));
+      BOOST_CHECK((boost::math::isfinite)(boost::math::ibeta_inv(static_cast<T>(12.125), m, static_cast<T>(0.125))));
+      BOOST_CHECK((boost::math::isfinite)(boost::math::ibeta_inv(m, m, static_cast<T>(0.125))));
+   }
 }