I don’t have a Julia session handy, but I think the right way to solve this is with asymptotic approximations. The beta function beta(p,q) = gamma(p)*gamma(q)/gamma(p+q). For large p and small q, the ratio gamma(p)/gamma(p+q) is given by an asymptotic expansion of Tricomi and Erdelyi (see equation 1). For really large p, you can use gamma(p)/gamma(p+q) roughly equal to p^-q. The remaining term gamma(q) poses no problem for 1/alpha1+1 in your example.