code 41010;
real procedure NEGBININV(PROB, K, P, LEFT);
value PROB, K, P, LEFT;
real PROB, K, P; Boolean LEFT;
begin integer X; real PX, PCUM;

    if PROB ≤ 0 ∨ PROB ≥ 1 then
    STATAL3 ERROR(“NEGBININV”, 1, PROB) else
    if K > ENTIER(K) ∨ K < 0 then
    STATAL3 ERROR(“NEGBININV”, 2, K) else
    if P ≤ 0 ∨ P > 1 then
    STATAL3 ERROR(“NEGBININV”, 3, P);
    
    if P = 1 ∨ K = 0 then
        NEGBININV:= (if LEFT then K - 1 else K + 1)
    else if LEFT then
    begin X:= (PHINV(PROB) ×
        SQRT(K × (1 - P)) + K - P / 2) / P;
        if X < K then X:= K;
        if PROB < P ⭡ K then NEGBININV:= K - 1 else
        begin PX:= NEGBINPROB(X, K, P);
            PCUM:= NEGBIN(X, K, P);
            if PCUM > PROB then
            begin for PCUM:= PCUM - PX
                while PCUM > PROB do
                begin PX:= PX × (X - K) / (1 - P)
                    / (X - 1);
                    X:= X - 1
                end; X:= X - 1
            end else
            begin for PX:= PX × (1- P) × X / (X - K + 1)
                while PCUM + PX < PROB do
                begin X:= X + 1; PCUM:= PCUM + PX end
            end;
            NEGBININV:= X
        end
    end else
    begin X:= (PHINV(1 - PROB) ×
        SQRT(K × (1 - P)) + K + P / 2) / P;
        if X > K then X:= K;
        PCUM:= 1 - NEGBIN(X - 1, K, P);
        PX:= NEGBINPROB(X, K, P);
        if PCUM < PROB then
        begin for PX:= PX × (X - K) / (1 - P) / (X - 1)
            while PCUM + PX < PROB do
            begin X:= X - 1; PCUM:= PCUM + PX end
        end else
        begin for PCUM:= PCUM - PX
            while PCUM > PROB do
            begin PX:= PX × (1 - P) × X / (X - K + 1);
                X:= X + 1
            end; X:= X + 1
        end
    end;
    NEGBININV:= X
end NEGBININV;
eop