# Gamma function

proc ln_gamma(z :: number) :: number is
  lz := reg(1.00000000019001, 76.1800917294715, -86.5053203294168, 24.0140982408309, -1.23173957245015, 0.0012086509738662, -0.000005395239385);
  if z < 0.5 then
    return ln(Pi / sin(Pi * z)) - ln_gamma(1.0 - z)
  else
    z -:= 1.0;
    b := z + 5.5;
    a := lz[1];
    for i from 2 to 7 do
      a +:= lz[i] / (z + i - 1)
    od;
    return (ln(sqrt(2 * Pi)) + ln(a) - b) + ln(b) * (z + 0.5)
  fi
end;

proc gamma(z :: number) :: number is
  return exp(ln_gamma(z));
end;

scope
  for x from 0.1 to 2.05 by 0.1 do
    printf("%3.1f\t%15.12f\n", x, gamma(x))
  od
epocs