function [r,m,b,infout] = btr( p, infoin )
%
% BTR   Computes the roots and the respective multiplicities of a polynomial.
%       (tested with Matlab 5.3).
%
% FORMAT:
%
%        [R,M,B,INFOUT] = BTR(P,INFOIN)
%
% where
%
% - P      = coefficients of the polynomial.
% - INFOIN = optional flag to be set to
%             0 : default (M set only in this case);
%             1 : if it is known that all the roots are simple;
%             2 : if only the distinct roots are required.
% - R      = contains the computed roots.
% - M      = contains the multiplicities of the found roots, if required.
% - B      = contains the centers of probable clusters of roots.
% - INFOUT = error flag set to
%              0 : all the roots and multiplicities,  if required,  have been
%                  properly computed;
%              1 : all the roots and multiplicities,  if required,  have been
%                  properly computed;  possible centers of clusters have been
%                  detected too;
%              2 : the distinct roots are in B,  and they are reliable approx-
%                  imations; less accurate approximations of all roots are in
%                  R; the multiplicities have not been completely computed: M
%                  contains a lower bound for the  exact multiplicity  of the
%                  corresponding root in B;  in particular, the highest multi-
%                  plicities are not  accurate,  while the others are correct:
%                  to this end,  the  multiplicity  of the zero root,  if any,
%                  must be considered apart (it is always correct);
%              3 : the computed roots, in R,  may be not very accurate, since
%                  it is not possible to recognize  any multiplicity;  approx-
%                  imations to the probable distinct roots are in B.
%
% REMARKS: when INFOUT=2,3, R contains the output of the Matlab function ROOTS.
%
% AUTHOR: Luigi Brugnano
%         Dipartimento di Matematica "U.Dini"
%         Viale Morgagni 67/A
%         50134 Firenze (Italy)
%         E-mail: brugnano@math.unifi.it
%
% DATE: January 30, 1999.
%
% 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.
%

% Needed functions: SCAL,TRIS,CHKMP,UPDATE,MAKEPI.
%
format long e, flops(0)
r = []; m = []; b = []; infout = 0;
if nargout<4 | nargin<1
   error('Not enough I/O arguments');
elseif nargin==1
   infoin = 0;        % default value for INFOIN
elseif min(abs([0 1 2]-infoin))>0
   infoin = 0;        % default value for INFOIN
elseif infoin==1                      % for all simple roots the Matlab func-
   r = roots( p ); return             % tion ROOTS gives comparable accuracy
end
if p(1)==0
   i = sum(cumprod(p==0));     % check for the first nonzero coefficient of P
   p = p(i+1:length(p));
end
n = length( p );
zero = 0;                                      % check for the zero root
for i = n:-1:1
   if p(i)~=0, break, else, zero = zero + 1; end
end
if zero
   n = n - zero;
   p = p(1:n);                                 % deflate the zero root
end
if n<=1                                  % the polynomial is a constant
   reverse = 0; go_out = 1;
elseif n==2                              % the polynomial as degree 1
   r = -p(2)/p(1); if infoin==0, m = 1; end
   reverse = 0; go_out = 1;
elseif n>2*fix(n/2) & all(p(2:2:n)==0)   % the polynomial is biquadratic
   [r,m,b,infout] = btr( p(1:2:n), infoin );
   if length( m )>0, m = kron( m(:), [1;1] ); end
   if length( r )>0, r = sqrt( r(:) ); r = [r;-r]; end
   if length( b )>0, b = sqrt( b(:) ); b = [b;-b]; end
   reverse = 0; go_out = 1;
else
   if ( norm( p(1:n-1), inf )/abs( p(n) ) )<( 1e-2/n )
     reverse = 1; p = p(n:-1:1);         % the reverse polynomial is processed
   else
     reverse = 0;
   end
   go_out = 0;
end
p  = p/p(1);
p0 = p;
left  = (n-1)*(go_out==0);                   % number of roots to be found
first = 1; err = 0;                              % first=1 until p==p0
while left>0                             % main loop
   [fat,d,cp,p1] = scal( p, ~first );           % normalization of p
   [t,mP,p2,brk] = tris( cp, p1 );       % basic procedure
   if max(size(t))>1 & first                            % in the case of all
      ll = eig( t(1:2,1:2) );                           % simple roots, notri
      if norm(ll,inf)==0 & n<20                         % ==1 prevents to
         notri = 0;                                     % process the tridia-
      else                                              % diagonal matrix
         notri = ( norm( ll, -inf )/norm( ll, inf )< ( eps^(.4) ) )|(n>=20);
      end
   else, notri = 0;
   end
   if brk(1)==2              % last found multiple roots not reliable
      if first==0            % the distinct roots have already been found
         infout = 2;                    % lower bound for the multiplicities
         b = r;                         % available
      else                   % no root found before
         infout = 3;                    % no information about multiplicities
         b = eig( t )/fat;              % available
      end
      r = roots( p0 );
      left = 0;
   elseif infoin==2                     % only the distinct roots needed
      if max(size(t))<left | ( notri==0  & brk(1)==0 )
         r = eig( t )/fat;
      else
         r = roots( p );
      end
      left = 0;
   else                      % infoin==0
      if brk(1)<0                       % a breakdown occurred at step -brk(1)
         ii = max( size( t ) );
         nt = norm( t(1:-brk(1),1:-brk(1)), inf );
         et = norm( t(-brk(1),1-brk(1):ii), inf );
         if first & -brk(3)==1       % no roots found since now and only 1 brk
            simple = ( ii==left );
            ii1 = chkmp( mP(-brk(2):length(mP)), simple );
         else, ii1 = 0;
         end
         if ii1<0 & first
            ii  = ii1 + brk(2) + 1;   % -ii probable distinct centers due to
            err = 2;                  %  a not reliable deflation recognized
         elseif nt>20*et                  % probable centers recognized
            b = eig( t(1:-brk(1),1:-brk(1)) )/fat;
            infout = 1;
         elseif -brk(1)==1
            b = t(1,1)/fat; infout = 1;
         elseif -brk(1)<=floor( left/3 )
            b = eig( t(1:-brk(1),1:-brk(1)) )/fat;
            nb = length( b );
            sumr = -p(2)/left;
            if nb==2
               alpha1 = polyfit( real([b(1) sumr]), imag([b(1) sumr]), 1 );
               alpha2 = polyfit( real(b), imag(b), 1 );
               if real(sumr)>max( real(b) ) | real(sumr)<min( real(b) ) | ...
                  imag(sumr)>max( imag(b) ) | imag(sumr)<min( imag(b) ) | ...
                  norm(alpha1-alpha2,inf)>eps^(.35)*norm(alpha2,inf)
                 b = []; end
            elseif nb>2
               fi = atan2( imag(b-sumr), real(b-sumr) );         % sorted
               fi = fi + (fi<0)*(2*pi);                          % arguments
               fi = sort(fi(:));                                 % of b-sumr
               deltafi = [fi(2:nb); 2*pi+fi(1)] - fi;
               if any( deltafi>=pi ), b = []; end % sumr is outside the convex
            end                                   % hull of the b(i)s
            if length(b)>1, infout = 1; end
         end
      else
         simple = max( size( t ) )==left & first;
         if left>1, ii = chkmp( mP, simple );      % number of roots found
         else,      ii = 1;                        % in this iterate
         end
         ii = min( [ii max(size(t))] );
      end
      if ii<0 & first
         err = 2;                              % not reliable deflation
         r = eig( t(1:-ii,1:-ii) )/fat;        % recognized
      elseif ii<0
         err = 1;                       % multiplicities no more reliable
      elseif max(size(t))==left & first
         if notri | brk(1)<0, r = roots( p );
         else,                r = eig( t )/fat;
         end                            % all simple roots
         m = ones(left,1); left = 0;
      elseif ii>1
         if first==0 & ii>length( r )
            err = 1;                    % multiplicities no more reliable
         else
            ri = eig( t(1:ii,1:ii) )/fat;         % found roots
            [r,m,first,err] = update( r, m, ri, first );           % update m
            if err==0
               left = left - ii;
               if left>0 & ii<length(mP)-1     % deflated problem
                  p = makepi( cp, p1, ii, d );
               elseif left>0
                  p = d(ii:ii+left).*p2(ii:ii+left); p = p/p(1);
               end
            end
         end
      else
         ri = t(1,1)/fat;                      % found root
         for i = 1:left
             [r,m,first,err] = update( r, m, ri, first );          % update m
             if err, break, end
         end
         left = 0;
      end
      if err
         if err==1, infout = 2;    % the multiplicities are no more reliable
         else,      infout = 3;    % no information on multiplicities
         end
         b = r;
         r = roots( p0 );
         left = 0;
      end
   end
end
if reverse, r = ones(length(r),1)./r; end    % roots of the reverse poly
if zero                                          % the zero root is inserted
   m1 = zero;                            % multiplicity of the zero root
   if infout==0 | infout==1
      r = [r(:); 0];  if infoin~=2, m = [m(:); m1]; end
      if infout==1, b = [b(:); 0]; end
   elseif infout==2
      r = [r(:); zeros(m1,1)];
      b = [b(:); 0];
      m = [m(:); m1];
   else % infout==3
      r = [r(:); zeros(m1,1)]; b = [b(:); 0];
   end
end
return



function [fat,d,cp,p1] = scal( p, notfirst )
%
%     Computes the diagonal scaling D for the companion matrix of the
%     polynomial P, which has already been deflated of the zero roots.
%     In output, CP contains the normalized coefficients and the ext-
%     radiagonal of the scaled companion matrix; P1 is the normalized
%     derivative. FAT is the scaling factor which normalizes CP.
%
n = length( p ) - 1;
tol = notfirst*sqrt( abs( eps*p(n+1) ) ); % in the first scaling, tol=0.
d = zeros(n,1); p = p(:)/p(1);
i0 = 1; d(1) = 1; last = 1; i1 = 1;
while i1<=n-last
   mm = min( find( abs( p(i0+1:n+1) )>tol*d(i0) ) );
   if isempty(mm), mm = 0; end, i1 = i0 + mm;
   if i1-i0==1
      d(i1) = abs( p(i1) ); i0 = i1; last = 1;
   elseif i1>i0
      m = abs( p(i1)/p(i0) )^( 1/(i1-i0) );
      for i = i0+1:i1-1; d(i) = d(i-1)*m; end
      i0 = i1 - 1; last = 0;
   else       % there are 2 scales and the problem is ill-conditioned
      for i = i0+1:n, d(i) = d(i-1); end
      i1 = n + 1 + last;
   end
end
cp = [-p(2:n+1)./d [d(2:n)./d(1:n-1); 0] ];
p1 = polyder( p ).'/( p(1)*n );
p1 = p1./d;
i1 = find( cp(:,1)~=0 );
if n==1
   ma = norm(cp(i1,1),inf);
   mi = norm(cp(i1,1),-inf);
else
   ma = max( [norm(cp(i1,1),inf) norm(cp(1:n-1,2),inf)] );
   mi = max( [norm(cp(i1,1),-inf) norm(cp(1:n-1,2),-inf)] );
end
fat = 2^(-round( log( sqrt( ma*mi ) )/log( 2 ) ));
cp = cp*fat;
return


function [t,mP,p2,brk] = tris( A, p2 )
%
%     Basic procedure of tridiagonalization of the companion matrix A stored
%     in compact form. P2 is the first column of the transformation matrix,
%     In output T is the tridiagonal matrix, MP is a vector with the
%     maxima on each column of the transformation matrix, and BRK is a
%     flag set to:
%        0 - if no trouble occurs,
%      k<0 - if a breakdown occurs,
%        2 - if an unreliable deflation has been detected;
%     P2 is the last nonzero column.
%
n = length( p2 );
beta = 1; p1 = zeros(n,1); t = zeros(n);
mP(1) = norm(p2,inf);
p2 = p2(:);
v = A(:,1)*p2(1) + A(:,2).*[p2(2:n);0];
if n==2, tol = eps;
else,    tol = n*eps^(.75);
end
ztol = sqrt( eps );
brk  = 0; esci = 0;
i = 1; firstbrk = []; ii = 0;
while i<=n-1
    flag = 1;  p0 = p1; p1 = p2;
    alfa = ( v(i) - beta*p0(i) )/p1(i);
    t(i,i) = alfa;
    p2 = v - beta*p0 - alfa*p1;
    p2(1:i) = zeros(i,1);
    mP(i+1) = norm( p2(i+1:n), inf );
    gamma = sqrt( abs( A(i,2)*p2(i+1)/p1(i) ) );
    if mP(i+1)==0 | ( gamma==0 & mP(i+1)<i*tol )
       esci = 1;                      % deflation
    elseif abs( p2(i+1) )<i*tol
       esci = 1;
       ratio = max([abs(p2(i+1)) eps]) / mP(i+1);
       e = [ratio<1e-3  ratio<6e-6 ...
            gamma > n*1e3*sqrt( abs(p2(i+1)) ) ...
            ( mP(i+1) > 1.234e-4 * gamma * mP(i) & gamma>0 ) |...
            ( mP(i+1) > i*tol & gamma==0 )];
       e1 = ( mP(i+1)<1e-5*gamma*mP(i) ) & gamma>0;
       se = sum( e );
       if length(firstbrk)>0 & e1     % deflation after a breakdown
       elseif e(2)
          esci = 0; brk = 1;          % breakdown
       elseif p2(i+1)~=0 & ( ( se & e(1)==0 ) | ( e(1) & e(4) ) )
          esci=0;                     % continue: the vector p2 is only small
       elseif ( e(1) & e(3) ) | ( e(1) & e(4)==0 ) | ( e(1)==0 & e(4) )
          brk = 2;                    % deflation not reliable
       end
    end
    if brk==1, brk = 0;            % breakdown
       for ii = i+1:n              % computes the first nonzero entry of p2
          if abs( p2(ii) )>ztol*max([1 mP(i+1)]), break, end
       end
       if length( firstbrk )==0               % firstbrk(1)=start first brk
          firstbrk = [i ii 1];                % firstbrk(2)=end last brk
       else                                   % firstbrk(3)=number brks
          firstbrk(2) = ii;
          firstbrk(3) = firstbrk(3) + 1;
       end
       p2(1:ii-1) = zeros(ii-1,1);
       k  = ii - i - 1;      % number of missed columns due to the breakdown
       gamma = sqrt( abs( p2(ii)*prod(A(i:ii-1,2))/p1(i) ) );
       t(i,ii) = gamma;
       p0 = p1; p1 = p2/gamma;
       mP(ii) = mP(i+1)/gamma;
       ausP = zeros(n,k);
       ausP(:,k) = A(:,2).*[p1(2:n);0];
       for j = k-1:-1:1
           ausP(:,j) = A(:,2).*[ausP(2:n,j+1);0];
       end
       p2 = A(:,2).*[ausP(2:n,1);0];
       t(i+1,i) = p2(i)/p0(i);
       p2 = p2 - p0*t(i+1,i);
       for j = 1:k
          t(i+1,i+j) = p2(i+j)/ausP(i+j,j);  t(i+1+j,i+j) = 1;
          p2 = p2 - ausP(:,j)*t(i+1,i+j);
          mP(i+j) = norm( ausP(i+j:n,j), inf );
       end
       t(i+1,i+k+1) = p2(i+k+1)/p1(i+k+1);
       p2 = p2 - p1*t(i+1,i+k+1);
       p2(1:ii) = zeros(ii,1);
       if norm( p2, inf )<10*ii*tol
          i = ii; mP(i+1) = 0; esci = 1;
       else
          flag = 0;
          i = ii + 1;
          if i<=n,
             v = A(:,2).*[p2(2:n);0];
             beta = v(ii)/p1(ii);
             gamma = sqrt( abs( beta ) );
             beta = beta/gamma;
             v(ii-1:n) = v(ii-1:n)/gamma;
             p2(ii:n) = p2(ii:n)/gamma;
             t(i,i-1) = beta;
             t(ii-k,ii+1) = gamma;
             mP(ii+1) = norm( p2(ii+1:n), inf );
          end
       end
    end
    if esci
       t = t(1:i,1:i);
       mP = mP(1:i+1);
       if brk==0, p2 = p1; end
       if length( firstbrk )>0, brk = -firstbrk; end
       return
    elseif flag & i<n
       t(i,i+1) = gamma;
       p2(i+1:n) = p2(i+1:n)/gamma;
       mP(i+1) = mP(i+1)/gamma;
       v = A(:,2).*[p2(2:n);0];
       beta = v(i)/p1(i);
       i = i + 1;
       t(i,i-1) = beta;
    end
end
if ii<n
   alfa = ( v(n) - beta*p1(n) )/p2(n);
   t(n,n) = alfa;
end
if brk==0, p2 = p1; end
if length( firstbrk )>0, brk = -firstbrk; end
return


function ii = chkmp( mP, simple )
%
%     Recover an eventual lost deflation.
%
n = length( mP ); c = mP;
if c(n)==0, c(n) = eps*min( c(1:n-1) ); end    % only c(n) may be 0
c = log( c );
c = c(1:n-1) - c(2:n);
[cmax,ii] = max( c );
if cmax>9 & mP(ii+1)<1e-3               % ii is the probable deflation
elseif cmax>7-simple, ii = -ii;           % the deflation is not reliable
else,                 ii = n;           % no deflation recognized
end
return


function [r,m,first,err] = update( r, m, ri, first )
%
%     Updates the multiplicities M of the roots R, if FIRST==0, with the
%     new found roots RI. If FIRST==1, R is set to RI and M is set to 1.
%     When the multiplicities can not be reliably recognized,  M  is not
%     updated, and the error flag ERR is set to 1.
%
err = 0;
if first
   r = ri(:); m = ones( length(r), 1 ); first = 0;
else
   rold = r; mold = m; f = eps^(.38);
   for i = 1:length( ri )
      [d,ind] = min( abs( r - ri(i) ) );
      if d > abs(r(ind))*f
         err = 1; r = rold; m = mold; break
      else
         m(ind) = m(ind) + 1;
      end
   end
end
return


function pi = makepi( A, p2, ii, d )
%
%     Computes the coefficients of the deflated problem.
%
n = length( p2 );
p2 = p2(:); p1 = zeros(n,1);
v = A(:,1)*p2(1) + A(:,2).*[p2(2:n);0];
beta = 0;
for i = 1:ii-1
    p0 = p1; p1 = p2;
    alfa = ( v(i) - beta*p0(i) )/p1(i);
    p2 = v - beta*p0 - alfa*p1;  p2(1:i) = zeros(i,1);
    gamma  = sqrt( abs( A(i,2)*p2(i+1)/p1(i) ) );
    p2(i+1:n) = p2(i+1:n)/gamma;
    v = A(:,2).*[p2(2:n);0];
    beta = v(i)/p1(i);
end
p2 = d.*p2;
pi = p2(ii:n)/p2(ii);
return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% END OF FILE %%%%%%%%%%%%%%%%%%%%%%%%%%%%%