Merge pull request #1259 from boostorg/scipy21725

Fix for inverse ibeta with large a,b.

Author: jzmaddock

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

@@ -683,11 +683,17 @@ BOOST_MATH_GPU_ENABLED T ibeta_inv_imp(T a, T b, T p, T q, const Policy& pol, T*
             }
             else
                y = 1;
-            if(y > 1e-5)
+            if((y > 1e-5) && (std::min)(a, b) < 1000)
             {
                x = temme_method_3_ibeta_inverse(a, b, p, q, pol);
                y = 1 - x;
             }
+            else if ((y > 1e-5) && (std::min)(a, b) > 1000)
+            {
+               // All options have failed, use the saddle point as a starting location:
+               x = (std::max)(a, b) / (a + b);
+               y = (std::min)(a, b) / (a + b);
+            }
          }
       }
    }

test/test_ibeta_inv.hpp

@@ -320,5 +320,14 @@ void test_spots(T)
       BOOST_CHECK((boost::math::isfinite)(boost::math::ibeta_inv(m, m, static_cast<T>(0.125))));
    }
    #endif
+   //
+   // scipy: https://github.com/scipy/scipy/issues/21725
+   //
+   BOOST_CHECK_CLOSE(
+      ::boost::math::ibeta_inv(
+         static_cast<T>(1.0e11L),
+         static_cast<T>(1.0e13L),
+         static_cast<T>(0.995L)),
+      static_cast<T>(0.0099010703473402885173268397418009652L), 5e-10);
 }