The computation of the roots of polynomials expressed in the Chebyshev basis has a lot of applications, for instance, it is useful in the computation of real roots of smooth functions.

We present an algorithm for the rootfinding of Chebyshev polynomials based on an improvement of the QR iteration presented in [Eidelman, Y., Gemignani, L., and Gohberg, I., Numer. Algorithms , 47.3 (2008): pp. 253-273]. We introduce an aggressive early deflation strategy, and we show that the rank-structure allows to parallelize the algorithm avoiding data dependencies which would be present in the unstructured QR. The method exploits the particular structure of the colleague linearization to achieve quadratic complexity and linear storage requirements. The (unbalanced) QR iteration used for Chebyshev rootfinding does not guarantee backward stability on the polynomial coefficients, unless the vector of coefficients satisfy ||p|| ~ 1, an hypothesis which is almost never verified for polynomials approximating smooth functions. Even though the presented method is mathematically equivalent to the latter algorithm, we show that exploiting the rank structure allows to guarantee a small backward error on the polynomial, up to an explicitly computable amplification factor ɣ(p), which depends on the polynomial under consideration. We show that this parameter is almost always of moderate size, making the method accurate on several numerical tests, in contrast with what happens in the unstructured unbalanced QR. We also discuss the connection between the size of this amplification factor and the existence of a good balancing. This provides some insight on why the accuracy of our method is often very close to balanced QR iteration.