The polynomial
Can be regarded as the characteristic polynomial of the matrix:
And if we compute the eigenvalues of this matrix, we get the roots of the polynomial.
I do this by the QR decomposition (see Eigenvalues QR decomposition) and it works fine. The algorithm converges quite reliable.
There is only one limitation remaining: B may not be 0. otherwise the algorithm, as it is, would not converge. But this limitation can be avoided by modifying the “shift” procedure of the QR algorithm a little. Usually for the shift value (phi in the algorithm) the last element of the main diagonal is used and in the first loop this is –B. If B is 0 this would not change anything. But in this case we can set phi to a random value and use this. I set phi = 0.2 if the last element of the main diagonal is 0 for each loop. That works fine.
Whit this the QR algorithm works really fine to compute all roots and it's advantage over a “Newton Maehly” (see Real roots of polynomials by Newton/Maehly) or a “Bairstow” (see Complex roots of polynomials by Bairstow) algorithm is that there is no deflation required and it can compute all roots, real and complex, at once. That makes it quite cool
