Fixes ibeta for large arguments

include/boost/math/special_functions/beta.hpp

@@ -19,7 +19,6 @@
 #include <boost/math/tools/cstdint.hpp>
 #include <boost/math/tools/tuple.hpp>
 #include <boost/math/tools/promotion.hpp>
-#include <boost/math/tools/cstdint.hpp>
 #include <boost/math/special_functions/gamma.hpp>
 #include <boost/math/special_functions/erf.hpp>
 #include <boost/math/special_functions/log1p.hpp>
@@ -475,7 +474,6 @@ BOOST_MATH_GPU_ENABLED T ibeta_power_terms(T a,
                         const char* = "boost::math::ibeta<%1%>(%1%, %1%, %1%)")
 {
    BOOST_MATH_STD_USING
-
    if(!normalised)
    {
       return prefix * pow(x, a) * pow(y, b);
@@ -818,25 +816,25 @@ struct ibeta_fraction2_t
 {
    typedef boost::math::pair<T, T> result_type;
 
-   BOOST_MATH_GPU_ENABLED ibeta_fraction2_t(T a_, T b_, T x_, T y_) : a(a_), b(b_), x(x_), y(y_), m(0) {}
+   BOOST_MATH_GPU_ENABLED ibeta_fraction2_t(T a_, T b_, T x_, T y_) : a(a_), b(b_), x(x_), y(y_), m(static_cast<T>(0)) {}
 
    BOOST_MATH_GPU_ENABLED result_type operator()()
    {
       T denom = (a + 2 * m - 1);
       T aN = (m * (a + m - 1) / denom) * ((a + b + m - 1) / denom) * (b - m) * x * x;
 
-      T bN = static_cast<T>(m);
+      T bN = m;
       bN += (m * (b - m) * x) / (a + 2*m - 1);
       bN += ((a + m) * (a * y - b * x + 1 + m *(2 - x))) / (a + 2*m + 1);
 
-      ++m;
+      m += 1;
 
       return boost::math::make_pair(aN, bN);
    }
 
 private:
    T a, b, x, y;
-   int m;
+   T m;
 };
 //
 // Evaluate the incomplete beta via the continued fraction representation:
@@ -1142,7 +1140,20 @@ BOOST_MATH_GPU_ENABLED T ibeta_large_ab(T a, T b, T x, T y, bool invert, bool no
 
    T x0 = a / (a + b);
    T y0 = b / (a + b);
-   T nu = x0 * log(x / x0) + y0 * log(y / y0);
+
+   // Expand nu about x0
+   T nu = 0;
+   for (int i=2; i<5; i++)
+   {
+      nu += pow(x-x0, i) / i * (pow(x0, -(i-1)) - pow(x0-1, -(i-1))) * pow(-1, i+1);
+   }
+   // Calculate the next term in the series
+   T remainder = pow(x-x0, 5) / 5 * (pow(x0, -4) - pow(x0-1, 4));
+
+   // If the remainder is large, then fall back to using the log formula
+   if (remainder >= tools::forth_root_epsilon<T>()){
+      nu = x0 * log(x / x0) + y0 * log(y / y0);
+   }
    //
    // Above compution is unstable, force nu to zero if
    // something went wrong:
@@ -1181,7 +1192,7 @@ BOOST_MATH_GPU_ENABLED T ibeta_large_ab(T a, T b, T x, T y, bool invert, bool no
    T mul = 1;
    if (!normalised)
       mul = boost::math::beta(a, b, pol);
-
+   // T log_erf_remainder = -0.5 * log(2 * constants::pi<T>() * (a+b)) + a * log(x / x0) + b * log((1-x) / (1-x0)) + log(abs(1/nu - sqrt(x0 * (1-x0)) / (x-x0)));
    return mul * ((invert ? (1 + boost::math::erf(-nu * sqrt((a + b) / 2), pol)) / 2 : boost::math::erfc(-nu * sqrt((a + b) / 2), pol) / 2));
 }
 
@@ -1571,37 +1582,47 @@ BOOST_MATH_GPU_ENABLED T ibeta_imp(T a, T b, T x, const Policy& pol, bool inv, b
       else
       {
          // a and b both large:
-         bool use_asym = false;
          T ma = BOOST_MATH_GPU_SAFE_MAX(a, b);
-         T xa = ma == a ? x : y;
          T saddle = ma / (a + b);
-         T powers = 0;
-         if ((ma > 1e-5f / tools::epsilon<T>()) && (ma / BOOST_MATH_GPU_SAFE_MIN(a, b) < (xa < saddle ? 2 : 15)))
+         // If ma large and x is close to saddle, use `erf` approximation in `ibeta_large_ab`. 
+         // In this case we don't need to invert
+         if ((ma > 1e-5f / tools::epsilon<T>()) && (fabs(saddle - x) < 1e-12) && (1 - saddle > 0.125))
          {
-            if (a == b)
-               use_asym = true;
-            else
-            {
-               powers = exp(log(x / (a / (a + b))) * a + log(y / (b / (a + b))) * b);
-               if (powers < tools::epsilon<T>())
-                  use_asym = true;
-            }
+            fract = ibeta_large_ab(a, b, x, y, invert, normalised, pol);
+            invert = false;
          }
-         if(use_asym)
+         else
          {
-            fract = ibeta_large_ab(a, b, x, y, invert, normalised, pol);
-            if (fract * tools::epsilon<T>() < powers)
+            // This is the same logic that is in `ibeta_fraction2`. We have repeated 
+            // the implementation here because when x is close to saddle, the continued
+            // fraction will not converge. If this occurs, we fall back to the `erf`
+            // approximation. Simply using `ibeta_fraction2` will cause an evaluation 
+            // error to be thrown and we won't be able to use the `erf` approximation.
+            typedef typename lanczos::lanczos<T, Policy>::type lanczos_type;
+            T local_result = ibeta_power_terms(a, b, x, y, lanczos_type(), normalised, pol);
+            if (p_derivative)
             {
-               // Erf approximation failed, correction term is too large, fall back:
-               fract = ibeta_fraction2(a, b, x, y, pol, normalised, p_derivative);
+               *p_derivative = local_result;
+               BOOST_MATH_ASSERT(*p_derivative >= 0);
+            }
+            if (local_result != 0)
+            {
+               ibeta_fraction2_t<T> f(a, b, x, y);
+               boost::math::uintmax_t max_terms = boost::math::policies::get_max_series_iterations<Policy>();
+               T local_fract = boost::math::tools::continued_fraction_b(f, boost::math::policies::get_epsilon<T, Policy>(), max_terms);
+               if (max_terms >= boost::math::policies::get_max_series_iterations<Policy>())
+               {
+                  // Continued fraction failed, fall back to asymptotic expansion:
+                  fract = ibeta_large_ab(a, b, x, y, invert, normalised, pol);
+                  invert = false;
+               }
+               else{
+                  fract = local_result / local_fract;
+               }
             }
             else
-               invert = false;
+               fract = 0;
          }
-         else
-            fract = ibeta_fraction2(a, b, x, y, pol, normalised, p_derivative);
-
-         BOOST_MATH_INSTRUMENT_VARIABLE(fract);
       }
    }
    if(p_derivative)

test/test_ibeta.cpp

@@ -380,18 +380,18 @@ BOOST_AUTO_TEST_CASE( test_main )
    gsl_set_error_handler_off();
 #endif
 #ifdef TEST_FLOAT
-   test_spots(0.0F);
+   test_spots(0.0F, "float");
 #endif
 #ifdef TEST_DOUBLE
-   test_spots(0.0);
+   test_spots(0.0, "double");
 #endif
 #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
 #ifdef TEST_LDOUBLE
-   test_spots(0.0L);
+   test_spots(0.0L, "long double");
 #endif
 #if !BOOST_WORKAROUND(BOOST_BORLANDC, BOOST_TESTED_AT(0x582))
 #if defined(TEST_REAL_CONCEPT) && !defined(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
 #endif

test/test_ibeta.hpp

@@ -152,12 +152,13 @@ void test_beta(T, const char* name)
 }
 
 template <class T>
-void test_spots(T)
+void test_spots(T, const char* name)
 {
    //
    // basic sanity checks, tolerance is 30 epsilon expressed as a percentage:
    // Spot values are from http://functions.wolfram.com/webMathematica/FunctionEvaluation.jsp?name=BetaRegularized
    // using precision of 50 decimal digits.
+   std::cout << "Testing spot values with type " << name << std::endl;
    T tolerance = boost::math::tools::epsilon<T>() * 3000;
    if (boost::math::tools::digits<T>() > 100)
       tolerance *= 2;
@@ -487,5 +488,35 @@ void test_spots(T)
       }
       BOOST_CHECK_EQUAL(boost::math::ibeta(static_cast<T>(2), static_cast<T>(1), static_cast<T>(0)), 0);
       BOOST_CHECK_EQUAL(boost::math::ibeta(static_cast<T>(1), static_cast<T>(2), static_cast<T>(0)), 0);
+
+      // Bug testing for large a,b and x close to a / (a+b). See PR 1363.
+      // The values for a,b are just too large for floats to handle. The tests here only
+      // check if ibeta is monotonically increasing for a,b fixed and x increasing. Both
+      // Mathematica and mpmath (in Python) are not able to evaluate ibeta for such large
+      // values of a,b so spot testing wasn't possible. The accuracy of the fix for 
+      // these large values is very sensitive to the computer architecture. 
+      // Macos with arm64 does the worst but linux and windows pass a larger range of tests. 
+      if (!std::is_same<T, float>::value)
+      {
+         // Larger values of a,b become more numerically unstable. Larger/smaller values
+         // of delta also become unstable.
+         T a_values[1] = {10000000272564224};
+         T b_values[1] = {9965820922822656};
+         T delta[7] = {0, 1e-15, 1e-14, 1e-13, 1e-12, 1e-11, 1e-10};
+         T a, b, x;
+         for (unsigned int i=0; i<1; i++){
+            a = a_values[i];
+            b = b_values[i];
+            x = a / (a+b); // roughly the median of ibeta
+
+            for (unsigned int j=0; j < 6; j++)
+            {
+               BOOST_CHECK_MESSAGE(boost::math::ibeta(a, b, x + delta[j+1]) > boost::math::ibeta(a, b, x + delta[j]), 
+                  "ibeta not monotonically increasing above a/(a+b) for ibeta(" << a << ", " << b << ", " << x << ") and delta=" << delta[j] << ": [" << boost::math::ibeta(a, b, x + delta[j+1]) << " >= " << boost::math::ibeta(a, b, x + delta[j]) << "]");
+               BOOST_CHECK_MESSAGE(boost::math::ibeta(a, b, x - delta[j]) > boost::math::ibeta(a, b, x - delta[j+1]), 
+                  "ibeta not monotonically increasing below a/(a+b) for ibeta(" << a << ", " << b << ", " << x << ") and delta=" << delta[j] << ": [" << boost::math::ibeta(a, b, x - delta[j]) << " >= " << boost::math::ibeta(a, b, x - delta[j+1]) << "]");  
+            }
+         }
+      }
 }