function value=eig_bisect(d,e,cycle);
% Syntax: value=eig_bisect(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.
%
% 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 precision of eigenvalue computation in eig_bisect
%     the max. error of eigenvalues is 4/2^cycle. 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,1);		% + one dummy element  
value=[-2; 2 * ones(dim,1)];	% + one dummy element
t_array=4 ./ ( 2 .^ [1:cycle+1]');
d=[d;0]; e=[e;0];		 % + one 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);
    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]
    
%--- nice piece of algorithm:
%    if si=j+1, just shift the cirrent eigenvalue
%       si<j+1, do nothing (except for the counter)
%       si>j+1, 

    if si >= (j+1),	% read: if the no. of ev. in the upper interval
    			% is greater or eq. than the index of actually 
    			% calculated ev.
      if value(si+1) > value(j+1),	
      	value(si+1) = value (j+1);
      	counter(si+1) = counter(j+1) + 1;
      end
      
      value(j+1) = value(j+1) + t_array(counter(j+1)+1); % this ev. in the
      			% upper interval, we can increment
    end    % de if
    
    counter(j+1) = counter(j+1) + 1;
  
  end  % de while
end  % de for j

value=value(1:dim);