53 lines
1 KiB
Text
53 lines
1 KiB
Text
(* Gamma function *)
|
|
(* Lanczos approximation *)
|
|
|
|
block
|
|
var
|
|
x: real;
|
|
|
|
unit gamma: function (z: real): real;
|
|
const
|
|
pi = 3.14159265358979;
|
|
var
|
|
lz: arrayof real;
|
|
|
|
unit ln_gamma: function (z: real): real;
|
|
(* Uses lz and pi *)
|
|
var
|
|
i: integer,
|
|
b, a: real;
|
|
begin
|
|
if z < 0.5 then
|
|
result := ln(pi / sin(pi * z)) - ln_gamma(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 + i)
|
|
od;
|
|
result := (ln(sqrt(2.0 * pi)) + ln(a) - b) + ln(b) * (z + 0.5)
|
|
fi;
|
|
end ln_gamma;
|
|
|
|
begin
|
|
array lz dim (0 : 6);
|
|
lz(0) := 1.00000000019001;
|
|
lz(1) := 76.1800917294715;
|
|
lz(2) := -86.5053203294168;
|
|
lz(3) := 24.0140982408309;
|
|
lz(4) := -1.23173957245015;
|
|
lz(5) := 1.2086509738662E-3;
|
|
lz(6) := -0.000005395239385;
|
|
result := exp(ln_gamma(z))
|
|
end gamma;
|
|
|
|
begin
|
|
x:= 0.1;
|
|
while x < 2.05
|
|
do
|
|
writeln(x: 4: 1, " ", gamma(x): 15: 12);
|
|
x:= x + 0.1
|
|
od;
|
|
end
|