A polynomial is factored into lower-degree distict-root polynomials with natual-order-integer powers. Roots with multiplicites may then be found with easy.
In process, the crucial concern is computating polynomial GCD. The classical Euclidean GCD algorithm was applied in the previous version, however, it was unstable numerically.
In the presented version, the process "Monic polynomial subtraction" will replace "Longhand polynomial division" of the Euclidean algorithm for computing polynomial GCD.
This novel process is very simple, fast, accurate, stable, and requires minimal arithematic computations.

Reference:
F C Chang, "Factorization of Multiple-Root Polynomials"
F C Chang, "GCD of two polynomials by Monic polynomial subtraction"

Amazingly, the simple routine gives expected results for the test polynomials of very high degree, such as
p(x) = (x + 1)^1000
p(x) = (x^4 - 1)^500
p(x) = (x - 123456789)^30
p(x) = (1234x + 56789)^60
p(x) = (x^4 -2x^3 +3x^2 -4x +5)^150

Example run in MATLAB
>>
>> % Create a test polynomial:
>> p=poly([ 1 1 1 1 1 -1 -1 -1 -1+2i -1+2i -1+2i ...
>> -1-2i -1-2i -1-2i 2 2 3 3 -3])
p =
1 -3 -8 -4 76
284 -536 -808 -2474 7486
5896 -3872 -27284 -15812 77848
2184 -82319 24045 28800 -13500
>> % Get lower-degree distinct-root polynomials
>> W = PolyFct(p); celldisp(W)
W{1} = 1 3
W{2} = 1 -5 6
W{3} = 1 3 7 5
W{4} = 1
W{5} = 1 -1
>> % End of Example
>>