function P()
    Sum=0
    Maxisp=5000
    Maxt=30000
    MatrixGammaBin=1:13
    MatrixExpBin=1:20
    MatrixErlangBin=1:25
    MatrixS(Maxt)=0
    MatrixResult(Maxt)=0
    for isp=1:Maxisp
        //Erlang
        TAU=-1
        for k=1:25
            while or([(TAU>30000), or(TAU<0)]) do
                R1=grand(1,1,'unf',0,1)
                R1=(1 - 0.13533528) * R1 + 0.13533528
                R2=grand(1,1,'unf',0,1)
                R2=(1 - 0.13533528) * R2 + 0.13533528
                R3=grand(1,1,'unf',0,1)
                R3=(1 - 0.13533528) * R3 + 0.13533528
                TAU=-(1 / 0.0001) * (log(R1*R2*R3))
            end
            MatrixErlang(k)=round(TAU)
            TAU=-1
        end
        //EndofErlang

        //Exp
        for k=1:20
            R=(grand(1,1,'unf',0,1))
            R = (1 - 0.04978707) * R + 0.04978707
            TAU = -(1 / 0.0001) * log(R)
            MatrixExp(k)=round(TAU)
        end
        //EndofExp

        //Gamma
        R1=9;
        R2=9;
        R3=9;
        TAU=-1;
        MatrixGamma=1:13
        for k=1:13
            while or([(TAU>30000), or(TAU<0)]) do
                R1=grand(1,1,'unf',0,1)
                R1=(1 - 0.13533528) * R1 + 0.13533528
                R2=grand(1,1,'unf',0,1)
                R2=(1 - 0.13533528) * R2 + 0.13533528
                R3=grand(1,1,'unf',0,1)
                R3=(1 - 0.13533528) * R3 + 0.13533528
                TAU=-(1 / 0.0001) * (log(R1) + log(R2) + log(R3))
            end
            MatrixGamma(k)=round(TAU)
            TAU=-1
        end
        //EndofGamma

        for t=1:Maxt
            for i=1:13
                if MatrixGamma(i)>t then
                    MatrixGammaBin(i)=1
                else
                    MatrixGammaBin(i)=0
                end
            end

            for ii=1:20
                if MatrixExp(ii)>t then
                    MatrixExpBin(ii)=1
                else
                    MatrixExpBin(ii)=0
                end
            end

            for iii=1:25
                if MatrixErlang(iii)>t then
                    MatrixErlangBin(iii)=1
                else
                    MatrixErlangBin(iii)=0
                end
            end

            S = 1 - ((1 - MatrixErlangBin(1) * MatrixExpBin(1)) * (1 - MatrixErlangBin(2) * MatrixExpBin(2)) * (1 - MatrixErlangBin(3) * MatrixExpBin(3) * MatrixGammaBin(1)) * (1 - MatrixErlangBin(4) * MatrixExpBin(4) * MatrixGammaBin(2)) * (1 - MatrixErlangBin(5) * MatrixExpBin(5) * MatrixGammaBin(3)) * (1 - MatrixErlangBin(6) * MatrixErlangBin(7) * MatrixExpBin(6) * MatrixGammaBin(4)) * (1 - MatrixErlangBin(8) * MatrixErlangBin(9) * MatrixExpBin(7) * MatrixGammaBin(5)) * (1 - MatrixErlangBin(10) * MatrixErlangBin(11) * MatrixExpBin(8) * MatrixGammaBin(6)) * (1 - MatrixErlangBin(12) * MatrixErlangBin(13) * MatrixExpBin(9) * MatrixGammaBin(7)) * (1 - MatrixErlangBin(14) * MatrixErlangBin(15) * MatrixExpBin(10) * MatrixGammaBin(8)) * (1 - MatrixErlangBin(16) * MatrixErlangBin(17) * MatrixExpBin(11) * MatrixExpBin(12) * MatrixGammaBin(9)) * (1 - MatrixErlangBin(18) * MatrixErlangBin(19) * MatrixExpBin(13) * MatrixExpBin(14) * MatrixGammaBin(10)) * (1 - MatrixErlangBin(20) * MatrixErlangBin(21) * MatrixExpBin(15) * MatrixExpBin(16) * MatrixGammaBin(11)) * (1 - MatrixErlangBin(22) * MatrixErlangBin(23) * MatrixExpBin(17) * MatrixExpBin(18) * MatrixGammaBin(12)) * (1 - MatrixErlangBin(24) * MatrixErlangBin(25) * MatrixExpBin(19) * MatrixExpBin(20) * MatrixGammaBin(13)))
            //TODO: исправить алгоритм суммирования
            MatrixS(t)=MatrixS(t)+S
            Sum=Sum+MatrixS(t)
            //Pt=((1/Maxisp)*Sum)
            //MatrixResult(t)=MatrixResult(t)+Pt
            Sum=0
        end
        print(%io(2),isp)
        if modulo(isp,100)==0 then
            clf()
            plot2d(MatrixS)
        end
    end
endfunction