Differences in MATLAB's versus Julia's beta(x, y) for large x or y?

General Usage
7

I think Julia’s answer is closer to correct, but with this large of a number it’s easy to blow up or hit a catastrophic cancellation. To take this back one step, a naive definition for beta is

\beta(p, q) = \frac{\Gamma(p) \Gamma(q)}{\Gamma(p+q)}

You wouldn’t want to use this implementation in the real world, though, because those intermediate values are huge. Like 10^{472392288679098} huge. So I’m quite certain this isn’t the implementation in either language.

julia> gamma(3.6e13)
Inf

But, hey, we can opt into higher precision floating point here:

julia> gamma(big(3.6e13))
7.882861104028365451673058070377064242970022825665897715265241972922672893761775e+472392288679098

julia> p = big(3.6e13)
3.6e+13

julia> q = 1/big(0.563483398)+1
2.774675178628776489804929213602311938188089959295862682297487944939957397047059

julia> gamma(p)*gamma(q)/gamma(p+q)
3.989149176788555366691087811957298867890005562542321569109120632926257477147039e-38

julia> beta(p, q)
3.989149176788555366691087811957298867890005562542321569109120643773484257104723e-38

So I’m not a numericist, but I think this is looking pretty good.

11 Likes