Hi everyone,

I was comparing a colleague’s MATLAB code to my version (we are coding in parallel, and I am using Julia. I am also trying to convince him to move to Julia. Making progress… ), and I noticed that the beta function returns something different for large x or y. Because of some other related functions these values are fed into, there are some differences that emerge in our results.

I elaborate on this in the following code:

```
using SpecialFunctions
beta1 = 3.6e13;
alpha1 = 0.563483398;
a13 = -1.27e29;
julia> beta1
3.6e13
julia> alpha1
0.563483398
julia> beta(beta1,1/alpha1+1)
3.9891491767885927e-38
# MATLAB's version
>> beta(beta1,1/alpha1+1)
ans = 4.47377930618112e-38
```

The differences emerge in our codes because we use the following expression, which result in different answers.

```
julia> a13
-1.27e29
julia> a13*beta1*beta(2*beta1,1/alpha1+1)
-26651.76910015019
# MATLAB's version
>> a13*beta1*beta(2*beta1,1/alpha1+1)
ans = -21558.4828041656
```

Would anyone have any insight about why these two values differ? I tried reading up on this and found some useful issues like #14349 and #14620 but I am wondering if someone has any additional comments or insights regarding what is likely happening here.