26 lines
633 B
Text
26 lines
633 B
Text
# 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
|