function r = mroots(p,flag)
%MROOTS   Computes polynomial roots and approximates the true value and
%         multiplicity of multiple roots.
%
%         R = MROOTS(P)  returns the roots of the polynomial P.

%      Author: P. Gahinet, 5-1-96
%      Copyright (c) 1986-97 by The MathWorks, Inc.
%      $Revision: 1.3 $


% Compute roots 
if nargin==1,
   r = roots(p);
   realflag = isreal(p);
else
   % Syntax MROOTS(R,'roots') where roots R already computed
   r = p(:);
   realflag = isequal(sort(r(imag(r)>0)),sort(conj(r(imag(r)<0))));
end

% Check for clustering around zero
tolabs = eps^0.5;
absr = abs(r);
izero = (absr<tolabs);
rzero = zeros(sum(izero),1);
r = r(~izero,1);
absr = absr(~izero,1);

% Compute relative separation of roots (SEPMAT)
lr = length(r);
repr = r(:,ones(1,lr));
repabs = absr(:,ones(1,lr));
sepmat = eye(lr) + (abs(repr-repr.')./min(repabs,repabs'));

% Sort the roots such that r(k) is the closest roots
% (in the SEPMAT measure) to r(1),r(2),...,r(k-1)
ilist = 1:min(1,lr);
irem = min(1,lr)+1:lr;
gaps = zeros(1,lr);      % stores d(r(k),{r(1),r(2),...,r(k-1)})
for k=2:lr,
   [gaps(k),imin] = min(min(sepmat(ilist,irem),[],1));
   ilist = [ilist irem(imin)];
   irem(imin) = [];
end
r = r(ilist,1);

% Recursively check & split the cluster R
if length(r)>1,
   r = CheckCluster(r,gaps);
end

% Put result together
r = [rzero ; r];
if realflag,
   r(imag(r)<0) = conj(r(imag(r)>0));
end

% end mroots.m

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function rc = CheckCluster(rc,gaps)
%CLUSTER   Checks if a given root cluster RC can be identified with 
%          a multiple root.  If not, the cluster is split into 
%          subclusters and CHECKCLUSTER is invoked on each subcluster.


% Ignore first entry of GAPS vector
gaps(1) = 0;  

% Get multiplicity and cluster center
m = length(rc); 
rmean = sum(rc)/m;
if isequal(sort(rc(imag(rc)>0)),sort(conj(rc(imag(rc)<0))))
   % Make sure RMEAN is real
   rmean = real(rmean);
end

% Tolerances
rtol = eps^(0.7/m);  
ptol = 10*sqrt(eps);       % max. rel. error on coefficients of POLY(rc)
maxgap = min(0.1,2*rtol);  % max. gap between roots

% Test if cluster can be assimilated to root RMEAN of 
% multiplicity M.  All roots should lie in disk of center RMEAN 
% and radius RTOL * ABS(RMEAN), and the worst gap between 
% the coefficient of POLY(RC) and POLY(RMEAN*ONES(1,M))
% should not exceed PTOL
if any(gaps>maxgap),
   % Some gaps exceed threshold: split cluster
   splitgap = maxgap;
elseif max(abs(rc-rmean)) > rtol*abs(rmean),
   % Some roots lie outside reference disk
   splitgap = max(gaps)/2;
else
   % Perform polynomial test
   pmult = poly(rmean*ones(1,m));
   prc = poly(rc);
   if max(abs(prc-pmult)./abs(pmult)) > ptol,
      % Split cluster
      splitgap = max(gaps)/2;
   else
      % Cluster can be assimilated to multiple root RMEAN
      rc = rmean * ones(m,1);
      return
   end
end

% Split cluster if fusion into multiple root RMEAN failed
igap = find(gaps>splitgap);
isub = [1 , igap ; igap-1 m];  % all subclusters
isub = isub(:,diff(isub)>0);   % discard subcluster with single roots
for ix=isub,
   rc(ix(1):ix(2)) = CheckCluster(rc(ix(1):ix(2)),gaps(ix(1):ix(2)));
end

   
% end CheckCluster


