require 'mathn'

L     = 5
M     = 4
K     = 3
N     = 2

Delta = 10**-5

data = <<MATRIX
0 1 0 0 0 0 0 0
0 0.5 0.5 0 0 0 0 0
0 0 0 1 0 0 0 0
0 0 0 0 1 0 0 0
0 0 0 0 0 0 1 0
0 0.333 0.333 0 0 0 0.333 0
0 0 0 0 0 0 0.5 0.5
0 0 0 0 0 0 0 1
MATRIX

transitions = data.split("\n").map {|row| row.split.map {|it| eval it }}

class FixedMatrix < Matrix
  def **(d)
    (d == 0 and square?) ? self.class.identity(column_size) : super
  end
end

class Program
  attr_reader :size, :k

  def initialize(data)
    @p      = FixedMatrix.rows data
    @size   = @p.row_size
    calculate_k
  end

  def calculate_k
    @k = 1
    @k += 1 while 1 - (@p**@k)[0, @size-1] > Delta
  end

  def probability_density_for(k)
    (@p**k)[0, @size-1] - (@p**(k-1))[0, @size-1]
  end

  def singular
    @f ||= Array.new(@k) {|k| probability_density_for(k + 1) }
  end

  def sequential(n)
    f  = singular
    fr = f.clone
    k  = fr.size * 2
    2.upto(n) do
      ff = Array.new(k, 0)
      0.upto(k-1) do |i|
        (i>=f.length ? (i-f.length+1) : 0).upto(i>=fr.length ? fr.length-1 : i) {|j| ff[i+1]+=fr[j]*f[i-j]}
      end
      fr = ff.clone
  		k += f.length
    end
    fr
  end

  def parallel_and(n)
    unless n.is_a?(Array)
      f  = singular
      fr = Array.new(f.length,0)
      0.upto(fr.length-1) do |i|
        proda,prodb=1,1
        1.upto(n) do
          suma,sumb=0,0
          0.upto(i) do |j|
            suma+=f[j]
            if j<i : sumb+=f[j]
            end
          end
          proda*=suma
          prodb*=sumb
        end
        fr[i]=proda-prodb
      end
      fr
    else
      fs = n
      fr = Array.new(k,0)
      (@k-1).times do |i|
        a = fs.inject(1) do |prod, f|
          prod * f.take(i+1).inject(0) {|sum, fi| sum + fi }
        end
        b = fs.inject(1) do |prod, f|
          f1 = f.clone
          f1.pop
          prod * f1.take(i).inject(0) {|sum, fi| sum + fi }
        end
        fr[i+1] = a - b
      end
      fr
    end
  end

  def parallel_or(n, f=nil)
    unless n.is_a?(Array)
      f  ||= singular
      fr = Array.new(f.length,0)
      0.upto(fr.length-1) do |i|
        proda,prodb=1,1
        1.upto(n) do
          suma,sumb=0,0
          0.upto(i) do |j|
            suma+=f[j]
            if j<i : sumb+=f[j]
            end
          end
          proda*=1-suma
          prodb*=1-sumb
        end
        fr[i]=prodb-proda
      end
      fr
    else
      f = n
      fr = Array.new(k,0)
      0.upto(fr.length-1) do |i|
        proda,prodb=1,1
        f.size.times do |f_index|
          suma,sumb=0,0
          0.upto(i) do |j|
            suma+=f[f_index][j]
            if j<i : sumb+=f[f_index][j]
            end
          end
          proda*=1-suma
          prodb*=1-sumb
        end
        fr[i]=prodb-proda
      end
      fr
    end
  end

  def binominal(m, l)
    return 0 if m > l
    c = 1
    1.upto(m) do |i|
      c *= l
      c /= i
      c -= 1
    end
    c
  end

  def out_of(m, l)
    return parallel_and l if m == l
    fr = parallel_and m
    bc = binominal m, l
    parallel_or bc, fr
  end
end

def dynamic_characteristics_for(f)
  m, d = 0, 0
  f.each_with_index {|fk, k| m += fk*(k+1) }
  f.each_with_index {|fk, k| d += fk*(k+1-m)**2 }
  e = 1 - f.inject {|sum, el| sum + el }
  [m, d, e]
end

def print_density(header, f)
  puts header
  puts "=================="
  f.each_with_index do |fk, i|
    k = i + 1
    puts "Fk(#{k})=#{fk}"
  end
  puts "==================\n\n"
end

def print_characteristics(header, data)
  m,d,e = *data
  puts header
  puts "Математическое ожидание: #{m},\nДисперсия:#{d},\nОшибка:#{e}\n\n"
end

program     = Program.new transitions
Singular    = program.singular
Out         = program.out_of M, L
And         = program.parallel_and K
Or          = program.parallel_or [Out, And]
Full        = program.parallel_and [Or, Singular]

puts "Параметры системы: L=#{L}, M=#{M}, K=#{K}, N=#{N}\n\n"
print_characteristics "Один процесс:", dynamic_characteristics_for(Singular)
print_characteristics "M из L:", dynamic_characteristics_for(Out)
print_characteristics "И:", dynamic_characteristics_for(And)
print_characteristics "ИЛИ:", dynamic_characteristics_for(Or)
print_characteristics "Полная программа:", dynamic_characteristics_for(Full)