1

# -*- coding: utf-8 -*-

def P(n):
    res = 1
    for i in range(0, n):
        if i <= N:
            lamda = lamda_const
        else:
            lamda = 2*lamda_const
        res *= lamda/mu_const
    return res * p0

def P0():
    res1 = 0
    res2 = 1
    for j in range(1, n+1):
        for i in range(n):
            if i <= N:
                lamda = lamda_const
            else:
                lamda = 2*lamda_const
            res2 *= lamda/mu_const
        res1 += res2
    return 1/(1+res1)

N = 4
n = 10
lamda_const = 1.4
mu_const = 9.5
p0 = P0()
#p0 = 1 - lamda_const/mu_const

if 2*lamda_const/mu_const < 1:
    p = [P(j) for j in range(1, n+1)]
    print("Начальное распределение:\n{}\n\nСтационарное распределение:\n{}\n\nСумма вероятностей (проверка):\n{}".format(p0, p, sum(p)))
else:
    print("Условие существования стационарного режима не выполнено")