Numerical decomposition of polynomials into second-order terms

Question

I am trying to numerically decompose real univariate polynomials into a product of second order polynomials. I have attempted to do this using Bairstow's method, and while this works fine for most even-order polynomials (odd-order polynomials obviously don't work as-is, but can easily be worked around), there are some where the method fails as the Jacobian becomes singular, e.g. $x^{16} + x^{12} + x^8 + x^4 + 1$ with the guess $x^2$. Which alternative methods work? As I actually want the second-order polynomials and not pairs of complex zeros (that I may first need to match up), I would strongly prefer methods specialized for real polynomials that directly extract the second-order terms. The polynomials I am looking at often have multiplicity $\ge 2$ of some of the roots. Parallel extraction of all terms would be nice to have, but not strictly necessary.