code 41014;
real procedure POISSONINV(PROB, MU, LEFT);
value PROB, MU, LEFT;
real PROB, MU; Boolean LEFT;
begin integer X; real PX, PCUM;

    if PROB ≤ 0 ∨ PROB ≥ 1 then
    STATAL3 ERROR(“POISSONINV”, 1, PROB) else
    if MU ≤ 0 then
    STATAL3 ERROR(“POISSONINV”, 2, MU);
      
    if LEFT then
    begin X:= (PHINV(PROB) / 2 + SQRT(MU)) ⭡ 2 - 1;
        if X < 0 then X:= 0;
        if PROB < EXP(-MU) then POISSONINV:= -1 else
        begin PX:= POISSONPROB(X, MU);
            PCUM:= POISSON(X, MU);
            if PCUM > PROB then
            begin for PCUM:= PCUM - PX
                while PCUM > PROB do
                begin PX:= PX × X / MU; X:= X - 1 end;
                X:=X - 1
            end else
            begin for PX:= PX × MU / (X + 1)
                while PCUM + PX < PROB do
                begin X:= X + 1; PCUM:= PCUM + PX end
            end;
            POISSONINV:= X
        end
    end else
    begin X:= (PHINV(1 - PROB) / 2 + SQRT(MU)) ⭡ 2 + 1;
       if X < 0 then X:= 0;
       PCUM:= 1 - POISSON(X - 1, MU);
       PX:= POISSONPROB(X, MU);
       if PCUM < PROB then
       begin for PX:= PX × X / MU
           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 × MU / (X + 1); X:= X + 1 end;
            X:= X + 1
        end;
        POISSONINV:= X
    end
end POISSONINV;
eop