code 41013;
real procedure POISSON(X, MU);
value X, MU; real X, MU;
begin integer IX;
    real procedure KSI(K, L);
    value K, L; real L; integer K;
    begin real U, U2, W; W:= SQRT(L);
         U:= 2 × (SQRT(K + 1) - W);
         U2:= U × U; KSI:=
         U + (U2 - 4) / 12 / W + (-U2 × U + 10 × U) / 72 /
         L + (21 × U2 × U2 - 371 × U2 - 52) / 6480 / L / W
    end;
    
    IX:= ENTIER(X);
    if IX < 0 then POISSON:= 0 else
    if MU ≤ 0 then
        POISSON:= STATAL3 ERROR(“POISSON”,2,MU)
    else
    if MU > 1000 then POISSON:= PHI(KSI(IX, MU))
    else
    begin integer I, MODUS; real MODUSPROB, PROB, CUM;
        MODUS:= ENTIER(MU) + 1; if IX < MODUS then
        begin PROB:= CUM:= POISSONPROB(IX, MU);
            for I:= IX step -1 until 1 do
            begin PROB:= PROB × I / MU;
                CUM:= CUM + PROB
            end
        end else
        begin MODUSPROB:= PROB:= CUM:=
            POISSONPROB(MODUS, MU);
            for I:= MODUS step -1 until 1 do
            begin PROB:= PROB × I / MU;
                CUM:= CUM + PROB
            end;  PROB:= MODUSPROB;
            for I:= MODUS + 1 step 1 until IX do
            begin PROB:= PROB × MU / I;
                CUM:= CUM + PROB
            end
        end;
        POISSON:= CUM
    end
end POISSON;
eop