def gamma_below(a, x):
    """Normalized lower incomplete gamma integral.

    ``(1/\Gamma(a)) \int_0^x e^{-t} t^{a - 1} dt``

    gamma_below is the complement of gamma_above:

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

    As x grows, gamma_below(a, x) converges to 1.

    For x <= max(1, a), this is computed by the power series[1]:

         x^a e^-x   /           x^2            x^3           \
        ----------- | 1 + x + ------- + -------------- + ... |
        a \Gamma(a) \         (a + 1)   (a + 1)(a + 2)       /

    For x > max(1, a), this is computed by ``1 - gamma_above(a, x)``.


    [1] NIST Digital Library of Mathematical Functions, Release 1.0.9
    of 2014-08-29, Eq. 8.7.1 <http://dlmf.nist.gov/8.7.E1>.

    In the NIST DLMF notation, ``gamma_below(a, x)`` is ``P(a, x)``,
    related to ``\gamma^*`` by ``\gamma^*(a, x) = x^{-a} P(a, x)``.

    To derive the power series, multiply (8.7.1) by ``x^a``, expand
    ``\Gamma(a + k + 1)`` into ``\Gamma(a) a(a + 1)...(a + k)``, and
    factor ``a \Gamma(a)`` out of the sum.
    """
    assert 0. < a               # XXX NaN?
    assert 0. <= x              # XXX NaN?
    if x == 0.:
        return 0.
    if x > max(1., a):
        return 1. - gamma_above(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 seq():
        k = 0
        t = 1.
        while True:
            yield t
            k += 1
            t *= x/(a + k)

    return (m/a)*abs_summation(seq())
dataset/ETHPy150Open probcomp/bayeslite/src/math_util.py/gamma_below/operator `<` -> `>=`