def Lgamma(x):
    '''return natural logarithm of gamma function of x
    raise ValueError if x is negative integer'''
    x = float(x)

    ##/* purge off +-inf, NaN, +-0, and negative arguments */
    if ((x == inf) or (x == -inf)):
        return inf
##    if (x is nan):
##        return nan
    
    e,ix = frexp(x)
    nadj = 0
    signgamp = 1

    if ((e==0.0) and (ix==0)):
        return inf

    if (ix>1020):
        return inf
    
    if ((e != 0.0) and (ix<-71)):
        if (x<0):
            return -log(-x)
        else:
            return -log(x)

    if (e<0):
        if (ix>52):
            return inf ##one/zero
        t = sin_pi(x)
        if (t==zero):
            ##return inf
            raise ValueError('gamma not defined for negative integer')
        nadj = log(pi/fabs(t*x))
        if (t<zero):
            signgamp = -1
        x = -x

    ##/* purge off 1 and 2 */
    if ((x==2.0) or (x==1.0)):
        r = 0.0
        
    ##/* for x < 2.0 */
    elif (ix<2):
        if(x<=0.9): ##/* lgamma(x) = lgamma(x+1)-log(x) */
            r = -log(x)
            if (x>=0.7316):
                y = one-x
                z = y*y
                p1 = a0+z*(a2+z*(a4+z*(a6+z*(a8+z*a10))))
                p2 = z*(a1+z*(a3+z*(a5+z*(a7+z*(a9+z*a11)))))
                p  = y*p1+p2
                r  += (p-0.5*y)
            elif (x>=0.23164):
                y= x-(tc-one)
                z = y*y
                w = z*y
                p1 = t0+w*(t3+w*(t6+w*(t9 +w*t12))) ##	/* parallel comp */
                p2 = t1+w*(t4+w*(t7+w*(t10+w*t13)))
                p3 = t2+w*(t5+w*(t8+w*(t11+w*t14)))
                p  = z*p1-(tt-w*(p2+y*p3))
                r += (tf + p)
            else:
                y = x
                p1 = y*(u0+y*(u1+y*(u2+y*(u3+y*(u4+y*u5)))))
                p2 = one+y*(v1+y*(v2+y*(v3+y*(v4+y*v5))))
                r += (-0.5*y + p1/p2)
        else:
            r = zero
            if(x>=1.7316):
                y=2.0-x ##/* [1.7316,2] */
                z = y*y
                p1 = a0+z*(a2+z*(a4+z*(a6+z*(a8+z*a10))))
                p2 = z*(a1+z*(a3+z*(a5+z*(a7+z*(a9+z*a11)))))
                p  = y*p1+p2
                r  += (p-0.5*y)
            elif(x>=1.23164):
                y=x-tc ##/* [1.23,1.73] */
                z = y*y
                w = z*y
                p1 = t0+w*(t3+w*(t6+w*(t9 +w*t12))) ##	/* parallel comp */
                p2 = t1+w*(t4+w*(t7+w*(t10+w*t13)))
                p3 = t2+w*(t5+w*(t8+w*(t11+w*t14)))
                p  = z*p1-(tt-w*(p2+y*p3))
                r += (tf + p)
            else:
                y=x-one
                p1 = y*(u0+y*(u1+y*(u2+y*(u3+y*(u4+y*u5)))))
                p2 = one+y*(v1+y*(v2+y*(v3+y*(v4+y*v5))))
                r += (-0.5*y + p1/p2)
    ##/* x < 8.0 */
    elif(ix<4):
        i = int(x)
        t = zero
        y = x-i
        p = y*(s0+y*(s1+y*(s2+y*(s3+y*(s4+y*(s5+y*s6))))))
        q = one+y*(r1+y*(r2+y*(r3+y*(r4+y*(r5+y*r6)))))
        r = half*y+p/q
        z = one ##/* lgamma(1+s) = log(s) + lgamma(s) */
        while (i>2):
            i -=1
            z *= (y+i)
        r += log(z)

    ##/* 8.0 <= x < 2**58 */
    elif (ix<58):
        t = log(x)
        z = one/x
        y = z*z
        w = w0+z*(w1+y*(w2+y*(w3+y*(w4+y*(w5+y*w6)))))
        r = (x-half)*(t-one)+w

    ##/* 2**58 <= x <= inf */
    else:
        r =  x*(log(x)-one)

    if (e<0):
        r = nadj - r
    return signgamp*r
return natural logarithm of gamma function of x
    raise ValueError if x is negative integer
recipes/Python/576393_gamma_functions/recipe-576393.py