Synopsis

#include <boost/math/special_functions/beta.hpp>

namespace boost{ namespace math{

template <class T1, class T2>
BOOST_MATH_GPU_ENABLED calculated-result-type beta(T1 a, T2 b);

template <class T1, class T2, class Policy>
BOOST_MATH_GPU_ENABLED calculated-result-type beta(T1 a, T2 b, const Policy&);

template <class T1, class T2>
BOOST_MATH_GPU_ENABLED calculated-result-type lbeta(T1 a, T2 b);

template <class T1, class T2, class Policy>
BOOST_MATH_GPU_ENABLED calculated-result-type lbeta(T1 a, T2 b, const Policy&);

}} // namespaces

Description

The beta function is defined by:

The final Policy argument is optional and can be used to control the behaviour of the function: how it handles errors, what level of precision to use etc. Refer to the policy documentation for more details.

The return type of these functions is computed using the result type calculation rules when T1 and T2 are different types.

lbeta returns the natural logarithm of the beta function, log B(a, b). It stays finite and accurate where beta underflows, for example when both a and b are large, and where lgamma(a) + lgamma(b) - lgamma(a + b) would lose every digit to cancellation, when one argument is much larger than the other.

Accuracy

The following table shows peak errors for various domains of input arguments, along with comparisons to the GSL-1.9 and Cephes libraries. Note that only results for the widest floating point type on the system are given as narrower types have effectively zero error.

Note that the worst errors occur when a or b are large, and that when this is the case the result is very close to zero, so absolute errors will be very small.

Testing

A mixture of spot tests of exact values, and randomly generated test data are used: the test data was computed using NTL::RR at 1000-bit precision.

Implementation

Traditional methods of evaluating the beta function either involve evaluating the gamma functions directly, or taking logarithms and then exponentiating the result. However, the former is prone to overflows for even very modest arguments, while the latter is prone to cancellation errors. As an alternative, if we regard the gamma function as a white-box containing the Lanczos approximation, then we can combine the power terms:

which is almost the ideal solution, however almost all of the error occurs in evaluating the power terms when a or b are large. If we assume that a > b then the larger of the two power terms can be reduced by a factor of b, which immediately cuts the maximum error in half:

This may not be the final solution, but it is very competitive compared to other implementation methods.

lbeta takes the logarithm of this form term by term. With a >= b, c = a + b and xgh = x + g - 1\/2, and since agh = cgh - b, the two power terms merge:

log B(a, b) = log(L(a) L(b) / L(c)) + (a - 1/2) log1p(-b / cgh) + b log(bgh / cgh) + (1 - log bgh) / 2

where L is the scaled Lanczos sum. None of these terms cancel when one argument is large, and none underflow when both are. In double precision the result is within 2 ulps where |log B| >= 1, and within 1 epsilon absolute error elsewhere, for arguments across the whole exponent range. Without a Lanczos approximation approximation, lbeta takes logarithms in the generic method below, except that when one argument is small it uses lgamma and tgamma_delta_ratio directly, and when both lie between 1 and the threshold for Stirling's approximation it takes the logarithm of beta itself, which cannot overflow or underflow there.

The generic implementation - where no Lanczos approximation approximation is available - is implemented in a very similar way to the generic version of the gamma function by means of Stirling's approximation. Again in order to avoid numerical overflow the power terms that prefix the series are collected together

There are a few special cases worth mentioning:

When a or b are less than one, we can use the recurrence relations:

to move to a more favorable region where they are both greater than 1.

In addition: