diff --git a/include/boost/math/special_functions/beta.hpp b/include/boost/math/special_functions/beta.hpp index f6327fa220..3f31bea408 100644 --- a/include/boost/math/special_functions/beta.hpp +++ b/include/boost/math/special_functions/beta.hpp @@ -475,7 +475,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,14 +817,14 @@ struct ibeta_fraction2_t { typedef boost::math::pair 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(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(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); @@ -836,7 +835,7 @@ struct ibeta_fraction2_t private: T a, b, x, y; - int m; + T m; }; // // Evaluate the incomplete beta via the continued fraction representation: @@ -1142,7 +1141,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()){ + nu = x0 * log(x / x0) + y0 * log(y / y0); + } // // Above compution is unstable, force nu to zero if // something went wrong: @@ -1181,7 +1193,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() * (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 +1583,48 @@ 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()) && (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()) && (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()) - 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() < 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::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 f(a, b, x, y); + ibeta_fraction2_t g(a, b, x, y); + boost::math::uintmax_t max_terms = boost::math::policies::get_max_series_iterations(); + T local_fract = boost::math::tools::continued_fraction_b(f, boost::math::policies::get_epsilon(), max_terms); + if (max_terms >= boost::math::policies::get_max_series_iterations()) + { + // 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) diff --git a/test/test_ibeta.cpp b/test/test_ibeta.cpp index 563844a06c..0059b46619 100644 --- a/test/test_ibeta.cpp +++ b/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 diff --git a/test/test_ibeta.hpp b/test/test_ibeta.hpp index 2f0ea47107..13ff157f9d 100644 --- a/test/test_ibeta.hpp +++ b/test/test_ibeta.hpp @@ -152,12 +152,13 @@ void test_beta(T, const char* name) } template -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() * 3000; if (boost::math::tools::digits() > 100) tolerance *= 2; @@ -487,5 +488,35 @@ void test_spots(T) } BOOST_CHECK_EQUAL(boost::math::ibeta(static_cast(2), static_cast(1), static_cast(0)), 0); BOOST_CHECK_EQUAL(boost::math::ibeta(static_cast(1), static_cast(2), static_cast(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::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]) << "]"); + } + } + } }