Complemented incomplete Gamma integral.

    SYNOPSIS:

    double a, x, y, igamc();

    y = igamc( a, x );

    DESCRIPTION:

    The function is defined by::

        igamc(a,x)   =   1 - igam(a,x)

                                inf.
                                   -
                          1       | |  -t  a-1
                    =   -----     |   e   t   dt.
                         -      | |
                        | (a)    -
                                    x

    In this implementation both arguments must be positive.
    The integral is evaluated by either a power series or
    continued fraction expansion, depending on the relative
    values of a and x.

def _igamc(a, x):
    """Complemented incomplete Gamma integral.

    SYNOPSIS:

    double a, x, y, igamc();

    y = igamc( a, x );

    DESCRIPTION:

    The function is defined by::

        igamc(a,x)   =   1 - igam(a,x)

                                inf.
                                   -
                          1       | |  -t  a-1
                    =   -----     |   e   t   dt.
                         -      | |
                        | (a)    -
                                    x

    In this implementation both arguments must be positive.
    The integral is evaluated by either a power series or
    continued fraction expansion, depending on the relative
    values of a and x.
    """
    # Compute  x**a * exp(-x) / Gamma(a)
    ax = math.exp(a * math.log(x) - x - math.lgamma(a))

    # Continued fraction
    y = 1.0 - a
    z = x + y + 1.0
    c = 0.0
    pkm2 = 1.0
    qkm2 = x
    pkm1 = x + 1.0
    qkm1 = z * x
    ans = pkm1 / qkm1
    while True:
        c += 1.0
        y += 1.0
        z += 2.0
        yc = y * c
        pk = pkm1 * z - pkm2 * yc
        qk = qkm1 * z - qkm2 * yc
        if qk != 0:
            r = pk/qk
            t = abs((ans - r) / r)
            ans = r
        else:
            t = 1.0;
        pkm2 = pkm1
        pkm1 = pk
        qkm2 = qkm1
        qkm1 = qk
        if abs(pk) > BIG:
                pkm2 *= BIGINV;
                pkm1 *= BIGINV;
                qkm2 *= BIGINV;
                qkm1 *= BIGINV;
        if t <= MACHEP:
            return ans * ax