function r = btr2( p, see )
%
% r = btr2( p )      Determines approximations to the roots of the
%                    polynomial p, by using the Brugnano-Trigiante
%                    method.
%
%                    An optional second input parameter, provides a
%                    trace of the execution.
%
% rel. 9.3.04.       
%
% References:
%  - L.Brugnano, D.Trigiante. Polynomial Roots: the Ultimate Answer?
%    Linear Algebra and its Applications, 225 (1995) 207-219.
%  - L.Brugnano.  Numerical Implementation of a New Algorithm for
%    Polynomials with Multiple Roots.
%    Journal of Difference Equations and Applications, 1 (1995) 187-207.
%  - L.Brugnano.  Errata to: Numerical Implementation of a New Algorithm for
%    Polynomials with Multiple Roots.
%    Journal of Difference Equations and Applications, 2 (1996) 101.
%
%
param = [eps^.6    % tolerace for 0
   eps^.45         % tolerance for deflation
   eps^.75         % tolerance for vdiff
   eps^.15];       % tolerance for scal
if nargin==1 
   see = 0; 
else
   disp('tolerance off-diag.entry = '),disp(param(1))   
   disp('tolerance for factor     = '),disp(param(2))   
   disp('tolerance for VDIFF      = '),disp(param(3))
   disp('tolerance for SCAL       = '),disp(param(4))
   pause
end
if sum(isnan(p))>0 | sum(isinf(p))>0, error('polynomial p not valid'), end
zin = sum(cumprod(p==0));
p(1:zin) = [];                   % remove possible leading zero coefficients.
if zin>0, disp(sprintf('Warning: %d leading zero coefficients removed',zin)), end
zero = flipud(p(:));
zero = cumprod(zero==0);
zero = sum(zero);
n    = length(p);
fatt = p(1:n-zero);             % deflate the zero roots, which
r    = zeros(zero,1);           % initialize the vector r.
n0   = inf;
ciclo = 0;
first = 1;
if see, disp('---------------------START---------------------'), end
while length(fatt)>1
   ciclo = ciclo +1;
   [T,fatt1] = base( fatt, param, see );
   n1 = max(size(T));
   if n0>=n1
      r1 = eig(full(T.'));
      if first
         r = [r; sort(r1)]; first = 0;
         param(1) = param(1)^(5/6);    %  relax the tolerances after
         param(2) = param(2)^(8/9);    %  the first factor is detected.
         if see
            disp('tolerance off-diag. entry updated = '),disp(param(1))   
            disp('tolerance for factors updated     = '),disp(param(2))   
         end
      else
         r = update_r( r, r1 );
      end   
      if see, disp(sprintf('BTR2-cycle %d - found roots = ',ciclo)),disp(r1), 
              disp('...............................................'),pause, end
      n0 = n1;
      fatt = fatt1;
   else
      disp(sprintf('Warning: loss of accuracy after root %d!',length(r)))
      if n1==length(fatt)-1
         r = [r; sort( eig( full(T.') ) )];
      else
         r = [r; sort(eig( compan(fatt).' ))]; 
      end   
      break
   end
end
if see, disp('----------------------END----------------------'), end
return

function [T,fatt] = base( p, param, see )
%
%  
%
tol  = param(1);
n    = length(p) -1;
tolv = n*param(2);
told = param(3);
tols = param(4);
p   = p(:)/p(1);    % monic polynomial.
C   = [-p(2:n+1) zeros(n,1)];
d = scal( p, tols );     % diagonal matrix of scaling.
if see, disp('polynomial  = '),disp(p),
   disp('diagonal scaling = '),disp(d),pause, end
for i = 1:n-1
   C(i+1,1)  = C(i+1,1)/d(i+1);
   C(i,2)    = d(i+1)/d(i);
end   
p1  = polyder(p)/n;
U   = zeros(n); 
U1  = p1(:)./d;

T = spdiags(ones(n-1,1),-1,n,n); 
i = 1;
while i <= n
   v = multp( C, U1, told );
   if i>1  
      v = vdiff( v, U0, told );   
      v(i-1) = 0;    % set to 0
   end
   T(i,i) = v(i)*d(i);
   if i==n, break, end
   v = vdiff( v, T(i,i)*U1, told ); 
   v(i) = 0;         % set to 0
   T(i,i+1) = v(i+1)*d(i+1); 
   if see, disp(sprintf('diagonal and off-diagonal entries on row %d = ',i))
      disp(full(T(i,i:i+1))),pause, end
   if abs(T(i,i+1)) <= tol
      if see, disp('vector norm ='),disp(norm(v)),pause, end
      if norm(v) <= tolv
         if see, disp('FACTOR FOUND')
         else, disp('factor')
         end
         T = T(1:i,1:i); break
      else
         if see, disp('BREAKDOWN')
         else, disp('breakdown')
         end
         i1 = min(find(abs(v)>tol));
         k = i1-i-1;
         if see, disp('dimension of the breakdown block = '),disp(k),pause,end
         T(i,i+k+1) = v(i1)*d(i1);
         if see, disp(sprintf('(%d,%d)-entry on the upper breakdown row = ',i,i+k+1))
            disp(full(T(i,i+k+1))),pause,end
         v(1:i1-1)=0;                 % set to 0
         U = zeros(n,k+1);
         U(:,k+1) = v/T(i,i+k+1);
         for j = k:-1:1
            U(:,j) = multp( C, U(:,j+1), told ); 
         end
         v = vdiff( multp(C,U(:,1),told),U1,told );
         for j = 1:k+1
            T(i+1,i+j) = v(i+j)*d(i+j);
            v  = vdiff( v, T(i+1,i+j)*U(:,j), told );
         end   
         if see, disp(sprintf('entries %d...%d on breakdown row %d = ',i+1,i+k+1,i+1))
            disp(full(T(i+1,i+1:i+k+1))),pause, end
         U1 = U(:,k+1);
         if i1==n 
            i = i1; break 
         else
            T(i+1,i1+1) = v(i1+1)*d(i1+1);
            if see, disp(sprintf('entry %d on breakdown row %d = ',i1+1,i+1))
               disp(full(T(i+1,i1+1))),pause, end
            if norm(v) <= tolv
               if see, disp('FACTOR FOUND DURING BREAKDOWN')
               else, disp('breakdown factor')
               end
               T = T(1:i1,1:i1); i = i1; break
            else
               U0 = U1;
               U1 = v/T(i+1,i1+1);
            end
         end 
         i = i1+1;
      end
   else   
      U0 = U1; 
      U1 = v/T(i,i+1);
      i = i+1;
   end   
end   
fatt = d(i:n).*U1(i:n);   
return

function d = scal( p, tols )
%
%     Computes the diagonal scaling D for the companion matrix of the
%     polynomial P, which has already been deflated of the zero roots.
%
n    = length( p ) -1;
d    = abs( p(:)/p(1) ); 
d(1) = 1;    %   set to 1
if d(n+1)==0
   error('zero roots not removed')
end
i   = 2;
ind = find( d(i:n+1)<=tols*d(n+1) );
while length(ind)>0
   k = sum(cumprod(diff(ind)==1))+1;  % number of consecutive zero entries.
   ind = min(ind)+i-1;                    % index of the first 0.
   M = abs( d(ind+k)/d(ind-1) )^(1/(k+1));   % multiplicative factor.
   for i = ind:ind+k-1
      d(i) = d(i-1)*M;
   end   
   i = ind+k;
   ind = find( d(i:n+1)<=tols*d(n+1) );
end
d = 2.^round(log2(d(1:n)));           % closest power of 2.
return

function v = vdiff( x, y, told )
%
%  Difference between the vectors x and y, with check for zero entries.
%
v = x-y;
ind = find( abs(v)<=told*( abs(x)+abs(y) ) );
v(ind) = 0; 
return

function v = multp( C, u, told );
%
%  Product C*u, with C the companion matrix, memorized in compact form (2 vectors).
%
n = length(u);
v = C(:,2).*[u(2:n);0];
if u(1)~=0
   v = vdiff( C(:,1)*u(1), -v, told );
end   
return

function rr = update_r( r, r1 );
%
%  Updates the vector with the roots.
%
rr = r;
n  = length(rr);
for i = 1:length(r1)
   [mm,j] = min( abs(rr-r1(i)) );
   rr = [rr(1:j); r1(i); rr(j+1:n)];
   n  = n +1;
end
return
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% EOF %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%