code 35061;
real procedure GAMMA(X); value X; real X;
begin real Y, S, F, G, ODD, EVEN;
    boolean INV;
    if X < .5 then 
    begin Y:= X - ENTIER(X / 2) * 2; S:= 3.14159 26535 8979;
        if Y >= 1 then begin S:= - S; Y:= 2 - Y end;
        if Y >= .5 then Y:= 1 - Y; INV:= true; X:= 1 - X;
        F:= S / SIN(3.14159 26535 8979 * Y)
    end 
    else INV:= false;
    if X > 22 then G:= EXP(LOG GAMMA(X)) else 
    begin S:= 1;
    NEXT: if X > 1.5 then 
        begin X:= X - 1; S:= S * X; goto NEXT end;
        G:= S / RECIP GAMMA(1 - X, ODD, EVEN)
    end;
    GAMMA:= if INV then F / G else G
end GAMMA

        eop