#include<iostream>

using namespace std;

struct Q
{
    double MS, Tp1, Tp2;

};

int n[2] = { 0, 0 }, Dis;
Q *Q1, *Q2;

double genExp(double lambda)
{
    double u, x;

    u = (double)rand() / (RAND_MAX + 1);
    x = -1 / lambda*log(1 - u);

    return x;
}

double genGauss(double medie, double sigma)
{
    double s = 0; int i;

    for (i = 1; i <= 12; i++)
        s += (double)rand() / RAND_MAX;

    return medie + sigma*(s - 6);
}


int evUrmator()
{
    if (n[1] == 0 && n[2] == 0)
        return 1;
    if (n[2] == 0)
    {
        if (Dis <= Q1[1].Tp1) return 1;
        if (Q1[1].Tp2 == 0) return 2;
        return 4;
    }
    if (n[1] == 0)
    {
        if (Dis <= Q2[1].Tp2) return 1;
        if (Q1[1].Tp1 == 0) return 3;
        return 5;
    }
    if (Dis<Q1[1].Tp1 && Dis<Q2[1].Tp2)
        return 1;
    if (Q1[1].Tp1<Q2[1].Tp2)
    {
        if (Q1[1].Tp2 == 0) return 2;
        return 4;
    }
    //if(Q2[1].Tp2<=Q1[1].Tp1)
    else
    {
        if (Q2[1].Tp1 == 0) return 3;
        return 5;
    }
    return 0;
}



int main()
{
    double lambda, DS, NS = 1000000;
    double miu[2] = { 10, 10 };
    double STR = 0, CT = 0, ceas = 0, Tp[10], k,TRm;

    cout << "Lambda = ";
    cin >> lambda;

    for (int i = 1; i <= 2; i++)
    {
        Tp[i] = genExp(miu[i]);
    }
    Q1[1] = { ceas, Tp[1], Tp[2] };

    Dis = genExp(lambda);
    do
    {

        k = evUrmator();
        if (k == 1)
        {
            ceas += Dis;

            if (n[1] > 0)
            {
                Q1[1].Tp1 -= Dis;
            }
            if (n[2] > 0)
            {
                Q2[1].Tp2 -= Dis;
            }

            Tp[1] = genExp(miu[1]);
            Tp[2] = genExp(miu[2]);

            if (n[1] <= n[2])
            {
                n[1]++;
                Q1[n[1]] = { ceas, Tp[1], Tp[2] };
            }
            else
            {
                n[2]++;
                Q2[n[2]] = { ceas, Tp[1], Tp[2] };
            }

            Dis = genExp(lambda);

        }
        else if (k == 2)
        {
            ceas += Q1[1].Tp1;
            Q2[1].Tp2 -= Q1[1].Tp1;
            Dis -= Q1[1].Tp1;

            STR += ceas - Q1[1].MS;
            CT++;

            n[1]--;
            for (int i = 1; i <= n[1]; i++)
            {
                Q1[i] = Q1[i + 1];
            }
        }
        else if (k == 3)
        {
            ceas += Q2[1].Tp2;
            Q1[1].Tp1 -= Q2[1].Tp2;
            Dis -= Q2[1].Tp2;

            STR += ceas - Q2[1].MS;
            CT++;

            n[2]--;
            for (int i = 1; i <= n[2]; i++)
            {
                Q2[i] = Q2[i + 1];
            }
        }
        else if (k == 4)
        {
            ceas += Q1[1].Tp1;
            Q2[1].Tp2 -= Q1[1].Tp1;
            Dis -= Q1[1].Tp1;

            Q1[1].Tp1 = 0;

            n[2]++;

            Q2[n[2]] = Q1[1];

            n[1]--;

            for (int i = 1; i <= n[1]; i++)
            {
                Q1[i] = Q1[i + 1];
            }
        }
        else if (k == 5)
        {
            ceas += Q2[1].Tp2;
            Q1[1].Tp1 = Q2[1].Tp2;
            Dis -= Q2[1].Tp2;

            Q2[2].Tp2 = 0;
            n[1]++;
            Q1[n[1]] = Q2[1];

            n[2]--;
            for (int i = 1; i <= n[2]; i++)
            {
                Q2[i] = Q2[i + 1];
            }
        }

    } while (ceas < DS);

TRm = STR / NS;

cout << "TRm = " << TRm;

system("Pause");
return 0;
}