function value=eig_jan(d,e,cycle);
% Syntax: value=eig_jan(d,e,cycle);
%
% searching eigenvalues of a tridiagonal matrice with d on the main diagonal
% unit vector on the upper one and e on the lower one, using iterative
% bisection algorithm.
%
% based on: Saoudi, Boucher, Guyader:
% "A new efficient algorithm to compute the LSP parameters for speech coding" 
% Signal Processing, no. 28, 1992, pp. 201-212.
% - algorithme reviewed by Jan, as the original did not work
% - the orig. algo. WORKS and is faster than this function... 
%
% e is the vector with elements of the main diagonal, dimension M
% d is the vector with elements of 1st subdiagonal, dimension M (1st element
%     is dummy)
% cycle (optional) is the no. of iterations. Default 9.
% value is the vector of eigenvalues, in decreasing order.
%
% - all vectors are COLUMN
%
% The implementation in Matlab is not "vectorized", so it may be very slow

if (nargin==2),
  cycle=9;		% default
end
dim=length(d);
  
%%% init %%%
counter=zeros(dim,1); 
value=[-2 * ones(dim,1)];
t_array=4 ./ ( 2 .^ [1:cycle+1]');   
d=[d;0]; e=[e;0];		 % adding zero because of k-cycle...

%%% cycle %%%
for j=0:dim-1,
  while counter(j+1) < cycle,
    given= value(j+1) + t_array(counter(j+1)+1);
    
  %--- how many eigenvalues are in the interval [given, 2] ? ---
    si=0;			  % no. of signs WITHOUT changes
    pa=1;			  % polynom index-2
    sc=1;			  % polynom index
    pb=d(1) - given;		  % 1st polynom
    pc=pb;
    for k=1:dim,		  % polynom "expansion"
      sb=sc;
      if pc < 0,
        sc = -1;		  % sign of c
      else
        sc = 1;
      end
            
      if (sb * sc) > 0,
        si = si + 1;		  % so if the sign of 2 polynoms was the SAME,
      end			  %  increment si !
      
      pc = (d(k+1) - given) * pb - e(k+1) * pa; %expansion of the polynom
      pa=pb; 
      pb=pc;			  % delaying of the polynoms
    end       % de for k - pol. expanded, si contains the no. of sign
              % NO-changes, hence the no. of eigenvalues in the interval
              % [given,2]

  %--- si<j+1: do nothing (except for counter of j+1)
  %    si=j+1: just shift value(j+1)
  %    si>j+1: shift value(j+1) and all other values from j+2 to si,
  %                modify their counters 
     
    if si>=j+1,
      value(j+1) = value(j+1) + t_array(counter(j+1)+1);   % the (j+1)th
      
      for k=j+1:si-1,
        value (k+1) = value(j+1);
        counter(k+1) = counter(j+1) + 1;
      end  
    end	% de if
    
    counter(j+1) = counter(j+1) + 1;    % toujours
  
  end   % de while
end % of for j