std.mathspecial.beta does not handle all edge cases correctly.

[tedgin]

The beta function B(x,y) has lots of edge cases when approximating real numbers with a floating point type. Assume x and y are real numbers, then B(x,y) can be defined as Γ(x)Γ(y)/Γ(x+y). The function Γ(z) isn't defined when z is a non-positive integer, but it is otherwise continuous, and 1/Γ(z) exists everywhere, being 0 when z is a non-positive integer. This means that whenever x or y is a non-positive integer, B(x,y) doesn't exist. Also, whenever x+y is a non-positive integer Γ(x+y) diverges, 1/Γ(x+y) is 0. In these case, as long as neither x nor y is a non-positive integer, if x+y is a non-positive integer, B(x,y) = 0. Approximating the set of real numbers using a IEEE 754 floating point type further complicates things because 0 isn't represented. Instead, +0.0 and -0.0 are represented where the first is a value infinitesimally larger than 0 and the second is a value infinitesimally smaller than 0. The edge cases are too many to enumerate here, but each of them need to be addressed in the implementation of std.mathspecial.beta. I've verified that many aren't.