code: def lgammln(xx):
      
  
      coeff = [76.18009173, -86.50532033, 24.01409822, -1.231739516,
               0.120858003e-2, -0.536382e-5]
      x = xx - 1.0
      tmp = x + 5.5
      tmp = tmp - (x+0.5)*math.log(tmp)
      ser = 1.0
      for j in range(len(coeff)):
          x = x + 1
          ser = ser + coeff[j]/x
      return -tmp + math.log(2.50662827465*ser)
docs: Returns the gamma function of xx.
      Gamma(z) = Integral(0,infinity) of t^(z-1)exp(-t) dt.
  (Adapted from: Numerical Recipies in C.)
  
  Usage:   lgammln(xx)