' Gamma function
DECLARE FUNCTION Gamma# (Z#)
DECLARE FUNCTION LnGamma# (BYVAL Z#)
FOR X# = .1 TO 2.05 STEP .1
  PRINT USING "##.#  "; X#;
  PRINT USING "#.#########"; Gamma#(X#)
NEXT X#
END

FUNCTION Gamma# (Z#)
  Gamma# = EXP(LnGamma#(Z#))
END FUNCTION

FUNCTION LnGamma# (BYVAL Z#)
  CONST PI# = 3.14159265358979#
  DIM LZ#(6)
  LZ#(0) = 1.00000000019001#
  LZ#(1) = 76.1800917294715#
  LZ#(2) = -86.5053203294168#
  LZ#(3) = 24.0140982408309#
  LZ#(4) = -1.23173957245015#
  LZ#(5) = 1.2086509738662D-03
  LZ#(6) = -.000005395239385#
  IF Z# < .5 THEN
    LnGamma# = LOG(PI# / SIN(PI# * Z#)) - LnGamma#(1# - Z#)
  ELSE
    Z# = Z# - 1#
    B# = Z# + 5.5
    A# = LZ#(0)
    FOR I% = 1 TO 6
      A# = A# + LZ#(I%) / (Z# + I%)
    NEXT I%
    LnGamma# = (LOG(SQR(2 * PI#)) + LOG(A#) - B#) + LOG(B#) * (Z# + .5)
  END IF
END FUNCTION