Obtaining Coefficients of Powers of polynomial (e.g., $P(x)^N$) for large N, becomes Numerically Unstable
I have a polynomial $P(x)$ where $-1\leq P(x) \leq 1$ for $-1 \leq x \leq 1$ and $-1\leq {\it coefficients} \leq 1$.
I want to find the coefficients of $P(x)^N$ where $N$ is an integer. In Matlab, I use coeffs(P(x)^N) to obtain the coefficients of $P(x)^N$. To check the results, I plot P(x)^N and also constructed polynomial using the obtained coefficients. However, as $N$ increases, the second polynomial becomes unstable.
Example:
P= 31/40 - 7/8*x^2 - 9/10*x^4
[Co,Vo]=coeffs(P^100);
P2=sum((Co).*Vo);
P2 doesn't match the P^100 and is numerically unstable (e.g. P2(x) $\rightarrow \inf$ for some $x \in [-1 1]$)
Does this related to "finite precision computation" of coefficients?
You can find this effect already with the simpler polynomial $$ (1-x)^N=\sum_{k=0}^N\binom{N}{k}(-x)^k $$ For $x\approx1$ the left side is very close to $0$. Evaluating the right side involves very large coefficients that have to cancel. The largest is the central one, $$ \binom{N}{\lfloor N/2\rfloor}\approx\frac{\sqrt{2\pi N}(N/e)^N}{\pi N(N/(2e))^N}=\frac{2^N}{\sqrt{\pi N/2}}. $$ As $2^{100}\sim 10^{30}$ one can expect a floating point error of size coefficient times machine epsilon, $\binom{100}{50}\epsilon\sim 2^{100-8/2-52}=2^{44}\sim 1.7\cdot 10^{13}$. Doing the actual numerical computation confirms this scale of the error.