Numerical evaluation of polynomials in Chebyshev basis I have high order (15 and higher) polynomials defined in Chebyshev basis and need to evaluate them (for plotting) on some intervals inside the canonical interval $[1,\,-1]$. A good accuracy near 1 and -1, where Chebyshev polynomials change rapidly, is also required. 
I know that there exists Clenshaw algorithm for evaluation of Chebyshev polynomials, which is somewhat similar to the Horner scheme.
And this is all that I know...
I also saw the question about similar problematics. The proposed solution is to use higher precision. 
I wonder if there are some other methods for accurate evaluation of polynomials (particularly in Chebyshev basis) on the given intervals that don't heavily rely on extended precision calculations? Is it possible at least to improve the problem as much as possible, so that extended precision is needed only from really high polynomial orders? 
I also thought about interpolation of a high order polynomial by lower order ones on the intervals of interest. However, I don't know if there's any systematic procedure for this approach.

Evaluating polynomials of arbitrarily large degree in a Chebyshev basis is practical, and provably numerically stable, using a barycentric interpolation formula. In this case, extended precision isn't needed, even for order 1,000,000 polynomials. See the first section of this paper and the references, or here (Myth #2) for more details. I'll summarize briefly. Let's say you have a polynomial $f$ in a Chebyshev basis, and you know its values at the Chebyshev nodes
$$
f_j = f(x_j)
$$
$$
x_j = \cos\left(\frac{\pi j}{N}\right), \;\; 0\leq j\leq N
$$
Then for any $x$ which isn't one of the Chebyshev nodes $x_j$, we have
$$
f(x) = \frac{\displaystyle \sum_{j=0}^N \frac{w_j}{x-x_j}f_j}{\displaystyle \sum_{j=0}^N \frac{w_j}{x-x_j}},
$$
where 
$$
w_j = \left\{
\begin{array}{cc}
(-1)^j/2, & j=0\text{ or }j=N,\\
(-1)^j, & \text{otherwise}
\end{array}
\right.
$$
I believe a naive implementation of the above for various values of $x$ provides a stable algorithm that does not suffer the numerical difficulties encountered when trying to sum up the Chebyshev polynomials directly. The key is to work with the representation of the function by its values, not by its coefficients in a Chebyshev basis.
