Merge pull request #1026 from boostorg/issue1006

Avoid spurious overflow and divide by zero in ibeta.

Author: jzmaddock

include/boost/math/special_functions/beta.hpp

@@ -230,7 +230,10 @@ T ibeta_power_terms(T a,
    T agh = static_cast<T>(a + Lanczos::g() - 0.5f);
    T bgh = static_cast<T>(b + Lanczos::g() - 0.5f);
    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 ((a < tools::min_value<T>()) || (b < tools::min_value<T>()))
+      result = 0;  // denominator overflows in this case
+   else
+      result = Lanczos::lanczos_sum_expG_scaled(c) / (Lanczos::lanczos_sum_expG_scaled(a) * Lanczos::lanczos_sum_expG_scaled(b));
    result *= prefix;
    // combine with the leftover terms from the Lanczos approximation:
    result *= sqrt(bgh / boost::math::constants::e<T>());
@@ -389,14 +392,15 @@ T ibeta_power_terms(T a,
          }
          else
          {
-            T p1 = pow(b1, a / b);
+            // This protects against spurious overflow in a/b:
+            T p1 = (b1 < 1) && (b < 1) && (tools::max_value<T>() * b < a) ? static_cast<T>(0) : static_cast<T>(pow(b1, a / b));
             T l3 = (p1 != 0) && (b2 != 0) ? (log(p1) + log(b2)) * b : tools::max_value<T>();  // arbitrary large value if the logs would fail!
             if((l3 < tools::log_max_value<T>())
                && (l3 > tools::log_min_value<T>()))
             {
                result *= pow(p1 * b2, b);
             }
-            else
+            else if(result != 0)  // we can elude the calculation below if we're already going to be zero
             {
                l2 += l1 + log(result);
                if(l2 >= tools::log_max_value<T>())
@@ -656,7 +660,10 @@ T ibeta_series(T a, T b, T x, T s0, const Lanczos&, bool normalised, T* p_deriva
       T agh = static_cast<T>(a + Lanczos::g() - 0.5f);
       T bgh = static_cast<T>(b + Lanczos::g() - 0.5f);
       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 ((a < tools::min_value<T>()) || (b < tools::min_value<T>()))
+         result = 0;  // denorms cause overflow in the Lanzos series, result will be zero anyway
+      else
+         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;
@@ -689,10 +696,13 @@ T ibeta_series(T a, T b, T x, T s0, const Lanczos&, bool normalised, T* p_deriva
          //
          // Oh dear, we need logs, and this *will* cancel:
          //
-         result = log(result) + l1 + l2 + (log(agh) - 1) / 2;
-         if(p_derivative)
-            *p_derivative = exp(result + b * log(y));
-         result = exp(result);
+         if (result != 0)  // elude calculation when result will be zero.
+         {
+            result = log(result) + l1 + l2 + (log(agh) - 1) / 2;
+            if (p_derivative)
+               *p_derivative = exp(result + b * log(y));
+            result = exp(result);
+         }
       }
    }
    else

include/boost/math/special_functions/detail/ibeta_inverse.hpp

@@ -825,7 +825,21 @@ T ibeta_inv_imp(T a, T b, T p, T q, const Policy& pol, T* py)
          std::swap(p, q);
          invert = !invert;
       }
-      if(pow(p, 1/a) < 0.5)
+      if (a < tools::min_value<T>())
+      {
+         // Avoid spurious overflows for denorms:
+         if (p < 1)
+         {
+            x = 1;
+            y = 0;
+         }
+         else
+         {
+            x = 0;
+            y = 1;
+         }
+      }
+      else if(pow(p, 1/a) < 0.5)
       {
 #ifndef BOOST_NO_EXCEPTIONS
          try 

include/boost/math/tools/roots.hpp

@@ -589,7 +589,8 @@ namespace detail {
    #ifdef BOOST_MATH_INSTRUMENT
          std::cout << "Second order root iteration, delta = " << delta << ", residual = " << f0 << "\n";
    #endif
-         T convergence = fabs(delta / delta2);
+         // We need to avoid delta/delta2 overflowing here:
+         T convergence = (fabs(delta2) > 1) || (fabs(tools::max_value<T>() * delta2) > fabs(delta)) ? fabs(delta / delta2) : tools::max_value<T>();
          if ((convergence > 0.8) && (convergence < 2))
          {
             // last two steps haven't converged.

test/Jamfile.v2

@@ -169,6 +169,7 @@ test-suite special_fun :
    [ run git_issue_810.cpp ../../test/build//boost_unit_test_framework  ]
    [ run git_issue_826.cpp ../../test/build//boost_unit_test_framework  ]
    [ run git_issue_961.cpp ]
+   [ run git_issue_1006.cpp ]
    [ run special_functions_test.cpp ../../test/build//boost_unit_test_framework  ]
    [ run test_airy.cpp test_instances//test_instances pch_light ../../test/build//boost_unit_test_framework  ]
    [ run test_bessel_j.cpp test_instances//test_instances pch_light ../../test/build//boost_unit_test_framework  ]

test/git_issue_1006.cpp

@@ -0,0 +1,60 @@
+//  (C) Copyright Matt Borland 2023.
+//  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)
+
+#include "math_unit_test.hpp"
+#include <cfenv>
+#include <iostream>
+#include <boost/math/distributions/beta.hpp>
+
+#pragma STDC FENV_ACCESS ON
+
+// Show and then clear the fenv flags
+void show_fpexcept_flags()
+{
+   bool fe = false;
+
+   if (std::fetestexcept(FE_OVERFLOW))
+   {
+      fe = true;
+      std::cerr << "FE_OVERFLOW" << std::endl;
+   }
+   if (std::fetestexcept(FE_UNDERFLOW))
+   {
+      //fe = true;
+      std::cerr << "FE_UNDERFLOW" << std::endl;
+   }
+   if (std::fetestexcept(FE_DIVBYZERO))
+   {
+      fe = true;
+      std::cerr << "FE_DIVBYZERO" << std::endl;
+   }
+   if (std::fetestexcept(FE_INVALID))
+   {
+      fe = true;
+      std::cerr << "FE_INVALID" << std::endl;
+   }
+
+   CHECK_EQUAL(fe, false);
+
+   std::feclearexcept(FE_ALL_EXCEPT);
+}
+
+int main()
+{
+   // Default Scipy policy
+   using my_policy = boost::math::policies::policy<boost::math::policies::promote_float<false>, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy>;
+
+   double a = 1e-307;
+
+   while (a)
+   {
+      const auto dist_a = boost::math::beta_distribution<double, my_policy>(1e-308, 5.0);
+      const auto dist_a_ppf = boost::math::quantile(dist_a, 0.2);
+      show_fpexcept_flags();
+      CHECK_MOLLIFIED_CLOSE(dist_a_ppf, 0.0, 1e-10);
+      a /= 2;
+   }
+   return boost::math::test::report_errors();
+}