Recover polynomial \$f(x)\$ from \$f^2(x)\$ Related: Calculate \$f^n(x)\$, Polynomialception
Challenge
Given a polynomial \$f(x) = a_0 + a_1 x + a_2 x^2 + \cdots + a_k x^k\$ of order \$k\$, we can compute its composition with itself \$f\left(f(x)\right) = f \circ f(x) =: f^2(x)\$, which evaluates to another polynomial of order \$k^2\$.
The challenge is to reverse the process. Assuming \$f(x)\$ is an integer polynomial (\$a_0, a_1, \cdots, a_k \in \mathbb{Z}\$), take the polynomial \$f^2(x)\$ and recover the polynomial \$f(x)\$. This is equivalent to finding a functional square root of the given polynomial.
The polynomials may be input and output as a list of coefficients (highest-order term going first or last) or a built-in polynomial object. You may assume \$k \ge 1\$, i.e. \$f(x)\$ is not a constant. It is guaranteed that a solution exists. If multiple distinct \$f(x)\$s have the same \$f^2(x)\$, your program may output any nonempty subset of all solutions.
Standard code-golf rules apply. The shortest code wins.
Test cases
Coefficient list is given in highest-order-term-first order.
Input: x / coefficients: [1, 0]
Output: x, -x, or -x + c for any integer c / coefficients: [1, 0] or [-1, c]

Input: x + 20 / coefficients: [1, 20]
Output: x + 10 / coefficients: [1, 10]

Input: 25x - 18 / coefficients: [25, -18]
Output: 5x - 3 / coefficients: [5, -3]

Input: x^4 - 4x^3 + 8x^2 - 8x + 6 / [1, -4, 8, -8, 6]
Output: x^2 - 2x + 3 / [1, -2, 3]

Input: -x^4 - 2x^2 - 2 / [-1, 0, -2, 0, -2]
Output: -x^2 - 1 / [-1, 0, -1]

Input: x^9 + 6x^8 + 21x^7 + 58x^6 + 119x^5 + 194x^4 + 262x^3 + 260x^2 + 201x + 112 / [1, 6, 21, 58, 119, 194, 262, 260, 201, 112]
Output: x^3 + 2x^2 + 3x + 4 / [1, 2, 3, 4]

Input: x^9 + 6x^8 + 21x^7 + 54x^6 + 103x^5 + 154x^4 + 182x^3 + 160x^2 + 105x + 40 / [1, 6, 21, 54, 103, 154, 182, 160, 105, 40]
Output: -x^3 - 2x^2 - 3x - 4 / [-1, -2, -3, -4]

Input: x^9 / [1, 0, 0, 0, 0, 0, 0, 0, 0, 0]
Output: x^3 or -x^3 / [1, 0, 0, 0] or [-1, 0, 0, 0]
Python 2, 638 634 608 603 572 560 508 470 bytes (without brute force as well)
Tthanks to @randomdude999 for the golfing down the code by like 100 bytes as well as everyone who has helped point out smaller golfing tricks!

from gmpy2 import*
R=range
r=reduce
p=lambda n,i=1:n/i and[l+[i]for l in p(n-i,i)]+p(n,i+1)or[[]][n:]
e=lambda f,x:r(lambda y,z:y*x+z,f)
def s(g):
 n=isqrt(len(g));N=n+1;t,_=iroot(abs(g[0]),N)
 for F in(t,-t):
&#9;f=[F]
&#9;for i in R(1,N):f+=[0];f[i]=(g[i]-sum(fac(n)/fac(n-len(j))/r(mul,[fac(j.count(k))for k in R(i+1)])*r(mul,[f[k]for k in j])*F**(N-len(j))for j in p(i))-(i==n)*f[1]*F**~-n)/n/F**n
&#9;if n&lt;2:f[1]=g[1]/-~F
&#9;if all(e(f,e(f,i))==e(g,i)for i in R(N*N)):return f

Try it online!
Explanation
In the code, the functions

m computes the multinomial coefficient \$\binom{n}{k_0,k_1,\dots,n-\sum k}\$
p list out all partitions of n (code from https://codegolf.stackexchange.com/a/84192/61039 (without this code was 601 bytes))
e evaluates our polynoimal
s does the actual functional square root

The way it works is by computing the coefficients of \$f\$ (of degree \$n\$) using the \$n+1\$th highest coefficients of \$f\circ f\$, then checks if the computation is correct as multiple polynomials can lead to the same highest coefficients. The way the computation works is by considering the following series of manipulations:
\begin{align*}
f(x)&amp;=\sum_{i=0}^na_ix^i\\
f(f(x))&amp;=\sum_{i=0}^na_i\left(\sum_{j=0}^na_jx^j\right)^i\\
&amp;=\sum_{i=0}^na_i\sum_{k_0+k_1+\dots+k_n=i}\binom{i}{k_0,k_1,\dots,k_n}\prod_{j=0}^n a_j^{k_j}x^{jk_j}\\
&amp;=\sum_{i=0}^na_i\sum_{k_0+k_1+\dots+k_n=i}\binom{i}{k_0,k_1,\dots,k_n}\left(\prod_{j=0}^na_j^{k_j}\right)x^{\sum_{j=0}^njk_j}\\
&amp;=\sum_{k_0+k_1+\dots+k_n\leq n}a_{k_0+k_1+\dots+k_n}\binom{k_0+k_1+\dots+k_n}{k_0,k_1,\dots,k_n}\left(\prod_{j=0}^na_j^{k_j}\right)x^{\sum_{j=0}^njk_j}\\
&amp;=\sum_{i=0}^{n^2}\left(\sum_{\sum_{j=0}^njk_j=i,k_i\geq0}a_{k_0+k_1+\dots+k_n}\binom{k_0+k_1+\dots+k_n}{k_0,k_1,\dots,k_n}\left(\prod_{j=0}^na_j^{k_j}\right)\right)x^i\\
\end{align*}
From here, if \$i&gt;n^2-n\$, then \$k_1+\dots+k_n=n\$, forcing \$k_0=0\$ since \$a_{&gt;n}=0\$. Furthermore,
$$\sum_{j=0}^njk_j=n\sum_{j=0}^nk_j-\sum_{j=0}^njk_{n-j}=n^2-\sum_{j=0}^njk_{n-j}$$
so the coefficient of \$x^{n^2-i}\$ for \$i&lt;n\$ simplifies to
\begin{align*}
&amp;\hphantom{=}a_n\sum_{\substack{\sum_{j=1}^njk_j=n^2-i\\k_0=0\\k_i\geq0}}\binom{n}{k_0,k_1,\dots,k_n}\left(\prod_{j=0}^na_j^{k_j}\right)\\
&amp;=a_n\sum_{\substack{\sum_{j=0}^njk_{n-j}=i\\\sum_{j=0}^nk_j=n\\k_i\geq0}}\binom{n}{k_0,k_1,\dots,k_n}\left(\prod_{j=0}^na_j^{k_j}\right)\\
&amp;=a_n\sum_{\substack{\sum_{j=0}^njk_j=i\\k_0=n-\sum_{j=1}^nk_j\\k_i\geq0}}\binom{n}{k_0,k_1,\dots,k_n}\left(\prod_{j=0}^na_{n-j}^{k_j}\right)
\end{align*}
and the \$\sum_{j=1}^njk_j=i\$ is really asking for partitions of \$i\$.
To extend this analysis to \$i=n\$ to get \$a_0\$, we simply notice that we need to include the case that \$1(n-1)+(n-1)n=n^2\$, i.e. we simply add \$a_0^{n-1}a_1\$.
The reason why this helps is this expansion gives us

\$x^{n^2}\$ coefficient of \$f\circ f\$ is \$a_n^{n+1}\$
\$x^{n^2-i}\$ coefficient is linear in \$a_i\$ for \$0&lt;i&lt;n\$ as the only term containing \$a_i\$ is \$na_0^na_i\$

As all terms in the coefficient of \$x^{n^2-i}\$ in \$f\circ f\$ only depends on \$a_{&lt;i}\$ other than \$na_0^na_i\$, we can easily compute \$a_i\$. The only issue is when \$n\$ is odd, which gives us \$2\$ solutions for \$a_n\$.