Julia: beta(-n,n) should return ±1/n, not NaN

[stevengj]

In #14165, @CarloLucibello noticed a bug in our beta function: beta(-n,n) is returning NaN for integers n, whereas the correct answer is (-1)^n / n, I think.

[stevengj]

scipy.special.beta gives the correct answer here.

[stevengj]

It would be easy to simply check for integer arguments. However, it is important to make sure that we also correctly handle any arguments _close_ to (-n,n).

[stevengj]

The scipy beta implementation is from Cephes, and includes some special-case code that we are missing:

• It has special handling for negative-integer arguments to avoid NaN.
• It has special handling of beta(a,b) for |a| >> |b| to avoid cancellation errors.
• For smallish arguments, it calls gamma directly rather than exponentiating lgamma; I suspect this is mostly a performance optimization.

[andreasnoack]

What is the defintion for the beta function with negative (integer) first argument? It appears to me that the usual integral definition doesn't converge for e.g. (-2,2). Negative integer arguments for x are allowed in Mathematica but I couldn't find any explanation or definition in their documentation. Abramowitz and Stegun as well as Wikipedia require positive real part of both arguments.

[CarloLucibello]

It is defined through the gamma functions I think, for negative x and y the integral representation doesn't hold anymore.

If you agree on this I could create a PR addressing this bug mimicking the scipy implementation

[stevengj]

If you take the limit of the gamma function definition as (z,n) approaches (-n,n), you get a finite and well-defined result. Besides being unambiguous, this seems to be totally standard; e.g. it is what Mathematica returns.

We should definitely implement the tricks from the Cephes implementation, which is well-respected.

[andreasnoack]

I have no problem with extending the beta function, but

...this seems to be totally standard

is not a fair description. Cephes doesn't support negative arguments. Support for that was a relatively recent modification. Neither does the Naval Surface Warfare Center library and, as mentioned, Abramowitz and Stegun doesn't define the function for negative arguments.

[stevengj]

@andreasnoack, thanks for the clarification, I didn't realize it was a recent modification to Cephes.

Maple also defines beta for negative integers by the limit.

[andreasnoack]

Thanks for the link. In contrast to Mathematica, Maple's documentation states the definition for negative integers. It would be great if we included a similar comment about negative integer x in our documentation.

[simonbyrne]

I mentioned the |a| >> |b| problem here.

I remember looking at the Cephes code before, and was a bit underwhelmed, but there probably are some good ideas there. Also, I never managed to find its licence status (though presumably whomever has copyright presumably doesn't care that much).

[stevengj]

@simonbyrne, apparently the Cephes author gave permission to distribute it under 3-clause BSD; see e.g. here. Also, apparently Debian inquired about this 11 years ago and eventually got 3-clause BSD permission: the /usr/share/doc/python-scipy/copyright file on my Debian machine lists:

Files: scipy/special/cephes/*
Copyright: 1984, 1987, 1989, 2000 Stephen L. Moshier
License: BSD-3-clause

So, as long as we use the Cephes code from SciPy, we should be okay.

[simonbyrne]

Good to know, thanks.

[Peiffap]

Perhaps this issue can be closed, as the functionality for the beta function is now in the SpecialFunctions.jl package.

There is also an issue there that seems to be about the very same problem, https://github.com/JuliaMath/SpecialFunctions.jl/issues/149. If this one does get closed, it might be worth mentioning there what was mentioned on this issue.