def gamma_above(a, x):
    """Normalized upper incomplete gamma integral.

    ``(1/\Gamma(a)) \int_x^\infty e^{-t} t^{a - 1} dt``

    gamma_above is the complement of gamma_below:

    ``gamma_above(a, x) = 1 - gamma_below(a, x)``

    As x goes to zero, gamma_above(a, x) converges to 1.

    For x > max(1, a), this is computed by the continued fraction[1]:

                   1
        -----------------------
                  1 - s
        x + -------------------
                       1
            1 + ---------------
                       2 - s
                x + -----------
                           2
                    1 + -------
                        x + ...

    which is then multiplied by ``x^a e^{-x} / \Gamma(a)``.

    For x <= max(1, a), this is computed by ``1 - gamma_below(a, x)``.

    [1] Abramowitz & Stegun, p. 263, 6.5.31
    """
    assert 0. < a               # XXX NaN?
    assert 0. <= x              # XXX NaN?
    if x <= max(1., a):
        return 1. - gamma_below(a, x)

    # m = \exp [a \log x - x - \log \Gamma(a)] = x^a e^{-x} / \Gamma(a)
    w = a*math.log(x) - x - math.lgamma(a)
    if w < -MAXLOG:
        return 0.
    m = ieee_exp(w)

    def contfrac():
        yield 1, x              # 1/(x + ...), without division.
        i = 0
        while True:
            if (i % 2) == 0:
                yield (i//2) + 1 - a, 1
            else:
                yield (i + 1)//2, x
            i += 1

    return m*limit(convergents(contfrac()))
dataset/ETHPy150Open probcomp/bayeslite/src/math_util.py/gamma_above/original