M_polyroots

function M_polyroots

    disp('  '),
    disp('**********************************************************************'),
    disp('The routine M_polyroots is to compute all the roots and multiplicities'),
    disp('  of any given polynomials.    '),
    disp('  developed by  F C Chang    Updated  01/07/11  ')
    disp('Ref: IEEE Antennas & Propagation Magazine, 51-6, p.151-155, Dec 2009  '),
    disp('  '),
    disp('**********************************************************************'),
    disp('Suppose a given polynomial p(x) is expressed as'),
    disp('  p(x) = x^9 +7x^8 +12x^7 -12x^6 -42x^5 -6x^4 +44x^3 +20x^2 -15x -9 '),
    disp('       = (x +3)^2 *(x +1)^4 *(x -1)^3 '),
    disp('       = (x^2 +2x -3)^2 *(x +1)^3 *(x^2 -1) '),
    disp('it does give '),
    disp('     c = [ +1 +7 +12 -12 -42 -6 +44 +20 -15 -9 ] '),
    disp('     r = [ -3 -3  -1 -1 -1 -1  +1 +1 +1 ] '),
    disp('     A = [  2  -3 +2 +1;  3  +1 +1 +0;  1  -1 +0 +1 ] '),
    disp('Then the polynomial coefficient vector p  may be created by either  '),
    disp('    (1)  p = c            '),
    disp('    (2)  p = poly(r)      '),
    disp('    (3)  p = polyget(A)   '),
    disp('The desired roots/multiplicities Z  are therefore computed from p  by '),
    disp('         Z = polyroots(p) '),
    disp('  '),
    disp('**********************************************************************'),
    disp('**********************************************************************'),
    format compact,
    format long g,
    warning off MATLAB:colon:operandsNotRealScalar
   
  for k = 1:100,
    disp('/  '),
    disp(' To create a coefficient vector by selecting the one from  '),
    disp('  (1) p = c, (2) p = poly(r), (3) p = polyget(A),  or (0) exit '),
           L = input(' Select 1, 2, 3, or 0  --->   ');
    if     L == 0,  disp('  Exit ');    break,   end;
    if     L == 2,  r = input(' Enter r =  ');  p = poly(r);  
    elseif L == 3,  A = input(' Enter A =  ');  p = polyget(A);
    else,           c = input(' Enter c =  ');  p = c;   
    end;
    if length(p) < 50,  p,   else,  disp('p = (omit) '),   end;
           Z = polyroots(p);
          
       fprintf('\n ');
    if norm(imag(Z(:,1))) == 0, 
       fprintf('          computed roots          multiplicities\n');
       fprintf('%25.15f              %3g \n', Z');
    else
       fprintf('                 computed roots                           multiplicities\n');
       fprintf('%22.15f   %+22.15f i           %6g \n', [real(Z(:,1)),imag(Z(:,1)),Z(:,2)]');
    end;
       disp('-------------------------------------------------------------------'),
           W{1,k} = p;  W{2,k} = Z;  
  end;
       YN = input(' Print out computed results now ?    1 or 0  --->    ');
    if YN == 1,
       disp('  ');
       disp('      *****   Poly coefficient vector :   p = W{1,*}  *****');
       disp('      *****   Roots /  Multiplicities :   Z = W{2,*}  *****');
       disp('  ');
       celldisp(W);   pause,   
    end;
       input('  Type "0" to terminate the program.  ');


         
function p = polyget(A)
% *** Polynomial coefficients expanded from powered factors ***
%     F C Chang  Updated  12/22/2010
    p = 1;
for i = 1:length(A(:,1)),
    q = 1;
for j = 1:A(i,1),
    q = conv(q,A(i,max(find(A(i,:))):-1:2));
end;
    p = conv(p,q);
end;
         


function Z = polyroots(p)
% *** Solve multiple-root polymonials ***
%     F C Chang  Updated  12/22/2010
    mz = length(p)-max(find(p));
    p0 = p(min(find(p)):max(find(p)));
 if length(p0) < 2, Z = [0,mz]; return, end;
    sr = abs(p0(end)/p0(1));
 if sr < 1,   p0 = p0(end:-1:1);     end;
    q0 = polyder(p0);
    g1 = p0/p0(1); 
    g2 = q0(1:max(find(q0)))/q0(1);
for k = 1:2*length(p0),
    z12 = zeros(1,length(g1)-length(g2));
    z21 = zeros(1,length(g2)-length(g1));
    g3 = [g2,z12]-[g1,z21];
    g3 = g3(min(find(abs(g3)>1.e-8)):max(find(abs(g3)>1.e-8))); 
 if norm(g3,inf)/norm(g2,inf) < 1.e-3,
                             break;  end;
 if length(g1) >= length(g2),  
    g1 = g2;                         end;
    g2 = g3/g3(1);
end;
    g0 = g1;
    u0 = deconv(p0,g0);
    v0 = deconv(q0,g0);
    w0 = polyder(u0);
    z0 = roots(u0);
    m0 = polyval(v0,z0)./polyval(w0,z0);
 if sr < 1,  z0 = z0.^-1;            end;
    Z = [z0,round(abs(m0))];
 if mz > 0,  Z = [Z; 0,mz];          end;