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 simple-root polynomial u(x)=0, instead of the original multiple-root 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, re-constructing 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 built-in 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 * (100x-1)^10
p(x) = (x+1)^40 * (x-2)^30 * (x+3)^20 * (x-4)^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(pz-p)/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(n-nc,m-mc);
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.e-6)): max(find(abs(p3)>1.e-6)))];
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.e-3, break; end;
end;
g = [p1,zeros(1,nz)];
function p = polyget(A)
% *** A polynomial coefficient vector from sub-polynomial 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 -1-2i -1-2i -1-2i 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