Code covered by the BSD License

# Solving multiple-root polynomials

### Feng Cheng Chang (view profile)

03 Aug 2010 (Updated )

Find roots and multiplicities of given polynomials using this short compact routine.

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;

```