code 35050;
real procedure INCBETA(X,P,Q,EPS);
value X,P,Q,EPS; real X,P,Q,EPS;
begin integer M,N; real G,F,FN,FN1,FN2,GN,GN1,GN2,DN,PQ;
    boolean N EVEN,RECUR;

    if X=0 or X=1 then INCBETA:= X else 
    begin if X>.5 then 
        begin F:= P; P:= Q; Q:= F; X:= 1-X; RECUR:= true end 
        else RECUR:= false;
        G:= FN2:= 0; M:= 0; PQ:= P+Q; F:= FN1:= GN1:= GN2:= 1;
        N EVEN:= false;
        for N:= 1,N+1 while ABS((F-G)/F) > EPS do 
        begin if N EVEN then 
            begin M:= M+1; DN:= M*X*(Q-M)/(P+N-1)/(P+N) end 
            else DN:= -X*(P+M)*(PQ+M)/(P+N-1)/(P+N);
            G:= F; FN:= FN1+DN*FN2; GN:= GN1+DN*GN2;
            N EVEN:= ^ N EVEN; F:= FN/GN;
            FN2:= FN1; FN1:= FN; GN2:= GN1; GN1:= GN
        end;
        F:= F*X**P*(1-X)**Q*GAMMA(P+Q)/GAMMA(P+1)/GAMMA(Q);
        if RECUR then F:= 1-F;
        INCBETA:= F
    end 
end INCBETA;
        eop