code 41513;
real procedure GAMMA(X, ALPHA, SCALE);
    value X, ALPHA, SCALE; real X, ALPHA, SCALE;
begin integer DELTA, UPP;
    real BETA, START, SUM, TERM;
    
      real procedure INCGAM(X, A, EPS);
      value X, A, EPS; real X, A, EPS;
      begin real C0, C1, C2, D0, D1, D2, X2, AX, P,
          Q, R, S, R1, R2, SCF; integer N;
          S:= EXP(-X + A × LN(X)); SCF:= 10+300;
          if X ≤ 1 then
          begin X2:= X × X; AX:= A × X;
               D0:= 1; P:= A; C0:= S;
               D1:= (A + 1) × (A + 2 - X);
               C1:= (D1 + AX + 2 × X) × S;
               R2:= C1 / D1;
               for N:= 1, N + 1
               while ABS((R2 - R1) / R2) > EPS do
               begin P:= P + 2;
                    Q:= (P + 1) × (P × (P + 2) - AX);
                    R:= N × (N + A) × (P + 2) × X2;
                    C2:= (Q × C1 + R × C0) / P;
                    D2:= (Q × D1 + R × D0) / P;
                    R1:= R2; R2:= C2 / D2;
                    C0:= C1; C1:= C2; D0:= D1; D1:= D2;
                    if ABS(C1) > SCF ∨ ABS(D1) > SCF then
                    begin C0:= C0 / SCF; C1:= C1 / SCF;
                        D0:= D0 / SCF; D1:= D1 / SCF
                    end
               end; INCGAM:= R2 / A / EXP(LOGGAMMA(A))
          end else
          begin C0:= A × S; C1:= (1 + X) × C0;
               Q:= X + 2 - A;
               D0:= X; D1:= X × Q; R2:= C1 / D1;
               for N:= 1, N + 1
               while ABS((R2 - R1) / R2) > EPS do
               begin Q:= Q + 2; R:=N × (N + 1 - A);
                   C2:= Q × C1 - R × C0; D2:= Q × D1 - R × D0;
                   R1:= R2; R2:= C2 / D2;
                   C0:= C1; C1:= C2; D0:= D1; D1:= D2;
                   if ABS(C1) > SCF ∨ ABS(D1) > SCF then
                   begin C0:= C0 / SCF; C1:= C1 / SCF;
                       D0:= D0 / SCF; D1:= D1 / SCF
                   end
               end; INCGAM:= 1 - R2 / A / EXP(LOGGAMMA(A))
          end
      end INCGAM;
      
      if ALPHA ≤ 0 then
          STATAL3 ERROR(“GAMMA”, 2, ALPHA)
      else if SCALE ≤ 0 then
          STATAL3 ERROR(“GAMMA”, 3, SCALE)
      else if X ≤ 0 then GAMMA:= 0 else
      if ALPHA ≥ 500 then
      begin
          GAMMA:= PHI(((X / SCALE / ALPHA) ⭡ .33333333353333
                   - 1 + 1 / (9 × ALPHA)) × 3 × SQRT(ALPHA))
      end else
      begin X:= X / SCALE; BETA:= ALPHA - ENTIER(ALPHA) + 1;
          START:= if X ≥ 40 then 1 else
                  INCGAM(X, BETA, 10-12);
          if ALPHA < 1 then
              GAMMA:= START + EXP(-X + ALPHA × LN(X)
                      - LOGGAMMA(ALPHA + 1))
          else if ALPHA < 2 then GAMMA:= START
          else if X > 700 then GAMMA:= 1
          else
          begin UPP:= ENTIER(ALPHA) - 2; SUM:= TERM:=
              EXP(-X + (ALPHA - 1) × LN(X) - LOGGAMMA(ALPHA));
              for DELTA:= 1 step 1 until UPP do
              begin TERM:= TERM × (ALPHA - DELTA) / X;
                  SUM:= SUM + TERM
              end;
              GAMMA:= START - SUM
          end
      end
end GAMMA;
eop