function [r,p,c,s] =genz(a,b,B)
% GENZ solves for the equilibrium distribution
% vector given its z-transform S(z).
% The matrix should be lower Hessenberg
% with identical diagonals except for the elements
% in the first column:
%
%      |a0  b0  0  0  0  ...  |
%      |a1  b1  b0 0  0  ...  |
%  P = |a2  b2  b1 b0 0  ...  |
%      |a3  b3  b2 b1 b0 ...  |
%      |.   .   .   .  .      |
%
% a = vector of first column of transition matrix.
% b = vector of other columns of transition matrix.
% B = number of components of distribution vector
% s =[s_0  s_1  s_2 ... s_B-1
% 

% first ensure that the matrix is column stochastic
if (sum(a) ~= 1.0000)
    fprintf('Error sum of a vector is not equal to 1 (sum(a) = %10.5e\n', sum(a));
end

if (sum(b) ~= 1.0000)
    fprintf('Error sum of b vector is not equal to 1 (sum(b) = %10.5e\n', sum(b));
end

a_prime(1) = 0;% Construct A'
a_prime(2) = 0;
for i=3:size(a,2)+1
   a_prime(i) = -(i-2)*a(i-1);
end
a_prime;

b_prime(1) = 0;% Construct B'
b_prime(2) = 0;
for i=3:size(b,2)+1
   b_prime(i) = -(i-2)*b(i-1);
end
b_prime;
a1_prime = sum(a_prime); % sum of terms when z=1
b1_prime = sum(b_prime); % sum of terms when z=1

s0 = (1+b1_prime)/(1+b1_prime-a1_prime);
% Construct S(z) but first ensure that size of [0,a] equals size of b:
a_size = size([0,a]);
b_size = size(b);
length_a = a_size(2);
length_b = b_size(2);
if length_a > length_b
    diff = length_a-length_b;
    padding = zeros(1,diff);
    Nom = s0*([0,a]-[b,padding]); % z^-1 X A is equivalent to a new polynomial [0,a]    
elseif length_b > length_a
    diff = length_b-length_a;
    padding = zeros(1,diff);
    Nom = s0*([0,a,padding]-b); % z^-1 X A is equivalent to a new polynomial [0,a]
else    
    Nom = s0*([0,a]-b); % z^-1 X A is equivalent to a new polynomial [0,a]
end
constant = zeros(1,size(b,2));
constant(2) = 1;
Denom = constant -b;
% Find residues and poles
[r,p,c] = residue(Nom,Denom);
%Find components of s
s(1) = s0;
for n=2:B
   temp=0;
   for i=1:size(p,1)
      temp = temp+r(i)*p(i)^(n-2) ;
   end
   s(n) = temp;
end
s;