*FLOAT64
      INSTALL @lib$+"FNUSING"

      FOR x = 0.1 TO 2.05 STEP 0.1
        PRINT FNusing("#.#",x), FNusing("##.############", FNgamma(x))
      NEXT
      END

      DEF FNgamma(z) = EXP(FNlngamma(z))

      DEF FNlngamma(z)
      LOCAL a, b, i%, lz()
      DIM lz(6)
      lz() = 1.000000000190015, 76.18009172947146, -86.50532032941677, \
      \ 24.01409824083091, -1.231739572450155, 0.0012086509738662, -0.000005395239385
      IF z < 0.5 THEN = LN(PI / SIN(PI * z)) - FNlngamma(1.0 - z)
      z -= 1.0
      b = z + 5.5
      a = lz(0)
      FOR i% = 1 TO 6
        a += lz(i%) / (z + i%)
      NEXT
      = (LNSQR(2*PI) + LN(a) - b) + LN(b) * (z+0.5)