Add support for very small and very large shape parameter values to betaIncomplete #10889
Description
std.mathspecial.betaIncomplete(a, b, x) computes the value of the regularized, incomplete beta function Ix(a,b). The domain variable x is restricted to the interval [0,1], and the shape parameters a and b must belong to ℝ+. The range of Ix is [0,1]. It is a non-decreasing function where I0(a,b) = 0 and I1(a,b) = 1 for all values of a and b.
When one or both of the shape parameters is very small or very large, betaIncomplete often fails to determine the value of Ix and returns real.nan. I propose that implementation of betaIncomplete be extended to support the degenerated cases of Ix(a,b) and then use interpolation to approximate the values of Ix it can't currently compute.
Most of the limits of Ix(a,b) as a or b ⟶ 0+ or ∞ exist and are called the degenerate cases. Here they are.
- Ix(+0,b) = lima⟶0+ Ix(a,b) = { 0, x=0; 1, 0<x≤1 }.
- Ix(∞,b) = lima⟶∞ Ix(a,b) = { 0, 0≤x<1; 1, x=1 }.
- Ix(a,+0) = limb⟶0+ Ix(a,b) = { 0, 0≤x<1; 1, x=1 }.
- Ix(a,∞) = limb⟶∞ Ix(a,b) = { 0, x=0; 1, 0<x≤1 }.
- Ix(+0,∞) = lima⟶0+,b⟶∞ Ix(a,b) = { 0, x=0; 1, 0<x≤1 }.
- Ix(∞,+0) = lima⟶∞,b⟶0+ Ix(a,b) = { 0, 0≤x<1; 1, x=1 }.
- Ix(+0,+0) = lima,b⟶0+ Ix(a,b) = { 0, x=0; 1, x=1; NaN, 0<x<1 }.
- Ix(∞,∞) = lima,b⟶∞ Ix(a,b) = { 0, x=0; 1, x=1; NaN, 0<x<1 }.