import mpmath


def interval_prob(x1, x2, k, theta):
    """
    Compute the probability of x in [x1, x2] for the log-gamma distribution.

    Mathematically, this is the same as

        loggamma.cdf(x2, k, theta) - loggamma.cdf(x1, k, theta)

    but when the two CDF values are nearly equal, this function will give
    a more accurate result.

    x1 must be less than or equal to x2.

    k is the shape parameter of the gamma distribution.
    theta is the scale parameter of the log-gamma distribution.
    """
    if x1 > x2:
        raise ValueError('x1 must not be greater than x2')

    with mpmath.extradps(5):
        x1 = mpmath.mpf(x1)
        x2 = mpmath.mpf(x2)
        k = mpmath.mpf(k)
        theta = mpmath.mpf(theta)
        z1 = x1/theta
        z2 = x2/theta
        return mpmath.gammainc(k, mpmath.exp(z1), mpmath.exp(z2),
                               regularized=True)