code 35062;
real procedure LOG GAMMA(X); value X; real X;
if X > 13 then 
begin real R, X2;
    R:= 1;
NEXT: if X <= 22 then 
    begin R:= R / X; X:= X + 1; goto NEXT end;
    X2:= - 1 / (X * X); R:= LN(R);
    LOG GAMMA:= LN(X) * (X - .5) - X + R + .91893 85332 04672 +
    (((.59523 80952 38095"-3  * X2 + .79365 07936 50794"-3) * X2 +
       .27777 77777 77778"-2) * X2 + .83333 33333 33333"-1) / X
end 
else 
begin real Y, F, U0, U1, U, Z;
    integer I;
    array B[1:18];
    F:= 1; U0:= U1:= 0;
    B[ 1]:= -.07611 41616 704358; B[ 2]:= +.00843 23249 659328;
    B[ 3]:= -.00107 94937 263286; B[ 4]:= +.00014 90074 800369;
    B[ 5]:= -.00002 15123 998886; B[ 6]:= +.00000 31979 329861;
    B[ 7]:= -.00000 04851 693012; B[ 8]:= +.00000 00747 148782;
    B[ 9]:= -.00000 00116 382967; B[10]:= +.00000 00018 294004;
    B[11]:= -.00000 00002 896918; B[12]:= +.00000 00000 461570;
    B[13]:= -.00000 00000 073928; B[14]:= +.00000 00000 011894;
    B[15]:= -.00000 00000 001921; B[16]:= +.00000 00000 000311;
    B[17]:= -.00000 00000 000051; B[18]:= +.00000 00000 000008;
    if X < 1 then 
    begin F:= 1 / X; X:= X + 1 end 
    else 
NEXT: if X > 2 then 
    begin X:= X - 1; F:= F * X; goto NEXT end;
    F:= LN(F); Y:= X + X - 3; Z:= Y + Y;
    for I:= 18 step - 1 until 1 do 
    begin U:= U0; U0:= Z * U0 + B[I] - U1; U1:= U end;
    LOG GAMMA:= (U0 * Y + .49141 53930 29387 - U1) * (X - 1) * (X - 2)
    + F
end LOG GAMMA

        eop