Note that the formula as written suffers from numerical instability when evaluated at finite precision. This problem can be avoided by using logarithms of the beta and gamma functions, and then exponentiating. See the source code of function p_half in file diff_exp.py for a complete, numerically stable implementation.
