To find poles/residues of the rational function, instead of roots/multiplicities of the polynomial
A given polynomial p(x) is transformed into a rational function r(x). The poles and residues of the derived rational function are found to be equivalent to the roots and multiplicities of the original polynomial.
p(x) = Given polynomial
= PROD[k=1:K]{(x  z_k)^m_k}
d(x) = (d/dx)p(x)
g(x) = GCD(p(x),d(x))
u(x) = p(x)/g(x)
w(x) = (d/dx)u(x)
v(x) = d(x)/g(x)
r(x) = v(x)/u(x)
= SUM[k=1:K]{m_k/(x  z_k)}
Thus, the roots z_k are computed from solving the simpleroot polynomial u(x)=0, instead of the original multipleroot polynomial p(x)=0; and the multiplicities m_k are determined as the partial fraction expansion coefficients of the derived rational function r(x)=v(x)/u(x),
z_k = Roots(u(x)), k=1,K
m_k = v(z_k)/w(z_k), k=1,K
In addition, reconstructing a polynomial pz(x) from the computed z_k and m_k, the overall deviation error of the original polynomial p(x) is calculated,
er = Norm(pz  p)/Norm(p)
The polynomial GCD is calculated from "Monic polynomial subtraction" derived from the longhand polynomial division in classical Euclidean GCD algorithm. It requirs only simple algebric operations without any high mathematics.
The source code contains total of only 43 lines, using merely basic builtin MATLAB functions, and applying only existing double precision. Amazingly, it gives the expected results of test polynomials of very high degree , such as
p(x) = (x  123456789)^30
p(x) = (x + 100)^20 * (100x1)^10
p(x) = (x+1)^40 * (x2)^30 * (x+3)^20 * (x4)^10
p(x) = (x + 1)^1000
______________________________________________________
The code is list here for reader's convenince. (only 43 lines)
function [zm,er] = polyroots(p)
% *** A polymonial with multiple roots ***
% Solved via partial fraction expansion
d = polyder(p);
g = polygcd(p,d);
u = deconv(p,g);
v = deconv(d,g);
w = polyder(u);
z = roots(u);
m = round(abs(polyval(v,z)./polyval(w,z)));
zm = [z,m]; % p,d,g,u,v,w,z,m,zm
pz = polyget([m,z,ones(length(z),1)])*p(1);
er = norm(pzp)/norm(p); % pz,er
function g = polygcd(p,q)
% *** GCD of a pair of polynomials ***
% by "Monic polynomial subtraction"
n = length(p)1; nc = max(find(p))1;
m = length(q)1; mc = max(find(q))1;
nz = min(nnc,mmc);
if nc*mc == 0, g = [1,zeros(1,nz)]; return, end;
p2 = [p(1:nc+1)];
p3 = [q(1:mc+1)];
for k = 1:nc+nc,
p3 = [p3(min(find(abs(p3)>1.e6)): max(find(abs(p3)>1.e6)))];
p1 = [p2/p2(1)]; % k,p1,
p2 = [p3/p3(1)];
p3 = [p2,zeros(1,length(p1)length(p2))][p1,zeros(1,length(p2)length(p1))];
if norm(p3)/norm(p2) < 1.e3, break; end;
end;
g = [p1,zeros(1,nz)];
function p = polyget(A)
% *** A polynomial coefficient vector from subpolynomial factors ***
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;
_____________________________________________________
Typical Numerical Example:
>> % Contruct a test polynomial:
>> p = poly([ 1 1 1 1 1 1 1 1 1 1 1+2i 1+2i 1+2i 12i 12i 12i 2 2 3 3 +i +i +i i i i 3 0 0 0 0 0 ])
p =
1 5 2 6 76 140 802 954 4251 13663 18740 28472 53504 45776 5212 77580 185243 220631 104794 52458 193356 248612 146266 9202 65791 87555 55800 13500 0 0 0 0 0
>> % Roots and multiplicities for the polynomial are computed.
>> zm = polyroots(p)
zm =
3.0000  : 2.0000
3.0000  : 1.0000
1.0000 + 2.0000i : 3.0000
1.0000  2.0000i : 3.0000
2.0000 : 2.0000
1.0000 : 3.0000
0.0000 + 1.0000i : 3.0000
0.0000  1.0000i : 3.0000
1.0000 : 7.0000
0.0000 : 5.0000
1.11  Correct typo in mfile 

1.10  Update the mfile, to include the overall deviation error of the original polynomial. 

1.7  Update the Polynomial GCD routine. 

1.4  Update the mfile to include both monichead and monictail for polynomial GCD computation. 

1.2  Update the source code. Add two numerical examples, including printout the complete PRS and related polynomials. 

1.1  Add numerical example to "Description". 
Inspired by: More Flexible Sorting and Multiplicity of Roots of a Polynomial, Symbolic Polynomial Manipulation, Variable Precision Integer Arithmetic
Inspired: Polynomials with multiple roots solved