std.mathspecial.beta(x,y) generates answer with wrong sign sometimes when x is negative and x+y is large

[tedgin]

When x+y is larger than some value std.mathspecial.beta(x,y), uses exp(logGamma(x) + logGamma(y) - logGamma(x+y)) to compute B(x,y) instead of gamma(x) * gamma(y) / gamma(x+y), because Γ(z) grows large very quickly as z increases. logGamma(z) is equivalent to ln|Γ(z)|. Notice that e^ln|Γ(z)| = |Γ(z)|, so beta loses the sign information of the result when it uses exp and logGamma to calculate B.

Here's an example of computing B(-0.5, 10,000). Looking a graph of B, we can see that the value should be negative, but when we use beta to compute it, we get a positive number.

? cat beta_example.d
import std.mathspecial;
import std.stdio;

void main() { writeln(beta(-0.5, 10_000)); }

? rdmd ./beta_example.d
354.477

std.mathspecial.sgnGamma can be used to fix this issue.