def distance(parallax, parallax_error):
    """Computes distance from Gaia parallaxes using the Bayesian method of Bailer-Jones 2015.
    Input: parallax [mas], parallax error [mas]
    Output: distance [pc], 1-sigma uncertainty in distance [pc]
    """
    import numpy as np
    L = 1350
    parallax, parallax_error = parallax / 1000.0, parallax_error / 1000.0
    coeff = np.array([1.0 / L, -2, parallax / parallax_error ** 2, -(1.0 / 
        parallax_error ** 2)])
    g = np.roots(coeff)
    reals = np.isreal(g)
    realsum = np.sum(reals)
    if realsum == 1:
        gd = np.real(g[np.where(reals)[0]])
    elif realsum == 3:
        if parallax >= 0:
            gd = np.min(g)
        elif parallax < 0:
            gd = g[np.where(g > 0)[0]]
    from scipy.optimize import brentq
    rmax = 1000000.0
    rmode = gd[0]
    M = rmode ** 2 * np.exp(-rmode / L) / parallax_error * np.exp(-1.0 / (2 *
        parallax_error ** 2) * (parallax - 1.0 / rmode) ** 2)
    lo = brentq(lambda x: 2 * np.log(x) - x / L - (parallax - 1.0 / x) ** 2 /
        (2 * parallax_error ** 2) + np.log(2) - np.log(M) - np.log(
        parallax_error), 0.001, rmode)
    hi = brentq(lambda x: 2 * np.log(x) - x / L - (parallax - 1.0 / x) ** 2 /
        (2 * parallax_error ** 2) + np.log(2) - np.log(M) - np.log(
        parallax_error), rmode, rmax)
    fwhm = hi - lo
    sigma = fwhm / 2.355
    return gd[0], sigma
