PROGRAM GammaTest;
{\$N+}
{\$E+}
VAR
Quit:BOOLEAN;
z:EXTENDED;

FUNCTION Gamma(x:EXTENDED):EXTENDED;
{ Compute gamma to a fair number of significant digits. Written by
Tom Hardy 16 March 1991 }
CONST
Big=1000.0;
MaxReal=1.0E100;

FUNCTION NearZero(x:EXTENDED):BOOLEAN;
BEGIN { NearZero }
NearZero:=(ABS(x)<1.0E-18)
END; { NearZero }

FUNCTION NearInteger(x:EXTENDED):BOOLEAN;
BEGIN { NearInteger }
NearInteger:=NearZero(FRAC(x)) OR NearZero(1-ABS(FRAC(x)))
END; { NearInteger }

FUNCTION Gamma2(z:EXTENDED):EXTENDED;
{ Taylor expansion for â(z+3). From Mathematical Functions and Their
Approximations by Yudell L. Luke, Academic Press 1975 }
CONST
Max=41;
b:ARRAY [0.. Max] OF EXTENDED=(2.00000000000000000000,
1.84556867019693427879,
1.24646499595134652897,
0.57499416892061222755,
0.23007494075411406302,
0.07371504661602386878,
0.02204110936751696733,
0.00544875407582030942,
0.00135522086023943520,
0.00026478566304549638,
0.00006120306281920073,
0.00000850557917488135,
0.00000240617724013144,
0.00000008802390990648,
0.00000011422276453422,
-0.00000001631475210083,
0.00000000862349738998,
-0.00000000244110402524,
0.00000000087291506387,
-0.00000000028390295131,
0.00000000009560889805,
-0.00000000003178270013,
0.00000000001061093470,
-0.00000000000353684281,
0.00000000000117938189,
-0.00000000000039318677,
0.00000000000013108214,
-0.00000000000004369853,
0.00000000000001456734,
-0.00000000000000485607,
0.00000000000000161876,
-0.00000000000000053961,
0.00000000000000017987,
-0.00000000000000005996,
0.00000000000000001999,
-0.00000000000000000666,
0.00000000000000000222,
-0.00000000000000000074,
0.00000000000000000025,
-0.00000000000000000008,
0.00000000000000000003,
-0.00000000000000000001);
VAR
n:INTEGER;
Sum,Powerz:EXTENDED;
BEGIN { Gamma2 }
Sum:=0.0;
Powerz:=1.0;
FOR n:=0 TO Max DO
BEGIN
Sum:=Sum+b[n]*Powerz;
Powerz:=Powerz*z
END;
Gamma2:=Sum
END; { Gamma2 }

FUNCTION Gamma1(z:EXTENDED):EXTENDED;
{ If z>=k then map z into the range [1-k,k] }
CONST
k=1.5;
VAR
n,i:INTEGER;
Product:EXTENDED;
BEGIN { Gamma1 }
IF (z>=k) THEN
BEGIN
Product:=1.0;
n:=TRUNC(z-k)+1;
FOR i:=-n TO -1 DO
Product:=Product*(i+z);
Gamma1:=Gamma2(z-n)*Product/((z-n)*(z-n+1)*(z-n+2))
END
ELSE
Gamma1:=Gamma2(z)/(z*(1+z)*(2+z))
END; { Gamma1 }

BEGIN { Gamma }
IF ((x>Big) OR (x<=-Big)) THEN
BEGIN
WRITELN('Invalid gamma function argument');
Halt(1)
END
ELSE
IF (x<=0.0) THEN
BEGIN
IF NearInteger(x) THEN
Gamma:=MaxReal
ELSE
Gamma:=Pi/(Gamma1(1-x)*SIN(Pi*x))
END
ELSE
Gamma:=Gamma1(x)
END; { Gamma }

PROCEDURE GetReal(Str:STRING;VAR x:EXTENDED;VAR Quit:BOOLEAN);
VAR
Valid:BOOLEAN;
Code:INTEGER;
s:STRING;
BEGIN { GetReal }
REPEAT
WRITE(Str);
Quit:=(s='');
Val(s,x,Code);
Valid:=(Code=0)
UNTIL (Valid OR Quit)
END; { GetReal }

BEGIN { GammaTest }
WRITELN('Compute the gamma function to a fair number of significant digits');
WRITELN;
GetReal('Enter z, [ENTER] to quit : ',z,Quit);
WHILE (NOT Quit) DO
BEGIN
WRITELN('â(z)= ',Gamma(z));
WRITELN;
GetReal('Enter z, [ENTER] to quit : ',z,Quit);
END
END. { GammaTest }