import mpmath


def logbeta(x, y):
    """
    Natural logarithm of beta(x, y).

    The beta function is

                     Gamma(x) Gamma(y)
        beta(x, y) = -----------------
                       Gamma(x + y)

    where Gamma(z) is the Gamma function.
    """
    with mpmath.extradps(5):
        return (mpmath.loggamma(x)
                + mpmath.loggamma(y)
                - mpmath.loggamma(mpmath.fsum([x, y])))