So I have a very very basic technique for converting shifted Chebyshev polynomials to regular polynomials but it seems to have a large issue with numerical stability and I don't completely understand why (but I partially do). The problem even seems to occur in a less exaggerated way for unshifted Chebyshev polynomials of higher degree.
First I'll start out be defining what I mean by the shifted Chebyshev polynomials. Let $L^{[a,b]} : [-1, 1] \to [a, b]$ be the affine map from [-1, 1] such that $L^{[a, b]}(-1) = a$ and $L^{[a, b]}(1) = b$. We define the shifted Chebyshev polynomials $T^{[a, b]}_n(x) : [a, b] \to \mathbb{R}$ by the following recurrence relation.
$$T^{[a, b]}_0(x) = 1$$ $$T^{[a, b]}_1(x) = L^{[a,b]}(x)$$ $$T^{[a, b]}_n(x) = 2 \cdot L^{[a,b]}(x) \cdot T^{[a, b]}_{n-1}(x) - T^{[a, b]}_{n-2}(x)$$
The regular Chebyshev polynomials $T_n(x) = T^{[-1, 1]}_n$ are a special case. One can also prove that $T_n(L^{[a,b]}(x)) = T^{[a, b]}_n(x)$
Now lets say I have an approximation of $f : [a, b] \to \mathbb{R}$ by the truncated Chebyshev expansion $P(x) = \sum_{k=0}^n c_k T^{[a, b]}_k(x)$. I can evaluate $P(x)$ using Clenshaw's algorithm with good success and we can see that $P(x) \approx f(x)$ for $x \in [a, b]$
How ever if I convert $P(x)$ to be represented as a standard polynomial $\sum_{k=0} a_k x^k$ and I then evaluate this polynomial with Horner's method I hit big issues. The results are off by a good bit and often more naive methods than using Chebyshev approximation produce better results than this converted polynomial. I use the following technique to do the conversion.
Let $L^{[a, b]}(x) = Ax + B$ then let $C^{[a, b]}_n(k)$ be the $k$th coefficient of $T^{[a, b]}_n(x)$, the following (rather complicated) recurrence defines it (you can ignore all but the last equation here to understand the remainder of my question)
$$C^{[a, b]}_0(0) = 1$$ $$C^{[a, b]}_1(0) = B$$ $$C^{[a, b]}_1(1) = A$$ $$C^{[a, b]}_n(0) = 2 B \cdot C^{[a, b]}_{n-1}(0) - C^{[a, b]}_{n-2}(0)$$ $$C^{[a, b]}_n(k) = 2 A \cdot 2 C^{[a, b]}_{n-1}(k-1) + B \cdot C^{[a, b]}_{n-1}(k) - C^{[a, b]}_{n-2}(k) \,\,\, 0 < k < n$$ $$C^{[a, b]}_n(n) = 2 A \cdot 2 C^{[a, b]}_{n-1}(n-1)$$
After we have this $a_k = \sum_{j=0}^n c_j \cdot C^{[a, b]}_j(k)$ and it kinda seems like this works in practice but the values of $C^{[a, b]}_n(n) = 2^{n-1}A^n$ are very large which seems to be causing the issue of numerical instability. It doesn't really seem like the computation of $C^{[a, b]}_n(n)$ itself is numerically stable frankly.
Is there a way to adapt this to improve stability? Is there a better technique for converting from Chebyshev approximations to regular polynomials? Is it just the case that these coefficients get very large and thus nothing can be done about this issue?