56 lines
1.1 KiB
Text
56 lines
1.1 KiB
Text
|
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.
|