Coefficients of $p\circ p$ for polynomial $p(x)=\sum_{k=0}^d a_kx^k$ Let $p(x)=\sum_{k=0}^d a_kx^k$ be a polynomial of degree $d$. Then
\begin{align*}
(p\circ p)(x) &= \sum_{k=0}^da_k\left(\sum_{j=0}^da_j x^j\right)^k\\[3pt]
&=\sum_{k=0}^d \sum_{\substack{k_0+\cdots+k_d=k\\[3pt]k_0,\dots,k_d\geqslant0}}a_k\frac{k!}{k_0!k_1!\cdots k_d!} x^{k_1+2k_2+3k_3+\cdots+d k_d}\prod_{t=0}^da_t^{k_t}
\end{align*}
by the multinomial theorem.
Is there a nicer way I can express the coefficients of $p\circ p$? Where can I read more about this polynomial and study its properties?
 A: TL;DR: There are several ways to organize the computation of these coefficients.  None of them give a fully closed form expression (as expected), so which one is best depends on the application.  The most efficient for implementation on a computer is the Fourier method explained below.

Let me provide multiple approaches to this problem.
Approach 1: Direct assault.
Let's see what happens if we just try to evaluate the coefficient of $x^m$  in $\sum_k a_k \left(\sum_j a_j x^j\right)^k$ directly, starting with small values of $m$.

*

*For $m=0$, only the $j=0$ term of the inner sum contributes, so the coefficient of $x^m$ is $\sum_k a_k (a_0)^k = p(a_0) = p(p(0))$.

*For $m=1$, in each term of the outer sum we are looking to have one power of $x$.  When expanding out $\left(\sum_j a_j x^j\right)^k$, there are $k$ ways to pick a term $a_1 x^1$, and the other $k-1$ factors in the power must be $a_0 x^0$.  Hence the coefficient of $x^1$ is $\sum_k a_k \left(k a_1 a_0^{k-1}\right)$.  This can be simplified by noting that $\sum_k k a_k y^{k-1}$ is precisely $p'(y)$, so that the coefficient of $x^1$ becomes $a_1 p'(a_0)$.  Since $a_1=p'(0)$, this can in turn be written as $p'(p(0)) p'(0)$.

Do you see the pattern yet?  All we are doing is computing Taylor series coefficients.  In general, the coefficient of $x^m$ will be $$\left.\frac{1}{m!} \frac{d^m p(p(x))}{dx^m}\right|_{x=0}$$
So e.g. for $m=2$ we will have $\frac{1}{2}\left(p''(p(0)) p'(0)^2 + p'(p(0)) p''(0)\right)$.
Approach 2: Fourier methods.
(Note: This is a recapitulation of this post, also linked in comments above.)
When you multiply two polynomials, their coefficients get "convolved".  So if $c_p(j)$ denotes the $j^{th}$ coefficient of polynomial $p$, then the coefficient of a product $pq$ is given by $c_{pq}(j) = \sum_k c_p(k) c_q(j-k) = (c_p * c_q)(j)$, where $*$ denotes convolution.  If you take a power of a polynomial, like $p^k$, then its coefficients get convolved with themselves $k$ times, which is denoted $c_p^{*k}$.  So with this notation we can write $c_{p\circ p}(j) = \sum_k c_p(k) c_p^{*k}(j)$.
The most efficient way to handle convolutions is via the Fourier transform.  If you are unfamiliar with Fourier analysis, this may seem mysterious, but if you actually want to compute coefficients of an iterated polynomial e.g. on a computer, this is almost certainly the best way.  The main result is (as shown in the linked post above) that $$c_{p\circ p} (j) = \mathcal{F}^{-1} \left[ p\circ p\left(e^{-2\pi i k/N}\right)\right](j)$$ where $\mathcal{F}^{-1}$ is an inverse discrete Fourier transform of length $N$.  This gives an efficient, and easy to implement, way to compute coefficients of iterated polynomials.
A pseudocode recipe for implementing this would be something like the following:
m = ... # set this to the degree of polynomial p
N = m^2 + 1 # chosen to be total number of coefficients of p(p(x))
y = [p(p(exp(-2*pi*i*k/N))) for k=0:(N-1)]
coefficients = InverseDFT(y) # these are the coefficients of p(p(x))

(Note that this assumes a conventional normalization of your DFT, so you may need to e.g. divide by $\sqrt{N}$ if your software uses a different normalization.)
Approach 3: Cauchy's formula.
Writing out the above Fourier method more explicitly yields a formula that might be easier to think about.  Let $\omega$ be a primitive $N^{th}$ root of unity, with $N \ge \deg(p\circ p)$.  Then the coefficient $c_{p\circ p}(j)$ can be written
$$c_{p\circ p}(j) = \frac{1}{N} \sum_{k=0}^{N-1} \omega^{-jk} p(p(\omega^k))$$
This is essentially Cauchy's integral formula specialized to the case at hand.  It follows directly from the Fourier formula above (cf. this related MO post).
A: Not quite the answer you're looking for, but for fun taking a cubic $p(z) = a_0 + a_1 z + a_2 z^2 + a_3 z^3$, $p(p(z)) = c_0 + \cdots + c_9 z^9$ where:
\begin{align}
c_0 &=  a_{0}^{3} a_{3} + a_{0}^{2} a_{2} + a_{0} a_{1} + a_{0} \\
c_1 &=  3 \, a_{0}^{2} a_{1} a_{3} + 2 \, a_{0} a_{1} a_{2} + a_{1}^{2} \\
c_2 &=  3 \, a_{0} a_{1}^{2} a_{3} + 3 \, a_{0}^{2} a_{2} a_{3} + a_{1}^{2} a_{2} + 2 \, a_{0} a_{2}^{2} + a_{1} a_{2} \\
c_3 &=  a_{1}^{3} a_{3} + 6 \, a_{0} a_{1} a_{2} a_{3} + 3 \, a_{0}^{2} a_{3}^{2} + 2 \, a_{1} a_{2}^{2} + 2 \, a_{0} a_{2} a_{3} + a_{1} a_{3} \\
c_4 &=  3 \, a_{1}^{2} a_{2} a_{3} + 3 \, a_{0} a_{2}^{2} a_{3} + 6 \, a_{0} a_{1} a_{3}^{2} + a_{2}^{3} + 2 \, a_{1} a_{2} a_{3} \\
c_5 &=  3 \, a_{1} a_{2}^{2} a_{3} + 3 \, a_{1}^{2} a_{3}^{2} + 6 \, a_{0} a_{2} a_{3}^{2} + 2 \, a_{2}^{2} a_{3} \\
c_6 &=  a_{2}^{3} a_{3} + 6 \, a_{1} a_{2} a_{3}^{2} + 3 \, a_{0} a_{3}^{3} + a_{2} a_{3}^{2} \\
c_7 &=  3 \, a_{2}^{2} a_{3}^{2} + 3 \, a_{1} a_{3}^{3} \\
c_8 &=  3 \, a_{2} a_{3}^{3} \\
c_9 &=  a_{3}^{4} \\
\end{align}
