How to find the coefficients of an iteration of a polynomial function? Let $f(x)$ be a polynomial of degree $N$ with real coefficients $c_n$:

$$ f(x) = \sum_{n=0}^{N} x^n c_n$$

Further, let $f^{(1)}(x) =f(x)$, $f^{(2)}(x) =f(f(x))$ and so on ...

It is obvious that the degree of $f^{(i)}(x)$ is $N^i$.

Here is my question:
How do you calculate the coefficients of $f^{(i)}(x)$?
I dont know how to tackle this.

If for example $i=2$ then 
$$f(f(x)) =  \sum_{n=0}^{N} c_n ( \sum_{k=0}^{N} x^k c_k)^n$$

But i dont get anywhere with this. Any help is appreciated.
Edit: See below for a much better way to derive the main result.
Represent the coefficients of a polynomial $p$ as a map $c_p$ from integers to real numbers, where $c_p(j)$ is the $j^{th}$ coefficient of $p$. Equivalently, $c_p$ is just the vector of coefficients of $p$.  Multiplying two polynomials $p$ and $q$ corresponds to taking a discrete convolution of the corresponding coefficient vectors:
$$
c_{pq}(i) = \sum_j c_p(i-j) c_q(j) = (c_p*c_q)(i)
$$
By the convolution theorem, the Fourier transform of a convolution is a product of the individual Fourier transforms.  Hence, the Fourier transform of the coefficient vector of $p^{k}$ (Caution! This is a regular power of $p$, not the same as $p^{(k)}$) is simply $\mathcal{F}[c_p^{*k}]=\mathcal{F}[c_p]^k$, where $\mathcal{F}$ is a discrete Fourier transform and $c_p^{*k}$ means $c_p$ convolved with itself $k$ times.
Now from the relation $p^{(k)}(x)=\sum_j c_p(j) \left(p^{(k-1)}(x)\right)^j$ we find that the coefficient vector $c_{p^{(k)}}$ satisfies $$c_{p^{(k)}}=\sum_j c_p(j) c_{p^{(k-1)}}^{*j}$$ hence the Fourier transform satisfies $$\mathcal{F}[c_{p^{(k)}}]=\sum_j c_p(j) \mathcal{F}[c_{p^{(k-1)}}^{*j}]=\sum_j c_p(j) \mathcal{F}[c_{p^{(k-1)}}]^j=p\left(\mathcal{F}[c_{p^{(k-1)}}]\right)$$
and inductively
$$\mathcal{F}[c_{p^{(k)}}]=p^{(k-1)}\left(\mathcal{F}[c_{p}]\right)$$
Hence $$c_{p^{(k)}}=\mathcal{F}^{-1}\left[ p^{(k-1)}\left(\mathcal{F}[c_{p}]\right) \right]$$
Using Fast Fourier Transforms, this is a very numerically efficient way to compute the coefficients you seek, as opposed to the brute force way of multiplying all the polynomials or convolving the coefficients directly.  (There may also be cases where the DFT can be computed analytically, so that this formula leads to closed form expressions for the coefficients, but this will not be the case in general.)

EDIT: A better derivation.
I recently found a much better way to derive the main result above.  As before, let $c_{p}(j)$ denote the $j^{th}$ coefficient of polynomial $p$.
By taking a (conventionally normalized) DFT we find $$\mathcal{F}[c_p](k) = \sum_j c_p(j) \left(e^{-2\pi i k/N}\right)^j = p\left(e^{-2\pi i k/N}\right)$$ where $N$ is the size of the DFT, and can be any integer larger than the degree of the polynomial $p$.  Hence, we have this expression for the coefficients of $p$:
$$c_p = \mathcal{F}^{-1}\left[ p\left(e^{-2\pi i k /N}\right) \right]$$
In particular, this works for iterated polynomials $p^{(n)}$:
$$\boxed{c_{p^{(n)}} = \mathcal{F}^{-1}\left[ p^{(n)}\left(e^{-2\pi i k /N}\right) \right]}$$
A few remarks:

A small variation on this works for analytic functions too.  You simply replace the DFT with a DTFT (i.e. let the size $N$ of the DFT $\rightarrow \infty$).  Then the expression for $c_p$ becomes $\mathcal{F}^{-1}\left[ p\left(e^{-2\pi i k}\right) \right]$, with $k\in [0,1]$ real valued.
The preceding expression can be viewed as Cauchy's integral formula.  Explicitly,
$$\begin{align*}
c_p(j) &amp; = \mathcal{F}^{-1}\left[ p\left(e^{-2\pi i k}\right) \right] \\
&amp; = \int_0^1 e^{2\pi i k j} p\left(e^{-2\pi i k}\right) dk \\
&amp; = \frac{1}{2\pi i}\oint \frac{p(z)}{z^{j+1}} dz \\
&amp; = \left. j! \frac{d^j p}{dz^j}\right|_{z=0}
\end{align*}
$$
where the third line follows from a change of variables and the last line is from Cauchy's derivative formula.
The boxed equation above can similarly be viewed as a discrete version of Cauchy's formula.


I implemented this in Mathematica to make sure it works, and sure enough it did. Something to watch out for if you try to implement this numerically is the normalization on your Fourier transforms.  We've assumed particular normalizations of the Fourier transform above, and if you use a different one (as is the default in Mathematica, it turns out), you need to multiply by a constant factor on every Fourier and inverse Fourier transform.  In Mathematica, that factor is $\sqrt{N}$, where $N$ is the length of vector you use to store the coefficients.  Define your own transform which is the default times $\sqrt{N}$, and your own inverse which is the default divided by $\sqrt{N}$, and everything will work out.  Also note that Mathematica's Fourier function is actually an inverse Fourier transform per standard conventions, and similarly InverseFourier is a Fourier transform.
The arbitrary polynomial can be presented via the vector of the coefficients
$\overrightarrow{C_n}$ in the form of
$$P_n\left(\overrightarrow{C_n},x\right) = 
\begin{pmatrix}c_{n}\\c_{n-1}\\ \dots \\ c_2 \\c_1\\c_0\end{pmatrix} 
\cdot
\begin{pmatrix}x^{n}\\x^{n-1}\\ \dots \\ x^2 \\ x\\ 1\end{pmatrix} 
= \overrightarrow{C_n}\cdot \overrightarrow{X_n}.\tag1$$

On the other hand, the arbitrary polynomial can be presented via the roots in the form of
$$Q_m(\overrightarrow {R_m},x) = q\prod\limits_{j=0}^m (x-r_j).\tag2$$

If to consider the presentation $(1)$ as an operator, the presentation $(2)$ as an operand, and the superposition as the resulted operator, then the operator transformation can be presented as the production of the modified polynomials in the form of

$$P_{mn}\left(\overrightarrow{C_n},Q_m(\overrightarrow {R_m},x)\right)
 = q\prod\limits_{j=0}^m \left(
\begin{pmatrix}c_{n}\\c_{n-1}\\ \dots \\ c_2 \\c_1\\c_0-r_j\end{pmatrix} 
\cdot
\begin{pmatrix}x^{n}\\x^{n-1}\\ \dots \\ x^2 \\ x\\ 1\end{pmatrix} 
\right).\tag3$$

For example, let 
$$f(x) = 2(x^2-3x+2),$$
then
$$P_2(x) = x^2-6x+4,\quad Q_2(x) = 2(x-1)(x-2),$$

$$f_2(x)=2(2x^2-6x+3)(2x^2-6x+2) = 8x^4 - 48x^3 + 92x^2 - 60x + 12,$$
$$f_3(x) = 2(8x^4 - 48x^3 + 92x^2 - 60x + 11)(8x^4 - 48x^3 + 92x^2 - 60x + 10)\\
= 4(32x^8 - 384x^7 + 1888x^6 - 4896x^5 + 7196x^4 - 6024x^3 + 2766x^2 - 630x + 55),$$
etc. (see also Wolfram Alpha calculations)

Therefore, the roots calculations looks the nesessary if we try to get useful closed form.

Now, let us use only the roots in the form $(2).$


  $$f_{(i)} = q^i\prod\limits_{j_1=0}^N (\prod\limits_{j_2=0}^N (\dots\prod\limits_{j_i=0}^N (x-r_{j_i})-\dots -r_{j_2})-r_{j_1}).\tag4$$


Ready

Calculations above gives result in the suitable form and does not require special algorithm. 

At the same time, applying of the FFT algoithm has some features:


as was pointed in the previous answer, using of the convolution for the polynomials mulplying requires to reorder one of them;
since the convolution algorithm calculates cyclic convolution, then the initial arrays should be filled zeros from the end to reserve the place for all coefficients of the resulting polynomial;
since DFT transform is a linear one, then, applying FFT templates for the constant, is possible to avoid the waste calls.