shylie/llvm-project
1
0
Fork 0
Code Issues Pull Requests Actions 6 Packages Projects Releases Wiki Activity
llvm-project/third-party/boost-math/include/boost/math/special_functions/beta.hpp
Louis Dionne 585da50d7d
[third-party] Add a snapshot of Boost.Math 1.89 standalone (#141508) 
This PR adds the code of Boost.Math as of version 1.89 into the
third-party directory, as discussed in a recent RFC [1].

The goal is for this code to be used as a back-end for the C++17
Math Special Functions.

As explained in third-paty/README.md, this code is cleared for
usage inside libc++ for the Math Special functions, however
the LLVM Foundation should be consulted before using this
code anywhere else in the LLVM project, due to the fact
that it is under the Boost Software License (as opposed
to the usual LLVM license). See the RFC [1] for more details.

[1]: https://discourse.llvm.org/t/rfc-libc-taking-a-dependency-on-boost-math-for-the-c-17-math-special-functions
2025-10-27 14:43:57 -07:00

1850 lines
63 KiB
C++
Raw Blame History

// (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/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>
#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 gh = Lanczos::g() - 0.5f;
T agh = static_cast<T>(a + gh);
T bgh = static_cast<T>(b + gh);
T cgh = static_cast<T>(c + gh);
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));
BOOST_MATH_INSTRUMENT_VARIABLE(result);
result *= prefix;
BOOST_MATH_INSTRUMENT_VARIABLE(result);
// combine with the leftover terms from the Lanczos approximation:
result *= sqrt(bgh / boost::math::constants::e<T>());
result *= sqrt(agh / cgh);
BOOST_MATH_INSTRUMENT_VARIABLE(result);
// l1 and l2 are the base of the exponents minus one:
T l1 = ((x * b - y * a) - y * gh) / agh;
T l2 = ((y * a - x * b) - x * gh) / 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
{
T l1 = Lanczos::lanczos_sum_expG_scaled(c);
T l2 = Lanczos::lanczos_sum_expG_scaled(a);
T l3 = Lanczos::lanczos_sum_expG_scaled(b);
if ((l2 > 1) && (l3 > 1) && (tools::max_value<T>() / l2 < l3))
result = (l1 / l2) / l3;
else
result = l1 / (l2 * l3);
}
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 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);
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);
boost::math::uintmax_t max_terms = boost::math::policies::get_max_series_iterations<Policy>();
T fract = boost::math::tools::continued_fraction_b(f, boost::math::policies::get_epsilon<T, Policy>(), max_terms);
boost::math::policies::check_series_iterations<T>("boost::math::ibeta", max_terms, pol);
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;
}
template <class T, class Policy>
BOOST_MATH_GPU_ENABLED T ibeta_large_ab(T a, T b, T x, T y, bool invert, bool normalised, const Policy& pol)
{
//
// Large arguments, symetric case, see https://dlmf.nist.gov/8.18
//
BOOST_MATH_STD_USING
T x0 = a / (a + b);
T y0 = b / (a + b);
T nu = x0 * log(x / x0) + y0 * log(y / y0);
//
// Above compution is unstable, force nu to zero if
// something went wrong:
//
if ((nu > 0) || (x == x0) || (y == y0))
nu = 0;
nu = sqrt(-2 * nu);
//
// As per https://dlmf.nist.gov/8.18#E10 we need to make sure we have the correct root:
//
if ((nu != 0) && (nu / (x - x0) < 0))
nu = -nu;
//
// The correction term in https://dlmf.nist.gov/8.18#E9 is badly unstable, and often
// makes the compution worse not better, we exclude it for now:
/*
T c0 = 0;
if (nu != 0)
{
c0 = 1 / nu;
T lim = fabs(10 * tools::epsilon<T>() * c0);
c0 -= sqrt(x0 * y0) / (x - x0);
if(fabs(c0) < lim)
c0 = (1 - 2 * x0) / (3 * sqrt(x0 * y0));
else
c0 *= exp(a * log(x / x0) + b * log(y / y0));
c0 /= sqrt(constants::two_pi<T>() * (a + b));
}
else
{
c0 = (1 - 2 * x0) / (3 * sqrt(x0 * y0));
c0 /= sqrt(constants::two_pi<T>() * (a + b));
}
*/
T mul = 1;
if (!normalised)
mul = boost::math::beta(a, b, pol);
return mul * ((invert ? (1 + boost::math::erf(-nu * sqrt((a + b) / 2), pol)) / 2 : boost::math::erfc(-nu * sqrt((a + b) / 2), pol) / 2));
}
//
// 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
{
// 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 (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;
}
}
if(use_asym)
{
fract = ibeta_large_ab(a, b, x, y, invert, normalised, pol);
if (fract * tools::epsilon<T>() < powers)
{
// Erf approximation failed, correction term is too large, fall back:
fract = ibeta_fraction2(a, b, x, y, pol, normalised, p_derivative);
}
else
invert = false;
}
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
Reference in New Issue View Git Blame Copy Permalink