MODULE GammaFun;
(* Gamma function *)
(* Lanczos approximation *)
FROM LongMath IMPORT
  exp, ln, pi, sin, sqrt;
FROM SLongIO IMPORT
  WriteFixed;
FROM STextIO IMPORT
  WriteLn, WriteString;

VAR
  X: LONGREAL;

PROCEDURE Gamma(Z: LONGREAL): LONGREAL;
TYPE
  LZArr = ARRAY [0 .. 6] OF LONGREAL;
VAR
  LZ: LZArr;

  PROCEDURE LnGamma(Z: LONGREAL): LONGREAL;
  (* Uses LZ *)
  VAR
    B, A: LONGREAL;
    I: CARDINAL;
  BEGIN
    IF Z < 0.5 THEN
      RETURN ln(pi / sin(pi * Z)) - LnGamma(1.0 - Z)
    ELSE
      Z:= Z - 1.0;
      B:= Z + 5.5;
      A:= LZ[0];
      FOR I:= 1 TO 6 DO
        A:= A + LZ[I] / (Z + LFLOAT(I))
      END;
      RETURN (ln(sqrt(2. * pi)) + ln(A) - B) + ln(B) * (Z + 0.5)
    END;
  END LnGamma;

BEGIN
  LZ:= LZArr{1.00000000019001, 76.1800917294715, -86.5053203294168,
             24.0140982408309, -1.23173957245015, 1.2086509738662E-3,
             -0.000005395239385};
  RETURN exp(LnGamma(Z))
END Gamma;

BEGIN
  X:= 0.1;
  WHILE X < 2.05 DO
    WriteFixed(X, 1, 4);
    WriteString("  ");
    WriteFixed(Gamma(X), 12, 15);
    WriteLn;
    X:= X + 0.1
  END;
END GammaFun.