boost/math/special_functions/beta.hpp
// (C) Copyright John Maddock 2006. // (C) Copyright Matt Borland 2024. // 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) #ifndef BOOST_MATH_SPECIAL_BETA_HPP #define BOOST_MATH_SPECIAL_BETA_HPP #ifdef _MSC_VER #pragma once #endif #include <boost/math/tools/config.hpp> #include <boost/math/tools/type_traits.hpp> #include <boost/math/tools/assert.hpp> #include <boost/math/tools/precision.hpp> #include <boost/math/tools/numeric_limits.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> #include <boost/math/special_functions/expm1.hpp> #include <boost/math/special_functions/trunc.hpp> #include <boost/math/special_functions/lanczos.hpp> #include <boost/math/policies/policy.hpp> #include <boost/math/policies/error_handling.hpp> #include <boost/math/constants/constants.hpp> #include <boost/math/special_functions/math_fwd.hpp> #include <boost/math/special_functions/binomial.hpp> #include <boost/math/special_functions/factorials.hpp> #include <boost/math/tools/roots.hpp> namespace boost{ namespace math{ namespace detail{ // // Implementation of Beta(a,b) using the Lanczos approximation: // template <class T, class Lanczos, class Policy> BOOST_MATH_GPU_ENABLED T beta_imp(T a, T b, const Lanczos&, const Policy& pol) { BOOST_MATH_STD_USING // for ADL of std names if(a <= 0) return policies::raise_domain_error<T>("boost::math::beta<%1%>(%1%,%1%)", "The arguments to the beta function must be greater than zero (got a=%1%).", a, pol); if(b <= 0) return policies::raise_domain_error<T>("boost::math::beta<%1%>(%1%,%1%)", "The arguments to the beta function must be greater than zero (got b=%1%).", b, pol); T result; // LCOV_EXCL_LINE T prefix = 1; T c = a + b; // Special cases: if((c == a) && (b < tools::epsilon<T>())) return 1 / b; else if((c == b) && (a < tools::epsilon<T>())) return 1 / a; if(b == 1) return 1/a; else if(a == 1) return 1/b; else if(c < tools::epsilon<T>()) { result = c / a; result /= b; return result; } /* // // This code appears to be no longer necessary: it was // used to offset errors introduced from the Lanczos // approximation, but the current Lanczos approximations // are sufficiently accurate for all z that we can ditch // this. It remains in the file for future reference... // // If a or b are less than 1, shift to greater than 1: if(a < 1) { prefix *= c / a; c += 1; a += 1; } if(b < 1) { prefix *= c / b; c += 1; b += 1; } */ if(a < b) { BOOST_MATH_GPU_SAFE_SWAP(a, b); } // Lanczos calculation: 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(a) * (Lanczos::lanczos_sum_expG_scaled(b) / Lanczos::lanczos_sum_expG_scaled(c)); T ambh = a - 0.5f - b; if((fabs(b * ambh) < (cgh * 100)) && (a > 100)) { // Special case where the base of the power term is close to 1 // compute (1+x)^y instead: result *= exp(ambh * boost::math::log1p(-b / cgh, pol)); } else { result *= pow(agh / cgh, a - T(0.5) - b); } if(cgh > 1e10f) // this avoids possible overflow, but appears to be marginally less accurate: result *= pow((agh / cgh) * (bgh / cgh), b); else result *= pow((agh * bgh) / (cgh * cgh), b); result *= sqrt(boost::math::constants::e<T>() / bgh); // If a and b were originally less than 1 we need to scale the result: result *= prefix; return result; } // template <class T, class Lanczos> beta_imp(T a, T b, const Lanczos&) // // Generic implementation of Beta(a,b) without Lanczos approximation support // (Caution this is slow!!!): // #ifndef BOOST_MATH_HAS_GPU_SUPPORT template <class T, class Policy> BOOST_MATH_GPU_ENABLED T beta_imp(T a, T b, const lanczos::undefined_lanczos& l, const Policy& pol) { BOOST_MATH_STD_USING if(a <= 0) return policies::raise_domain_error<T>("boost::math::beta<%1%>(%1%,%1%)", "The arguments to the beta function must be greater than zero (got a=%1%).", a, pol); if(b <= 0) return policies::raise_domain_error<T>("boost::math::beta<%1%>(%1%,%1%)", "The arguments to the beta function must be greater than zero (got b=%1%).", b, pol); const T c = a + b; // Special cases: if ((c == a) && (b < tools::epsilon<T>())) return 1 / b; else if ((c == b) && (a < tools::epsilon<T>())) return 1 / a; if (b == 1) return 1 / a; else if (a == 1) return 1 / b; else if (c < tools::epsilon<T>()) { T result = c / a; result /= b; return result; } // Regular cases start here: const T min_sterling = minimum_argument_for_bernoulli_recursion<T>(); long shift_a = 0; long shift_b = 0; if(a < min_sterling) shift_a = 1 + ltrunc(min_sterling - a); if(b < min_sterling) shift_b = 1 + ltrunc(min_sterling - b); long shift_c = shift_a + shift_b; if ((shift_a == 0) && (shift_b == 0)) { return pow(a / c, a) * pow(b / c, b) * scaled_tgamma_no_lanczos(a, pol) * scaled_tgamma_no_lanczos(b, pol) / scaled_tgamma_no_lanczos(c, pol); } else if ((a < 1) && (b < 1)) { return boost::math::tgamma(a, pol) * (boost::math::tgamma(b, pol) / boost::math::tgamma(c)); } else if(a < 1) return boost::math::tgamma(a, pol) * boost::math::tgamma_delta_ratio(b, a, pol); else if(b < 1) return boost::math::tgamma(b, pol) * boost::math::tgamma_delta_ratio(a, b, pol); else { T result = beta_imp(T(a + shift_a), T(b + shift_b), l, pol); // // Recursion: // for (long i = 0; i < shift_c; ++i) { result *= c + i; if (i < shift_a) result /= a + i; if (i < shift_b) result /= b + i; } return result; } } // template <class T>T beta_imp(T a, T b, const lanczos::undefined_lanczos& l) #endif // // Compute the leading power terms in the incomplete Beta: // // (x^a)(y^b)/Beta(a,b) when normalised, and // (x^a)(y^b) otherwise. // // Almost all of the error in the incomplete beta comes from this // function: particularly when a and b are large. Computing large // powers are *hard* though, and using logarithms just leads to // horrendous cancellation errors. // template <class T, class Lanczos, class Policy> BOOST_MATH_GPU_ENABLED T ibeta_power_terms(T a, T b, T x, T y, const Lanczos&, bool normalised, const Policy& pol, T prefix = 1, const char* function = "boost::math::ibeta<%1%>(%1%, %1%, %1%)") { BOOST_MATH_STD_USING if(!normalised) { // can we do better here? return pow(x, a) * pow(y, b); } T result; // LCOV_EXCL_LINE T c = a + b; // combine power terms with Lanczos approximation: 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); 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>()); result *= sqrt(agh / cgh); // l1 and l2 are the base of the exponents minus one: T l1 = (x * b - y * agh) / agh; T l2 = (y * a - x * bgh) / bgh; if((BOOST_MATH_GPU_SAFE_MIN(fabs(l1), fabs(l2)) < 0.2)) { // when the base of the exponent is very near 1 we get really // gross errors unless extra care is taken: if((l1 * l2 > 0) || (BOOST_MATH_GPU_SAFE_MIN(a, b) < 1)) { // // This first branch handles the simple cases where either: // // * The two power terms both go in the same direction // (towards zero or towards infinity). In this case if either // term overflows or underflows, then the product of the two must // do so also. // *Alternatively if one exponent is less than one, then we // can't productively use it to eliminate overflow or underflow // from the other term. Problems with spurious overflow/underflow // can't be ruled out in this case, but it is *very* unlikely // since one of the power terms will evaluate to a number close to 1. // if(fabs(l1) < 0.1) { result *= exp(a * boost::math::log1p(l1, pol)); BOOST_MATH_INSTRUMENT_VARIABLE(result); } else { result *= pow((x * cgh) / agh, a); BOOST_MATH_INSTRUMENT_VARIABLE(result); } if(fabs(l2) < 0.1) { result *= exp(b * boost::math::log1p(l2, pol)); BOOST_MATH_INSTRUMENT_VARIABLE(result); } else { result *= pow((y * cgh) / bgh, b); BOOST_MATH_INSTRUMENT_VARIABLE(result); } } else if(BOOST_MATH_GPU_SAFE_MAX(fabs(l1), fabs(l2)) < 0.5) { // // Both exponents are near one and both the exponents are // greater than one and further these two // power terms tend in opposite directions (one towards zero, // the other towards infinity), so we have to combine the terms // to avoid any risk of overflow or underflow. // // We do this by moving one power term inside the other, we have: // // (1 + l1)^a * (1 + l2)^b // = ((1 + l1)*(1 + l2)^(b/a))^a // = (1 + l1 + l3 + l1*l3)^a ; l3 = (1 + l2)^(b/a) - 1 // = exp((b/a) * log(1 + l2)) - 1 // // The tricky bit is deciding which term to move inside :-) // By preference we move the larger term inside, so that the // size of the largest exponent is reduced. However, that can // only be done as long as l3 (see above) is also small. // bool small_a = a < b; T ratio = b / a; if((small_a && (ratio * l2 < 0.1)) || (!small_a && (l1 / ratio > 0.1))) { T l3 = boost::math::expm1(ratio * boost::math::log1p(l2, pol), pol); l3 = l1 + l3 + l3 * l1; l3 = a * boost::math::log1p(l3, pol); result *= exp(l3); BOOST_MATH_INSTRUMENT_VARIABLE(result); } else { T l3 = boost::math::expm1(boost::math::log1p(l1, pol) / ratio, pol); l3 = l2 + l3 + l3 * l2; l3 = b * boost::math::log1p(l3, pol); result *= exp(l3); BOOST_MATH_INSTRUMENT_VARIABLE(result); } } else if(fabs(l1) < fabs(l2)) { // First base near 1 only: T l = a * boost::math::log1p(l1, pol) + b * log((y * cgh) / bgh); if((l <= tools::log_min_value<T>()) || (l >= tools::log_max_value<T>())) { l += log(result); if(l >= tools::log_max_value<T>()) return policies::raise_overflow_error<T>(function, nullptr, pol); // LCOV_EXCL_LINE we can probably never get here, probably! result = exp(l); } else result *= exp(l); BOOST_MATH_INSTRUMENT_VARIABLE(result); } else { // Second base near 1 only: T l = b * boost::math::log1p(l2, pol) + a * log((x * cgh) / agh); if((l <= tools::log_min_value<T>()) || (l >= tools::log_max_value<T>())) { l += log(result); if(l >= tools::log_max_value<T>()) return policies::raise_overflow_error<T>(function, nullptr, pol); // LCOV_EXCL_LINE we can probably never get here, probably! result = exp(l); } else result *= exp(l); BOOST_MATH_INSTRUMENT_VARIABLE(result); } } else { // general case: T b1 = (x * cgh) / agh; T b2 = (y * cgh) / bgh; l1 = a * log(b1); l2 = b * log(b2); BOOST_MATH_INSTRUMENT_VARIABLE(b1); BOOST_MATH_INSTRUMENT_VARIABLE(b2); BOOST_MATH_INSTRUMENT_VARIABLE(l1); BOOST_MATH_INSTRUMENT_VARIABLE(l2); if((l1 >= tools::log_max_value<T>()) || (l1 <= tools::log_min_value<T>()) || (l2 >= tools::log_max_value<T>()) || (l2 <= tools::log_min_value<T>()) ) { // Oops, under/overflow, sidestep if we can: if(a < b) { T p1 = pow(b2, b / a); T l3 = (b1 != 0) && (p1 != 0) ? (a * (log(b1) + log(p1))) : 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 * b1, a); } else { l2 += l1 + log(result); if(l2 >= tools::log_max_value<T>()) return policies::raise_overflow_error<T>(function, nullptr, pol); // LCOV_EXCL_LINE we can probably never get here, probably! result = exp(l2); } } else { // 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 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>()) return policies::raise_overflow_error<T>(function, nullptr, pol); // LCOV_EXCL_LINE we can probably never get here, probably! result = exp(l2); } } BOOST_MATH_INSTRUMENT_VARIABLE(result); } else { // finally the normal case: result *= pow(b1, a) * pow(b2, b); BOOST_MATH_INSTRUMENT_VARIABLE(result); } } BOOST_MATH_INSTRUMENT_VARIABLE(result); if (0 == result) { if ((a > 1) && (x == 0)) return result; // true zero LCOV_EXCL_LINE we can probably never get here if ((b > 1) && (y == 0)) return result; // true zero LCOV_EXCL_LINE we can probably never get here return boost::math::policies::raise_underflow_error<T>(function, nullptr, pol); } return result; } // // Compute the leading power terms in the incomplete Beta: // // (x^a)(y^b)/Beta(a,b) when normalised, and // (x^a)(y^b) otherwise. // // Almost all of the error in the incomplete beta comes from this // function: particularly when a and b are large. Computing large // powers are *hard* though, and using logarithms just leads to // horrendous cancellation errors. // // This version is generic, slow, and does not use the Lanczos approximation. // #ifndef BOOST_MATH_HAS_GPU_SUPPORT template <class T, class Policy> BOOST_MATH_GPU_ENABLED T ibeta_power_terms(T a, T b, T x, T y, const boost::math::lanczos::undefined_lanczos& l, bool normalised, const Policy& pol, T prefix = 1, const char* = "boost::math::ibeta<%1%>(%1%, %1%, %1%)") { BOOST_MATH_STD_USING if(!normalised) { return prefix * pow(x, a) * pow(y, b); } T c = a + b; const T min_sterling = minimum_argument_for_bernoulli_recursion<T>(); long shift_a = 0; long shift_b = 0; if (a < min_sterling) shift_a = 1 + ltrunc(min_sterling - a); if (b < min_sterling) shift_b = 1 + ltrunc(min_sterling - b); if ((shift_a == 0) && (shift_b == 0)) { T power1, power2; bool need_logs = false; if (a < b) { BOOST_MATH_IF_CONSTEXPR(boost::math::numeric_limits<T>::has_infinity) { power1 = pow((x * y * c * c) / (a * b), a); power2 = pow((y * c) / b, b - a); } else { // We calculate these logs purely so we can check for overflow in the power functions T l1 = log((x * y * c * c) / (a * b)); T l2 = log((y * c) / b); if ((l1 * a > tools::log_min_value<T>()) && (l1 * a < tools::log_max_value<T>()) && (l2 * (b - a) < tools::log_max_value<T>()) && (l2 * (b - a) > tools::log_min_value<T>())) { power1 = pow((x * y * c * c) / (a * b), a); power2 = pow((y * c) / b, b - a); } else { need_logs = true; } } } else { BOOST_MATH_IF_CONSTEXPR(boost::math::numeric_limits<T>::has_infinity) { power1 = pow((x * y * c * c) / (a * b), b); power2 = pow((x * c) / a, a - b); } else { // We calculate these logs purely so we can check for overflow in the power functions T l1 = log((x * y * c * c) / (a * b)) * b; T l2 = log((x * c) / a) * (a - b); if ((l1 * a > tools::log_min_value<T>()) && (l1 * a < tools::log_max_value<T>()) && (l2 * (b - a) < tools::log_max_value<T>()) && (l2 * (b - a) > tools::log_min_value<T>())) { power1 = pow((x * y * c * c) / (a * b), b); power2 = pow((x * c) / a, a - b); } else need_logs = true; } } BOOST_MATH_IF_CONSTEXPR(boost::math::numeric_limits<T>::has_infinity) { if (!(boost::math::isnormal)(power1) || !(boost::math::isnormal)(power2)) { need_logs = true; } } if (need_logs) { // // We want: // // (xc / a)^a (yc / b)^b // // But we know that one or other term will over / underflow and combining the logs will be next to useless as that will cause significant cancellation. // If we assume b > a and express z ^ b as(z ^ b / a) ^ a with z = (yc / b) then we can move one power term inside the other : // // ((xc / a) * (yc / b)^(b / a))^a // // However, we're not quite there yet, as the term being exponentiated is quite likely to be close to unity, so let: // // xc / a = 1 + (xb - ya) / a // // analogously let : // // 1 + p = (yc / b) ^ (b / a) = 1 + expm1((b / a) * log1p((ya - xb) / b)) // // so putting the two together we have : // // exp(a * log1p((xb - ya) / a + p + p(xb - ya) / a)) // // Analogously, when a > b we can just swap all the terms around. // // Finally, there are a few cases (x or y is unity) when the above logic can't be used // or where there is no logarithmic cancellation and accuracy is better just using // the regular formula: // T xc_a = x * c / a; T yc_b = y * c / b; if ((x == 1) || (y == 1) || (fabs(xc_a - 1) > 0.25) || (fabs(yc_b - 1) > 0.25)) { // The above logic fails, the result is almost certainly zero: power1 = exp(log(xc_a) * a + log(yc_b) * b); power2 = 1; } else if (b > a) { T p = boost::math::expm1((b / a) * boost::math::log1p((y * a - x * b) / b)); power1 = exp(a * boost::math::log1p((x * b - y * a) / a + p * (x * c / a))); power2 = 1; } else { T p = boost::math::expm1((a / b) * boost::math::log1p((x * b - y * a) / a)); power1 = exp(b * boost::math::log1p((y * a - x * b) / b + p * (y * c / b))); power2 = 1; } } return prefix * power1 * power2 * scaled_tgamma_no_lanczos(c, pol) / (scaled_tgamma_no_lanczos(a, pol) * scaled_tgamma_no_lanczos(b, pol)); } T power1 = pow(x, a); T power2 = pow(y, b); T bet = beta_imp(a, b, l, pol); if(!(boost::math::isnormal)(power1) || !(boost::math::isnormal)(power2) || !(boost::math::isnormal)(bet)) { int shift_c = shift_a + shift_b; T result = ibeta_power_terms(T(a + shift_a), T(b + shift_b), x, y, l, normalised, pol, prefix); if ((boost::math::isnormal)(result)) { for (int i = 0; i < shift_c; ++i) { result /= c + i; if (i < shift_a) { result *= a + i; result /= x; } if (i < shift_b) { result *= b + i; result /= y; } } return prefix * result; } else { T log_result = log(x) * a + log(y) * b + log(prefix); if ((boost::math::isnormal)(bet)) log_result -= log(bet); else log_result += boost::math::lgamma(c, pol) - boost::math::lgamma(a, pol) - boost::math::lgamma(b, pol); return exp(log_result); } } return prefix * power1 * (power2 / bet); } #endif // // Series approximation to the incomplete beta: // template <class T> struct ibeta_series_t { typedef T result_type; BOOST_MATH_GPU_ENABLED ibeta_series_t(T a_, T b_, T x_, T mult) : result(mult), x(x_), apn(a_), poch(1-b_), n(1) {} BOOST_MATH_GPU_ENABLED T operator()() { T r = result / apn; apn += 1; result *= poch * x / n; ++n; poch += 1; return r; } private: T result, x, apn, poch; int n; }; template <class T, class Lanczos, class Policy> BOOST_MATH_GPU_ENABLED T ibeta_series(T a, T b, T x, T s0, const Lanczos&, bool normalised, T* p_derivative, T y, const Policy& pol) { BOOST_MATH_STD_USING T result; BOOST_MATH_ASSERT((p_derivative == 0) || normalised); if(normalised) { T c = a + b; // incomplete beta power term, combined with the Lanczos approximation: 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); 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; // LCOV_EXCL_LINE we can probably never get here, covered already above? T l1 = log(cgh / bgh) * (b - 0.5f); T l2 = log(x * cgh / agh) * a; // // Check for over/underflow in the power terms: // if((l1 > tools::log_min_value<T>()) && (l1 < tools::log_max_value<T>()) && (l2 > tools::log_min_value<T>()) && (l2 < tools::log_max_value<T>())) { if(a * b < bgh * 10) result *= exp((b - 0.5f) * boost::math::log1p(a / bgh, pol)); else result *= pow(cgh / bgh, T(b - T(0.5))); result *= pow(x * cgh / agh, a); result *= sqrt(agh / boost::math::constants::e<T>()); if(p_derivative) { *p_derivative = result * pow(y, b); BOOST_MATH_ASSERT(*p_derivative >= 0); } } else { // // Oh dear, we need logs, and this *will* cancel: // 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 { // Non-normalised, just compute the power: result = pow(x, a); } if(result < tools::min_value<T>()) return s0; // Safeguard: series can't cope with denorms. ibeta_series_t<T> s(a, b, x, result); boost::math::uintmax_t max_iter = policies::get_max_series_iterations<Policy>(); result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter, s0); policies::check_series_iterations<T>("boost::math::ibeta<%1%>(%1%, %1%, %1%) in ibeta_series (with lanczos)", max_iter, pol); return result; } // // Incomplete Beta series again, this time without Lanczos support: // #ifndef BOOST_MATH_HAS_GPU_SUPPORT template <class T, class Policy> BOOST_MATH_GPU_ENABLED T ibeta_series(T a, T b, T x, T s0, const boost::math::lanczos::undefined_lanczos& l, bool normalised, T* p_derivative, T y, const Policy& pol) { BOOST_MATH_STD_USING T result; BOOST_MATH_ASSERT((p_derivative == 0) || normalised); if(normalised) { const T min_sterling = minimum_argument_for_bernoulli_recursion<T>(); long shift_a = 0; long shift_b = 0; if (a < min_sterling) shift_a = 1 + ltrunc(min_sterling - a); if (b < min_sterling) shift_b = 1 + ltrunc(min_sterling - b); T c = a + b; if ((shift_a == 0) && (shift_b == 0)) { result = pow(x * c / a, a) * pow(c / b, b) * scaled_tgamma_no_lanczos(c, pol) / (scaled_tgamma_no_lanczos(a, pol) * scaled_tgamma_no_lanczos(b, pol)); } else if ((a < 1) && (b > 1)) result = pow(x, a) / (boost::math::tgamma(a, pol) * boost::math::tgamma_delta_ratio(b, a, pol)); else { T power = pow(x, a); T bet = beta_imp(a, b, l, pol); if (!(boost::math::isnormal)(power) || !(boost::math::isnormal)(bet)) { result = exp(a * log(x) + boost::math::lgamma(c, pol) - boost::math::lgamma(a, pol) - boost::math::lgamma(b, pol)); } else result = power / bet; } if(p_derivative) { *p_derivative = result * pow(y, b); BOOST_MATH_ASSERT(*p_derivative >= 0); } } else { // Non-normalised, just compute the power: result = pow(x, a); } if(result < tools::min_value<T>()) return s0; // Safeguard: series can't cope with denorms. ibeta_series_t<T> s(a, b, x, result); boost::math::uintmax_t max_iter = policies::get_max_series_iterations<Policy>(); result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter, s0); policies::check_series_iterations<T>("boost::math::ibeta<%1%>(%1%, %1%, %1%) in ibeta_series (without lanczos)", max_iter, pol); return result; } #endif // // Continued fraction for the incomplete beta: // template <class T> 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 result_type operator()() { T aN = (a + m - 1) * (a + b + m - 1) * m * (b - m) * x * x; T denom = (a + 2 * m - 1); aN /= denom * denom; T bN = static_cast<T>(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; return boost::math::make_pair(aN, bN); } private: T a, b, x, y; int m; }; // // Evaluate the incomplete beta via the continued fraction representation: // template <class T, class Policy> BOOST_MATH_GPU_ENABLED inline T ibeta_fraction2(T a, T b, T x, T y, const Policy& pol, bool normalised, T* p_derivative) { typedef typename lanczos::lanczos<T, Policy>::type lanczos_type; BOOST_MATH_STD_USING T result = ibeta_power_terms(a, b, x, y, lanczos_type(), normalised, pol); if(p_derivative) { *p_derivative = result; BOOST_MATH_ASSERT(*p_derivative >= 0); } if(result == 0) return result; ibeta_fraction2_t<T> f(a, b, x, y); T fract = boost::math::tools::continued_fraction_b(f, boost::math::policies::get_epsilon<T, Policy>()); BOOST_MATH_INSTRUMENT_VARIABLE(fract); BOOST_MATH_INSTRUMENT_VARIABLE(result); return result / fract; } // // Computes the difference between ibeta(a,b,x) and ibeta(a+k,b,x): // template <class T, class Policy> BOOST_MATH_GPU_ENABLED T ibeta_a_step(T a, T b, T x, T y, int k, const Policy& pol, bool normalised, T* p_derivative) { typedef typename lanczos::lanczos<T, Policy>::type lanczos_type; BOOST_MATH_INSTRUMENT_VARIABLE(k); T prefix = ibeta_power_terms(a, b, x, y, lanczos_type(), normalised, pol); if(p_derivative) { *p_derivative = prefix; BOOST_MATH_ASSERT(*p_derivative >= 0); } prefix /= a; if(prefix == 0) return prefix; T sum = 1; T term = 1; // series summation from 0 to k-1: for(int i = 0; i < k-1; ++i) { term *= (a+b+i) * x / (a+i+1); sum += term; } prefix *= sum; return prefix; } // // This function is only needed for the non-regular incomplete beta, // it computes the delta in: // beta(a,b,x) = prefix + delta * beta(a+k,b,x) // it is currently only called for small k. // template <class T> BOOST_MATH_GPU_ENABLED inline T rising_factorial_ratio(T a, T b, int k) { // calculate: // (a)(a+1)(a+2)...(a+k-1) // _______________________ // (b)(b+1)(b+2)...(b+k-1) // This is only called with small k, for large k // it is grossly inefficient, do not use outside it's // intended purpose!!! BOOST_MATH_INSTRUMENT_VARIABLE(k); BOOST_MATH_ASSERT(k > 0); T result = 1; for(int i = 0; i < k; ++i) result *= (a+i) / (b+i); return result; } // // Routine for a > 15, b < 1 // // Begin by figuring out how large our table of Pn's should be, // quoted accuracies are "guesstimates" based on empirical observation. // Note that the table size should never exceed the size of our // tables of factorials. // template <class T> struct Pn_size { // This is likely to be enough for ~35-50 digit accuracy // but it's hard to quantify exactly: #ifndef BOOST_MATH_HAS_NVRTC static constexpr unsigned value = ::boost::math::max_factorial<T>::value >= 100 ? 50 : ::boost::math::max_factorial<T>::value >= ::boost::math::max_factorial<double>::value ? 30 : ::boost::math::max_factorial<T>::value >= ::boost::math::max_factorial<float>::value ? 15 : 1; static_assert(::boost::math::max_factorial<T>::value >= ::boost::math::max_factorial<float>::value, "Type does not provide for 35-50 digits of accuracy."); #else static constexpr unsigned value = 0; // Will never be called #endif }; template <> struct Pn_size<float> { static constexpr unsigned value = 15; // ~8-15 digit accuracy #ifndef BOOST_MATH_HAS_GPU_SUPPORT static_assert(::boost::math::max_factorial<float>::value >= 30, "Type does not provide for 8-15 digits of accuracy."); #endif }; template <> struct Pn_size<double> { static constexpr unsigned value = 30; // 16-20 digit accuracy #ifndef BOOST_MATH_HAS_GPU_SUPPORT static_assert(::boost::math::max_factorial<double>::value >= 60, "Type does not provide for 16-20 digits of accuracy."); #endif }; template <> struct Pn_size<long double> { static constexpr unsigned value = 50; // ~35-50 digit accuracy #ifndef BOOST_MATH_HAS_GPU_SUPPORT static_assert(::boost::math::max_factorial<long double>::value >= 100, "Type does not provide for ~35-50 digits of accuracy"); #endif }; template <class T, class Policy> BOOST_MATH_GPU_ENABLED T beta_small_b_large_a_series(T a, T b, T x, T y, T s0, T mult, const Policy& pol, bool normalised) { typedef typename lanczos::lanczos<T, Policy>::type lanczos_type; BOOST_MATH_STD_USING // // This is DiDonato and Morris's BGRAT routine, see Eq's 9 through 9.6. // // Some values we'll need later, these are Eq 9.1: // T bm1 = b - 1; T t = a + bm1 / 2; T lx, u; // LCOV_EXCL_LINE if(y < 0.35) lx = boost::math::log1p(-y, pol); else lx = log(x); u = -t * lx; // and from from 9.2: T prefix; // LCOV_EXCL_LINE T h = regularised_gamma_prefix(b, u, pol, lanczos_type()); if(h <= tools::min_value<T>()) return s0; if(normalised) { prefix = h / boost::math::tgamma_delta_ratio(a, b, pol); prefix /= pow(t, b); } else { prefix = full_igamma_prefix(b, u, pol) / pow(t, b); } prefix *= mult; // // now we need the quantity Pn, unfortunately this is computed // recursively, and requires a full history of all the previous values // so no choice but to declare a big table and hope it's big enough... // T p[ ::boost::math::detail::Pn_size<T>::value ] = { 1 }; // see 9.3. // // Now an initial value for J, see 9.6: // T j = boost::math::gamma_q(b, u, pol) / h; // // Now we can start to pull things together and evaluate the sum in Eq 9: // T sum = s0 + prefix * j; // Value at N = 0 // some variables we'll need: unsigned tnp1 = 1; // 2*N+1 T lx2 = lx / 2; lx2 *= lx2; T lxp = 1; T t4 = 4 * t * t; T b2n = b; for(unsigned n = 1; n < sizeof(p)/sizeof(p[0]); ++n) { /* // debugging code, enable this if you want to determine whether // the table of Pn's is large enough... // static int max_count = 2; if(n > max_count) { max_count = n; std::cerr << "Max iterations in BGRAT was " << n << std::endl; } */ // // begin by evaluating the next Pn from Eq 9.4: // tnp1 += 2; p[n] = 0; T mbn = b - n; unsigned tmp1 = 3; for(unsigned m = 1; m < n; ++m) { mbn = m * b - n; p[n] += mbn * p[n-m] / boost::math::unchecked_factorial<T>(tmp1); tmp1 += 2; } p[n] /= n; p[n] += bm1 / boost::math::unchecked_factorial<T>(tnp1); // // Now we want Jn from Jn-1 using Eq 9.6: // j = (b2n * (b2n + 1) * j + (u + b2n + 1) * lxp) / t4; lxp *= lx2; b2n += 2; // // pull it together with Eq 9: // T r = prefix * p[n] * j; sum += r; // r is always small: BOOST_MATH_ASSERT(tools::max_value<T>() * tools::epsilon<T>() > fabs(r)); if(fabs(r / tools::epsilon<T>()) < fabs(sum)) break; } return sum; } // template <class T, class Lanczos>T beta_small_b_large_a_series(T a, T b, T x, T y, T s0, T mult, const Lanczos& l, bool normalised) // // For integer arguments we can relate the incomplete beta to the // complement of the binomial distribution cdf and use this finite sum. // template <class T, class Policy> BOOST_MATH_GPU_ENABLED T binomial_ccdf(T n, T k, T x, T y, const Policy& pol) { BOOST_MATH_STD_USING // ADL of std names T result = pow(x, n); if(result > tools::min_value<T>()) { T term = result; for(unsigned i = itrunc(T(n - 1)); i > k; --i) { term *= ((i + 1) * y) / ((n - i) * x); result += term; } } else { // First term underflows so we need to start at the mode of the // distribution and work outwards: int start = itrunc(n * x); if(start <= k + 1) start = itrunc(k + 2); result = static_cast<T>(pow(x, T(start)) * pow(y, n - T(start)) * boost::math::binomial_coefficient<T>(itrunc(n), itrunc(start), pol)); if(result == 0) { // OK, starting slightly above the mode didn't work, // we'll have to sum the terms the old fashioned way. // Very hard to get here, possibly only when exponent // range is very limited (as with type float): // LCOV_EXCL_START for(unsigned i = start - 1; i > k; --i) { result += static_cast<T>(pow(x, static_cast<T>(i)) * pow(y, n - i) * boost::math::binomial_coefficient<T>(itrunc(n), itrunc(i), pol)); } // LCOV_EXCL_STOP } else { T term = result; T start_term = result; for(unsigned i = start - 1; i > k; --i) { term *= ((i + 1) * y) / ((n - i) * x); result += term; } term = start_term; for(unsigned i = start + 1; i <= n; ++i) { term *= (n - i + 1) * x / (i * y); result += term; } } } return result; } // // The incomplete beta function implementation: // This is just a big bunch of spaghetti code to divide up the // input range and select the right implementation method for // each domain: // template <class T, class Policy> BOOST_MATH_GPU_ENABLED T ibeta_imp(T a, T b, T x, const Policy& pol, bool inv, bool normalised, T* p_derivative) { constexpr auto function = "boost::math::ibeta<%1%>(%1%, %1%, %1%)"; typedef typename lanczos::lanczos<T, Policy>::type lanczos_type; BOOST_MATH_STD_USING // for ADL of std math functions. BOOST_MATH_INSTRUMENT_VARIABLE(a); BOOST_MATH_INSTRUMENT_VARIABLE(b); BOOST_MATH_INSTRUMENT_VARIABLE(x); BOOST_MATH_INSTRUMENT_VARIABLE(inv); BOOST_MATH_INSTRUMENT_VARIABLE(normalised); bool invert = inv; T fract; T y = 1 - x; BOOST_MATH_ASSERT((p_derivative == 0) || normalised); if(!(boost::math::isfinite)(a)) return policies::raise_domain_error<T>(function, "The argument a to the incomplete beta function must be finite (got a=%1%).", a, pol); if(!(boost::math::isfinite)(b)) return policies::raise_domain_error<T>(function, "The argument b to the incomplete beta function must be finite (got b=%1%).", b, pol); if (!(0 <= x && x <= 1)) return policies::raise_domain_error<T>(function, "The argument x to the incomplete beta function must be in [0,1] (got x=%1%).", x, pol); if(p_derivative) *p_derivative = -1; // value not set. if(normalised) { if(a < 0) return policies::raise_domain_error<T>(function, "The argument a to the incomplete beta function must be >= zero (got a=%1%).", a, pol); if(b < 0) return policies::raise_domain_error<T>(function, "The argument b to the incomplete beta function must be >= zero (got b=%1%).", b, pol); // extend to a few very special cases: if(a == 0) { if(b == 0) return policies::raise_domain_error<T>(function, "The arguments a and b to the incomplete beta function cannot both be zero, with x=%1%.", x, pol); if(b > 0) return static_cast<T>(inv ? 0 : 1); } else if(b == 0) { if(a > 0) return static_cast<T>(inv ? 1 : 0); } } else { if(a <= 0) return policies::raise_domain_error<T>(function, "The argument a to the incomplete beta function must be greater than zero (got a=%1%).", a, pol); if(b <= 0) return policies::raise_domain_error<T>(function, "The argument b to the incomplete beta function must be greater than zero (got b=%1%).", b, pol); } if(x == 0) { if(p_derivative) { *p_derivative = (a == 1) ? (T)1 : (a < 1) ? T(tools::max_value<T>() / 2) : T(tools::min_value<T>() * 2); } return (invert ? (normalised ? T(1) : boost::math::beta(a, b, pol)) : T(0)); } if(x == 1) { if(p_derivative) { *p_derivative = (b == 1) ? T(1) : (b < 1) ? T(tools::max_value<T>() / 2) : T(tools::min_value<T>() * 2); } return (invert == 0 ? (normalised ? 1 : boost::math::beta(a, b, pol)) : 0); } if((a == 0.5f) && (b == 0.5f)) { // We have an arcsine distribution: if(p_derivative) { *p_derivative = 1 / (constants::pi<T>() * sqrt(y * x)); } T p = invert ? asin(sqrt(y)) / constants::half_pi<T>() : asin(sqrt(x)) / constants::half_pi<T>(); if(!normalised) p *= constants::pi<T>(); return p; } if(a == 1) { BOOST_MATH_GPU_SAFE_SWAP(a, b); BOOST_MATH_GPU_SAFE_SWAP(x, y); invert = !invert; } if(b == 1) { // // Special case see: http://functions.wolfram.com/GammaBetaErf/BetaRegularized/03/01/01/ // if(a == 1) { if(p_derivative) *p_derivative = 1; return invert ? y : x; } if(p_derivative) { *p_derivative = a * pow(x, a - 1); } T p; // LCOV_EXCL_LINE if(y < 0.5) p = invert ? T(-boost::math::expm1(a * boost::math::log1p(-y, pol), pol)) : T(exp(a * boost::math::log1p(-y, pol))); else p = invert ? T(-boost::math::powm1(x, a, pol)) : T(pow(x, a)); if(!normalised) p /= a; return p; } if(BOOST_MATH_GPU_SAFE_MIN(a, b) <= 1) { if(x > 0.5) { BOOST_MATH_GPU_SAFE_SWAP(a, b); BOOST_MATH_GPU_SAFE_SWAP(x, y); invert = !invert; BOOST_MATH_INSTRUMENT_VARIABLE(invert); } if(BOOST_MATH_GPU_SAFE_MAX(a, b) <= 1) { // Both a,b < 1: if((a >= BOOST_MATH_GPU_SAFE_MIN(T(0.2), b)) || (pow(x, a) <= 0.9)) { if(!invert) { fract = ibeta_series(a, b, x, T(0), lanczos_type(), normalised, p_derivative, y, pol); BOOST_MATH_INSTRUMENT_VARIABLE(fract); } else { fract = -(normalised ? 1 : boost::math::beta(a, b, pol)); invert = false; fract = -ibeta_series(a, b, x, fract, lanczos_type(), normalised, p_derivative, y, pol); BOOST_MATH_INSTRUMENT_VARIABLE(fract); } } else { BOOST_MATH_GPU_SAFE_SWAP(a, b); BOOST_MATH_GPU_SAFE_SWAP(x, y); invert = !invert; if(y >= 0.3) { if(!invert) { fract = ibeta_series(a, b, x, T(0), lanczos_type(), normalised, p_derivative, y, pol); BOOST_MATH_INSTRUMENT_VARIABLE(fract); } else { fract = -(normalised ? 1 : boost::math::beta(a, b, pol)); invert = false; fract = -ibeta_series(a, b, x, fract, lanczos_type(), normalised, p_derivative, y, pol); BOOST_MATH_INSTRUMENT_VARIABLE(fract); } } else { // Sidestep on a, and then use the series representation: T prefix; // LCOV_EXCL_LINE if(!normalised) { prefix = rising_factorial_ratio(T(a+b), a, 20); } else { prefix = 1; } fract = ibeta_a_step(a, b, x, y, 20, pol, normalised, p_derivative); if(!invert) { fract = beta_small_b_large_a_series(T(a + 20), b, x, y, fract, prefix, pol, normalised); BOOST_MATH_INSTRUMENT_VARIABLE(fract); } else { fract -= (normalised ? 1 : boost::math::beta(a, b, pol)); invert = false; fract = -beta_small_b_large_a_series(T(a + 20), b, x, y, fract, prefix, pol, normalised); BOOST_MATH_INSTRUMENT_VARIABLE(fract); } } } } else { // One of a, b < 1 only: if((b <= 1) || ((x < 0.1) && (pow(b * x, a) <= 0.7))) { if(!invert) { fract = ibeta_series(a, b, x, T(0), lanczos_type(), normalised, p_derivative, y, pol); BOOST_MATH_INSTRUMENT_VARIABLE(fract); } else { fract = -(normalised ? 1 : boost::math::beta(a, b, pol)); invert = false; fract = -ibeta_series(a, b, x, fract, lanczos_type(), normalised, p_derivative, y, pol); BOOST_MATH_INSTRUMENT_VARIABLE(fract); } } else { BOOST_MATH_GPU_SAFE_SWAP(a, b); BOOST_MATH_GPU_SAFE_SWAP(x, y); invert = !invert; if(y >= 0.3) { if(!invert) { fract = ibeta_series(a, b, x, T(0), lanczos_type(), normalised, p_derivative, y, pol); BOOST_MATH_INSTRUMENT_VARIABLE(fract); } else { fract = -(normalised ? 1 : boost::math::beta(a, b, pol)); invert = false; fract = -ibeta_series(a, b, x, fract, lanczos_type(), normalised, p_derivative, y, pol); BOOST_MATH_INSTRUMENT_VARIABLE(fract); } } else if(a >= 15) { if(!invert) { fract = beta_small_b_large_a_series(a, b, x, y, T(0), T(1), pol, normalised); BOOST_MATH_INSTRUMENT_VARIABLE(fract); } else { fract = -(normalised ? 1 : boost::math::beta(a, b, pol)); invert = false; fract = -beta_small_b_large_a_series(a, b, x, y, fract, T(1), pol, normalised); BOOST_MATH_INSTRUMENT_VARIABLE(fract); } } else { // Sidestep to improve errors: T prefix; // LCOV_EXCL_LINE if(!normalised) { prefix = rising_factorial_ratio(T(a+b), a, 20); } else { prefix = 1; } fract = ibeta_a_step(a, b, x, y, 20, pol, normalised, p_derivative); BOOST_MATH_INSTRUMENT_VARIABLE(fract); if(!invert) { fract = beta_small_b_large_a_series(T(a + 20), b, x, y, fract, prefix, pol, normalised); BOOST_MATH_INSTRUMENT_VARIABLE(fract); } else { fract -= (normalised ? 1 : boost::math::beta(a, b, pol)); invert = false; fract = -beta_small_b_large_a_series(T(a + 20), b, x, y, fract, prefix, pol, normalised); BOOST_MATH_INSTRUMENT_VARIABLE(fract); } } } } } else { // Both a,b >= 1: T lambda; // LCOV_EXCL_LINE if(a < b) { lambda = a - (a + b) * x; } else { lambda = (a + b) * y - b; } if(lambda < 0) { BOOST_MATH_GPU_SAFE_SWAP(a, b); BOOST_MATH_GPU_SAFE_SWAP(x, y); invert = !invert; BOOST_MATH_INSTRUMENT_VARIABLE(invert); } if(b < 40) { if((floor(a) == a) && (floor(b) == b) && (a < static_cast<T>((boost::math::numeric_limits<int>::max)() - 100)) && (y != 1)) { // relate to the binomial distribution and use a finite sum: T k = a - 1; T n = b + k; fract = binomial_ccdf(n, k, x, y, pol); if(!normalised) fract *= boost::math::beta(a, b, pol); BOOST_MATH_INSTRUMENT_VARIABLE(fract); } else if(b * x <= 0.7) { if(!invert) { fract = ibeta_series(a, b, x, T(0), lanczos_type(), normalised, p_derivative, y, pol); BOOST_MATH_INSTRUMENT_VARIABLE(fract); } else { fract = -(normalised ? 1 : boost::math::beta(a, b, pol)); invert = false; fract = -ibeta_series(a, b, x, fract, lanczos_type(), normalised, p_derivative, y, pol); BOOST_MATH_INSTRUMENT_VARIABLE(fract); } } else if(a > 15) { // sidestep so we can use the series representation: int n = itrunc(T(floor(b)), pol); if(n == b) --n; T bbar = b - n; T prefix; // LCOV_EXCL_LINE if(!normalised) { prefix = rising_factorial_ratio(T(a+bbar), bbar, n); } else { prefix = 1; } fract = ibeta_a_step(bbar, a, y, x, n, pol, normalised, static_cast<T*>(nullptr)); fract = beta_small_b_large_a_series(a, bbar, x, y, fract, T(1), pol, normalised); fract /= prefix; BOOST_MATH_INSTRUMENT_VARIABLE(fract); } else if(normalised) { // The formula here for the non-normalised case is tricky to figure // out (for me!!), and requires two pochhammer calculations rather // than one, so leave it for now and only use this in the normalized case.... int n = itrunc(T(floor(b)), pol); T bbar = b - n; if(bbar <= 0) { --n; bbar += 1; } fract = ibeta_a_step(bbar, a, y, x, n, pol, normalised, static_cast<T*>(nullptr)); fract += ibeta_a_step(a, bbar, x, y, 20, pol, normalised, static_cast<T*>(nullptr)); if(invert) fract -= 1; // Note this line would need changing if we ever enable this branch in non-normalized case fract = beta_small_b_large_a_series(T(a+20), bbar, x, y, fract, T(1), pol, normalised); if(invert) { fract = -fract; invert = false; } BOOST_MATH_INSTRUMENT_VARIABLE(fract); } else { fract = ibeta_fraction2(a, b, x, y, pol, normalised, p_derivative); BOOST_MATH_INSTRUMENT_VARIABLE(fract); } } else { fract = ibeta_fraction2(a, b, x, y, pol, normalised, p_derivative); BOOST_MATH_INSTRUMENT_VARIABLE(fract); } } if(p_derivative) { if(*p_derivative < 0) { *p_derivative = ibeta_power_terms(a, b, x, y, lanczos_type(), true, pol); } T div = y * x; if(*p_derivative != 0) { if((tools::max_value<T>() * div < *p_derivative)) { // overflow, return an arbitrarily large value: *p_derivative = tools::max_value<T>() / 2; // LCOV_EXCL_LINE Probably can only get here with denormalized x. } else { *p_derivative /= div; } } } return invert ? (normalised ? 1 : boost::math::beta(a, b, pol)) - fract : fract; } // template <class T, class Lanczos>T ibeta_imp(T a, T b, T x, const Lanczos& l, bool inv, bool normalised) template <class T, class Policy> BOOST_MATH_GPU_ENABLED inline T ibeta_imp(T a, T b, T x, const Policy& pol, bool inv, bool normalised) { return ibeta_imp(a, b, x, pol, inv, normalised, static_cast<T*>(nullptr)); } template <class T, class Policy> BOOST_MATH_GPU_ENABLED T ibeta_derivative_imp(T a, T b, T x, const Policy& pol) { constexpr auto function = "ibeta_derivative<%1%>(%1%,%1%,%1%)"; // // start with the usual error checks: // if (!(boost::math::isfinite)(a)) return policies::raise_domain_error<T>(function, "The argument a to the incomplete beta function must be finite (got a=%1%).", a, pol); if (!(boost::math::isfinite)(b)) return policies::raise_domain_error<T>(function, "The argument b to the incomplete beta function must be finite (got b=%1%).", b, pol); if (!(0 <= x && x <= 1)) return policies::raise_domain_error<T>(function, "The argument x to the incomplete beta function must be in [0,1] (got x=%1%).", x, pol); if(a <= 0) return policies::raise_domain_error<T>(function, "The argument a to the incomplete beta function must be greater than zero (got a=%1%).", a, pol); if(b <= 0) return policies::raise_domain_error<T>(function, "The argument b to the incomplete beta function must be greater than zero (got b=%1%).", b, pol); // // Now the corner cases: // if(x == 0) { return (a > 1) ? 0 : (a == 1) ? 1 / boost::math::beta(a, b, pol) : policies::raise_overflow_error<T>(function, nullptr, pol); } else if(x == 1) { return (b > 1) ? 0 : (b == 1) ? 1 / boost::math::beta(a, b, pol) : policies::raise_overflow_error<T>(function, nullptr, pol); } // // Now the regular cases: // typedef typename lanczos::lanczos<T, Policy>::type lanczos_type; T y = (1 - x) * x; T f1; if (!(boost::math::isinf)(1 / y)) { f1 = ibeta_power_terms<T>(a, b, x, 1 - x, lanczos_type(), true, pol, 1 / y, function); } else { return (a > 1) ? 0 : (a == 1) ? 1 / boost::math::beta(a, b, pol) : policies::raise_overflow_error<T>(function, nullptr, pol); } return f1; } // // Some forwarding functions that disambiguate the third argument type: // template <class RT1, class RT2, class Policy> BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<RT1, RT2>::type beta(RT1 a, RT2 b, const Policy&, const boost::math::true_type*) { BOOST_FPU_EXCEPTION_GUARD typedef typename tools::promote_args<RT1, RT2>::type result_type; typedef typename policies::evaluation<result_type, Policy>::type value_type; typedef typename lanczos::lanczos<value_type, Policy>::type evaluation_type; typedef typename policies::normalise< Policy, policies::promote_float<false>, policies::promote_double<false>, policies::discrete_quantile<>, policies::assert_undefined<> >::type forwarding_policy; return policies::checked_narrowing_cast<result_type, forwarding_policy>(detail::beta_imp(static_cast<value_type>(a), static_cast<value_type>(b), evaluation_type(), forwarding_policy()), "boost::math::beta<%1%>(%1%,%1%)"); } template <class RT1, class RT2, class RT3> BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<RT1, RT2, RT3>::type beta(RT1 a, RT2 b, RT3 x, const boost::math::false_type*) { return boost::math::beta(a, b, x, policies::policy<>()); } } // namespace detail // // The actual function entry-points now follow, these just figure out // which Lanczos approximation to use // and forward to the implementation functions: // template <class RT1, class RT2, class A> BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<RT1, RT2, A>::type beta(RT1 a, RT2 b, A arg) { using tag = typename policies::is_policy<A>::type; using ReturnType = tools::promote_args_t<RT1, RT2, A>; return static_cast<ReturnType>(boost::math::detail::beta(a, b, arg, static_cast<tag*>(nullptr))); } template <class RT1, class RT2> BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<RT1, RT2>::type beta(RT1 a, RT2 b) { return boost::math::beta(a, b, policies::policy<>()); } template <class RT1, class RT2, class RT3, class Policy> BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<RT1, RT2, RT3>::type beta(RT1 a, RT2 b, RT3 x, const Policy&) { BOOST_FPU_EXCEPTION_GUARD typedef typename tools::promote_args<RT1, RT2, RT3>::type result_type; typedef typename policies::evaluation<result_type, Policy>::type value_type; typedef typename policies::normalise< Policy, policies::promote_float<false>, policies::promote_double<false>, policies::discrete_quantile<>, policies::assert_undefined<> >::type forwarding_policy; return policies::checked_narrowing_cast<result_type, forwarding_policy>(detail::ibeta_imp(static_cast<value_type>(a), static_cast<value_type>(b), static_cast<value_type>(x), forwarding_policy(), false, false), "boost::math::beta<%1%>(%1%,%1%,%1%)"); } template <class RT1, class RT2, class RT3, class Policy> BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<RT1, RT2, RT3>::type betac(RT1 a, RT2 b, RT3 x, const Policy&) { BOOST_FPU_EXCEPTION_GUARD typedef typename tools::promote_args<RT1, RT2, RT3>::type result_type; typedef typename policies::evaluation<result_type, Policy>::type value_type; typedef typename policies::normalise< Policy, policies::promote_float<false>, policies::promote_double<false>, policies::discrete_quantile<>, policies::assert_undefined<> >::type forwarding_policy; return policies::checked_narrowing_cast<result_type, forwarding_policy>(detail::ibeta_imp(static_cast<value_type>(a), static_cast<value_type>(b), static_cast<value_type>(x), forwarding_policy(), true, false), "boost::math::betac<%1%>(%1%,%1%,%1%)"); } template <class RT1, class RT2, class RT3> BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<RT1, RT2, RT3>::type betac(RT1 a, RT2 b, RT3 x) { return boost::math::betac(a, b, x, policies::policy<>()); } template <class RT1, class RT2, class RT3, class Policy> BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<RT1, RT2, RT3>::type ibeta(RT1 a, RT2 b, RT3 x, const Policy&) { BOOST_FPU_EXCEPTION_GUARD typedef typename tools::promote_args<RT1, RT2, RT3>::type result_type; typedef typename policies::evaluation<result_type, Policy>::type value_type; typedef typename policies::normalise< Policy, policies::promote_float<false>, policies::promote_double<false>, policies::discrete_quantile<>, policies::assert_undefined<> >::type forwarding_policy; return policies::checked_narrowing_cast<result_type, forwarding_policy>(detail::ibeta_imp(static_cast<value_type>(a), static_cast<value_type>(b), static_cast<value_type>(x), forwarding_policy(), false, true), "boost::math::ibeta<%1%>(%1%,%1%,%1%)"); } template <class RT1, class RT2, class RT3> BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<RT1, RT2, RT3>::type ibeta(RT1 a, RT2 b, RT3 x) { return boost::math::ibeta(a, b, x, policies::policy<>()); } template <class RT1, class RT2, class RT3, class Policy> BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<RT1, RT2, RT3>::type ibetac(RT1 a, RT2 b, RT3 x, const Policy&) { BOOST_FPU_EXCEPTION_GUARD typedef typename tools::promote_args<RT1, RT2, RT3>::type result_type; typedef typename policies::evaluation<result_type, Policy>::type value_type; typedef typename policies::normalise< Policy, policies::promote_float<false>, policies::promote_double<false>, policies::discrete_quantile<>, policies::assert_undefined<> >::type forwarding_policy; return policies::checked_narrowing_cast<result_type, forwarding_policy>(detail::ibeta_imp(static_cast<value_type>(a), static_cast<value_type>(b), static_cast<value_type>(x), forwarding_policy(), true, true), "boost::math::ibetac<%1%>(%1%,%1%,%1%)"); } template <class RT1, class RT2, class RT3> BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<RT1, RT2, RT3>::type ibetac(RT1 a, RT2 b, RT3 x) { return boost::math::ibetac(a, b, x, policies::policy<>()); } template <class RT1, class RT2, class RT3, class Policy> BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<RT1, RT2, RT3>::type ibeta_derivative(RT1 a, RT2 b, RT3 x, const Policy&) { BOOST_FPU_EXCEPTION_GUARD typedef typename tools::promote_args<RT1, RT2, RT3>::type result_type; typedef typename policies::evaluation<result_type, Policy>::type value_type; typedef typename policies::normalise< Policy, policies::promote_float<false>, policies::promote_double<false>, policies::discrete_quantile<>, policies::assert_undefined<> >::type forwarding_policy; return policies::checked_narrowing_cast<result_type, forwarding_policy>(detail::ibeta_derivative_imp(static_cast<value_type>(a), static_cast<value_type>(b), static_cast<value_type>(x), forwarding_policy()), "boost::math::ibeta_derivative<%1%>(%1%,%1%,%1%)"); } template <class RT1, class RT2, class RT3> BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<RT1, RT2, RT3>::type ibeta_derivative(RT1 a, RT2 b, RT3 x) { return boost::math::ibeta_derivative(a, b, x, policies::policy<>()); } } // namespace math } // namespace boost #include <boost/math/special_functions/detail/ibeta_inverse.hpp> #include <boost/math/special_functions/detail/ibeta_inv_ab.hpp> #endif // BOOST_MATH_SPECIAL_BETA_HPP