beta function returns imprecise results

[rawls238]

I'm noticing that the invocation of the beta function is returning imprecise results at higher decimal places. For instance, invoking beta(0.2, 0.25) returns 8.45704009347458 in Julia, while on Wolfram Alpha invoking the same thing it returns 8.4570400934745837. I've noticed this on several other invocations. @kmsquire did some digging into this and noted there appears to be a discrepancy in the last few digits of the invocation of lgamma_r though the lgamma_r function can return irrational numbers so in some cases that may be expected behavior. In any case, this issue is for fixing the results of the beta function

[yuyichao]

This difference is 2ulp

julia> (8.4570400934745837 - 8.45704009347458) / eps(8.45704009347458)
2.0

Do we (or libm or openspecfun) provide this level of precision garentee?

[simonbyrne]

beta is implemented in Julia: it is possible to get greater error by choosing more extreme values, see #4301 (comment).

There is an open PR #14349.

[CarloLucibello]

The problem for the case considered, beta(0.2,0.25), is in the gamma function, on which the beta relies. I did some comparison with the GSL implementation of the gamma:

julia> x=0.2; (gamma(x)-4.5908437119988025898)/eps(gamma(x))
0.0
julia> x=0.4; (gamma(x)-2.218159543757686869)/eps(gamma(x))
3.0
julia> x=0.6; (gamma(x)-1.4891922488128179758)/eps(gamma(x))
-4.0
julia> x=1.6;(gamma(x)-0.89351534928769116295)/eps(gamma(x))
-8.0
julia> x=5.6;(gamma(x)-61.553915006289244616 )/eps(gamma(x))
-1.0

A three way comparison between julia, GSL and Mathematica gives

In: Gamma[1.6]-0.89351534928769116295 (*GSL*)
Out: -9.99201*10^-16

In:Gamma[1.6]-0.8935153492876903 (*Julia*)
Out:-1.11022*10^-16

Since julia's gamma calls cmath's gamma, and since the "error" is so small, I think this issue could be closed

[CarloLucibello]

beta is implemented in Julia: it is possible to get greater error by choosing more extreme values, see #4301 (comment).

That issue, for beta(1e10,1e-10), is of a different nature and it will be fixed in #14349

[andreasnoack]

I agree. There is only an issue for large arguments and that is covered by #14349 so I'll close this one.