[Question]
  [
Related: [Calculate \$f^n(x)\$](https://codegolf.stackexchange.com/q/137527/78410), [Polynomialception](https://codegolf.stackexchange.com/q/78334/78410)
## 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](https://en.wikipedia.org/wiki/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](/questions/tagged/code-golf "show questions tagged '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]
```
  
[Answer]
# [Python 2](https://docs.python.org/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):
	f=[F]
	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
	if n<2:f[1]=g[1]/-~F
	if all(e(f,e(f,i))==e(g,i)for i in R(N*N)):return f
```
[Try it online!](https://tio.run/##dVDBbuMgED3XX8GRsaE23jZKneWaow@9IrRyE@MlsbFLsNT00F/PDmlTZSNVgpnhzbw3M0zH8Hd05elk/DiQbpiOJbHDNPqQJs/SN65rEy99u503bTLJvhletg1xzEpRudySxm1VnymrzehJT6wjE3XcMgs6w4jZTMDoldJauUon7UXBsLfK06/Hkb1Xx/Qte2cGkm1ryIF2UCXESXt49YH2rUMAVrV0mVgF9kdaP46BNi9YqAoNrIaExAnWOAENjAek3xmp1hod4jZO9kwFFlYmk8hZGRxaIt1qfpgHapoNdZCfHY8NdwC5p8PcMxXB3f1mnF2ge4AouP8UjOtpSC91aq@/kzvE12lK64tazOw@f8gCcGqldJAaJXSs@@DY3eUYueTOGuJ@l1XMyQ5Nzj/WZ7Tpe9pSw@JFESlb2mFwtWKd1gCVb8PsHTGnyVsX8D@VYAQ/Krl@l/8B5SMjXCxvavgDI0t0eBfXKR71EC@/3A1tgeroHpEmxBOaJ9QpF2U0RewsYuInFtaK4heac7RElogsUeCID7dbFD8dDad/ "Python 2 – 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)&=\sum\_{i=0}^na\_ix^i\\
f(f(x))&=\sum\_{i=0}^na\_i\left(\sum\_{j=0}^na\_jx^j\right)^i\\
&=\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}\\
&=\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}\\
&=\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}\\
&=\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>n^2-n\$, then \$k\_1+\dots+k\_n=n\$, forcing \$k\_0=0\$ since \$a\_{>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<n\$ simplifies to
\begin{align\*}
&\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)\\
&=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)\\
&=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<i<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\_{<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\$.
[Answer]
# Python + numpy (758 bytes; non-brute force)
I thought it would be interesting to give a non-brute force answer to this question. Without further ado, here is the code:
```
import numpy as np
import numpy.polynomial.polynomial as Poly
from math import comb
def expand_bell(n_deriv, derivs, bell_arr):
  """Fill in the n_deriv'th diagonal of bell_arr
  Note that n=1 is the main diagonal, n=2 is the second diagonal etc.
  
  The result will have bell_arr[n, k] + B_(n, k)(...derivs) correct up to the n_deriv'th
  diagonal, where B_(n, k) is the partial bell polynomial"""
  for n in range(n_deriv, len(bell_arr)):
    k = n - n_deriv + 1
    for i in range(n_deriv):
      bell_arr[n,k] += comb(n - 1, i) * derivs[i] * bell_arr[n-i-1,k-1]
def get_nth_derivative(n, ff_derivs, bell_polys):
  """Get the nth derivative of f(f) in terms of the nth derivative of f
  given derivatives 0...n-1 of f evaluated at the fixed point
  For this we use the Faa di Bruno formula
  The return value is (m, b), coefficients such that m * f^(n)(x0) + b = ff^(n)(x0)
  
  Note that for n=1, ff'(x0) = f'(x0)^2, which is not linear; but for n>1 it is linear"""
  assert n >= 2
  b = 0
  for k in range(2, n):
    b += ff_derivs[k] * bell_polys[n, k]
  # k=1 and k=n
  m = ff_derivs[1] + ff_derivs[1]**n
  return m, b
x = Poly.Polynomial([0, 1])
def get_fixed_points(f):
  return (f - x).roots()
def poly_from_fixed_pt_and_first_deriv(ff, x0, x1):
  """Given the fixed point x0 of f, the first derivative x1, and ff(x) = f(f(x))
  All of the other derivatives of f are forced
  
  Compute these other derivatives and return the corresponding polynomial"""
  f_deriv = ff.deriv(2)
  derivs = [x0, x1]
  d = round(np.sqrt(ff.degree()))
  bell_arr = np.zeros((d+1, d+1), dtype=complex)
  ar = np.arange(d+1)
  bell_arr[(ar, ar)] = x1 ** ar
  for i in range(2, d+1):
    m, b = get_nth_derivative(i, derivs, bell_arr)
    deriv = (f_deriv(x0) - b) / m
    derivs.append(deriv)
    expand_bell(i, derivs[1:], bell_arr)
    f_deriv = f_deriv.deriv()
  return get_poly_from_derivs(x0, np.array(derivs))
  
def get_poly_from_derivs(x0, derivs):
  d = len(derivs)
  # Factorials
  fact = np.concatenate(([1], np.cumprod(np.arange(1, d))))
  coeffs = derivs / fact
  p = Poly.Polynomial(coeffs)
  return p(x - x0)
def round_polynomial(p):
  return Poly.Polynomial(np.round(np.real(p.coef)))
def functional_square_root(ff):
  if ff.degree() == 1 and ff.coef[1] == 1:
    # Special case because there is no fixed point for ff (or every point is a fixed point)
    return x + ff(0) / 2
  fixed_pts = get_fixed_points(ff)
  ff_deriv = ff.deriv()
  # Guess a fixed point
  for x0 in fixed_pts:
    first_deriv_st = ff_deriv(x0)**0.5
    # There are two possibilities for the first derivative
    for x1 in (first_deriv_st, -first_deriv_st):
      guess = round_polynomial(poly_from_fixed_pt_and_first_deriv(ff, x0, x1))
      if (ff.has_samecoef(guess(guess))):
        return guess
```
Try it on repl.it: <https://replit.com/@praneethkolichala/Functional-root-of-polynomials>
Golfed version:
```
import numpy as N
import numpy.polynomial.polynomial as Q
from math import comb
P,z,C=Q.Polynomial,N.zeros,complex
def G(f,x,y):
 D,d=f.deriv(2),round(N.sqrt(f.degree()))
 B,ar,s=z((d+1,d+1),dtype=C),N.arange(d+1),z(d+1,dtype=C)
 B[(ar,ar)],s[:2]=y**ar,[x,y]
 for i in range(2,d+1):
  m,b=s[1]+s[1]**i,sum(s[k]*B[i,k] for k in range(2, i))
  s[i]=(D(x)-b)/m
  for n in range(i,len(B)):B[n,n-i+1]=sum(comb(n-1,l)*s[1+l]*B[n-l-1,n-i] for l in range(i))
  D=D.deriv()
 F=N.hstack(([1],N.cumprod(N.arange(1,d+1))))
 return P(s/F)(P([0,1])-x)
def R(f):
 if f.degree()==1 and f.coef[1]==1:
  return P([0,1])+f(0)/2
 p,D = (f-P([0,1])).roots(),f.deriv()
 for x in p:
  t=D(x)**0.5
  for y in (t,-t):
   g=P(N.round(N.real((G(f,x,y)).coef)))
   if (f==g(g)):return g
```
# Explanation
*(Aside):*
At first, I was looking at [polynomial decomposition](https://en.wikipedia.org/wiki/Polynomial_decomposition). This *almost* solves the problem, because it decomposes a polynomial as $$h = u\_1 \circ u\_2 \circ \dots \circ u\_n$$
where \$u\_1, \dots u\_n\$ are not linear and cannot be further decomposed. Moreover, it turns out that these \$ u\_i \$'s are essentially unique up to linear factors. For example, for any linear polynomial \$\ell(x) = ax + b\$, clearly \$ u\_1 \circ u\_2 = (u\_1 \circ \ell) \circ (\ell^{-1} \circ u\_2)\$.
However, sometimes \$u\_1 \circ u\_2 = v\_1 \circ v\_2\$ are distinct decompositions. For example, \$T\_m \circ T\_n = T\_n \circ T\_m\$ where \$ T\_i \$ are the [Chebyshev polynomials](https://en.wikipedia.org/wiki/Chebyshev_polynomials) or even $$ x^n \circ x^sh(x^n) = x^sh(x)^n \circ x^n = x^{sn}h(x^n)^n $$
According to [this paper](https://arxiv.org/pdf/0807.3578.pdf), those are the only two cases, and all decompositions can be generated by taking some decomposition and applying these transformations or changing a polynomial by a linear factor. I couldn't figure out how to use these decompositions and combine them in a way to make \$f\circ f\$, although it seems very possible, and I'm sure someone else with some more mathematical ability could do it.
---
So how did I end up doing it? The first thing to notice was that composing a function with itself can only expand the set of fixed points, where a fixed point is an \$x\$ such that \$ f(x) = x \$. By the fundamental theorem of algebra, there exists a fixed point of \$ f(x) \$ in \$ \mathbb{C} \$.
Now, suppose we had \$ f\circ f \$ and a fixed point \$x\_0\$ of \$ f \$. We can compute the derivatives of \$ f\circ f \$ using [Faà di Bruno's formula](https://en.wikipedia.org/wiki/Fa%C3%A0_di_Bruno%27s_formula). For example, we have
$$ \frac{d}{dx}f(f(x)) = f'(f(x))f'(x) $$
Thus, evaluating at \$ x=x\_0 \$, we get \$ f'(x\_0)^2 \$, since \$ f(x\_0) = x\_0 \$. Since we can compute the derivative of \$ f(f(x)) \$ this gives us two possibilities for \$ f'(x\_0) \$.
In fact, it turns out that for \$ n \geq 2 \$, the value of \$ \frac{d^n}{dx^n} f(f(x)) \$ evaluated at \$ x = x\_0 \$ will be a linear polynomial in \$ f^{(n)}(x\_0) \$ (here, \$ f^{(n)} \$ is the \$n\$th derivative) in the Faà di Bruno formula, assuming \$ f'(x\_0), f''(x\_0), \dots f^{(n-1)}(x\_0) \$ are known and considered constants. Only for \$ n=1 \$ is it a quadratic; thus, we can solve exactly for the rest of the derivatives.
In sum, our procedure is to find the roots of \$ f(f(x)) = x \$ and iterate through them. Eventually, one of those roots will also have \$ f(x) = x \$. For each one, calculate the derivatives of \$ f(x) \$ assuming it were a fixed point of \$f(x)\$. Then, check if the polynomial with these derivatives actually matches \$ f(f(x)) \$ when it is composed with itself.
In principle, this method works even if the coefficients aren't integers. You would have to remove the `round_polynomial` method, however, and change `has_samecoef` to something that checks for approximately matching coefficients, so that floating point errors don't throw everything off.
[Answer]
# [Haskell](https://www.haskell.org/), ~~110~~… 97 bytes
* -3 bytes + a provably correct solution by shamelessly stealing the idea of trying increasing values of `b` from [Arnauld's answer](https://codegolf.stackexchange.com/a/225241/82619).
```
k#p=[q|b<-[1..],q<-mapM([-b..b]<$id)p,all((==).e p<*>e q.e q)[0..k*k]]!!0
e q x=foldl1((+).(*x))q
```
[Try it online!](https://tio.run/##XUvNjsIgEL73KTB6gO4MKY1uGlO8m6zxAQgHalGbIoL24GGffVns0Uzyzfd7Nc/ROpfSuAxSxd@uRSU41xBbvJlwoAo7zjvdroaeBTDOUSol45aEttxZEjOLTFWcj@Wo9WJRFdkgL3m@u94JSr8Yp@WLsZhuZvCyvxd5iQ9r@h@/3e6Pau8nnb3H4CcagDrrL9OVBIJEsGVgBXnvkhJQ6SJj/X71BlA0s8Y1NIANfGeFuQRYzzCH1efpv9PZmcsz4bH@Bw "Haskell – Try It Online")
The relevant function is `(#)`, which takes as input the degree `k` of the polynomial and the list `p` of coefficients, highest to lowest order.
## How?
```
e q x=foldl1((+).(*x))q
```
`e` is a simple function that evaluates a polynomial `q` at point `x`.
---
```
k#p=[q|b<-[1..],q<-mapM([-b..b]<$id)p,all((==).e p<*>e q.e q)[0..k*k]]!!0
```
`(#)` bruteforces all the polynomials `q` of degree at most `k`. We try all the possible coefficients between `-b` and `b`, starting from `b=1` and working our way up. As soon as we find a valid polynomial, we stop and return it. To check whether `q` is valid (i.e. `q(q(x))==p(x)`), instead of computing the coefficients of the composition, we try whether the equality holds for every `x` between `0` and `k*k`. Since both `q(q(x))` and `p(x)` have degree at most `k*k`, if they are equal for `k*k+1` values of `x` then they are the same polynomial.
[Answer]
# [J](http://jsoftware.com/), 80 67 bytes
Takes in the polynomials with lowest term first. Bruteforce (in the "generate all possible answers and filter them" way) and thus runs out of memory for the bigger test cases.
```
((-:[:+/[*1 0,+//.@:(*/)^:(]1<@-~#)~)"1#])>.@%:@#(>@,@{@#<@i:)>.&|/
```
[Try it online!](https://tio.run/##XY/BTsMwEETv@YpRI6jdpo7XTarUlMqiEifEgWsVfECJgAsfAOqvh0kCNFTWruTZfTu7793MzFvcesyRwcIzVgaHp4f7TqmVP/plflwIbLbMcxO8WuT62atadmF1SvVJzySt9d6EKx9StQ9Z@AzpLrx5atdfeaeT5uX1Ay02iBUqxAKSjJqFQNU3jdEs8zOq8k92f3pcozzrUSq4cqysEd2UmTqNKNfv05l2vTAkSR7vDIa2gpP6QY1hi4ijucBtLIPyluNki5K@QosJR6M47ECLFgU3tiWEnFTkSnKWuxe/3M/t9vL@yyfdNw "J – Try It Online")
Rather straight forward:
* `>.@%:@#(>@,@{@#<@i:)>.&|/` is a maybe too verbose way to generate `square(highest exponent)`-length lists with elements in the range of `-(max (abs coefficients)) … (max (abs coefficients))`. For each polynomial `p`:
* `+//.@:(*/)` multiplies two polynomials (with one side fixed to `p`)
* `^:(]1<@-~#)` do this `length of list - 1` times, keeping all the results. So we get `p^1, p^2, p^3, …`.
* `1 0,` prepend `p^0`: `1 + 0x + 0x^2 + 0x^3 …`.
* `[*` multiply by the coefficients, `a_0*p^0, a_1*p^1, a_2*p^2, a_3*p^3, …`
* `[:+/` as the `p`s are actually list, sum the rows together to get the final result.
* `-: … #` keep only the lists that have the original polynomial as the final result.
[Answer]
# JavaScript (ES7), ~~192~~ 178 bytes
Expects the coefficients from lowest to highest order and returns a list in the same format.
Pretty long, but very fast for the test cases.
```
p=>eval("r=g=(p,x)=>p.reduce((s,v,i)=>s+v*x**i);for(m=0;r;)for(a=Array(-~(p.length**.5)).fill(++m);p.some((_,x)=>g(p,x)-g(a,g(a,x)))&&!a.every((v,j)=>r=v+m?a[j]--&0:a[j]=m););a")
```
[Try it online!](https://tio.run/##jVPLbsIwELzzFS4HZCd2iENA0MhUPfQrUFRZwQlBgUQOjeDSX6e7Ie6hVSmWH6v1zox3be91p9vMls1JHOutuebq2qi16XRFx1YVijb8zNS6CazZfmSG0pZ3vARP63fe2fNKluS1pQcVJjZhaGr1aq2@UPFJm6Ayx@K087xgzliQl1VFff/AkiZo6wOQvffkRS8iCqo5jjNjbDJ50oHpjL1Q2vE9BFnV@YcXvdmnQkzCZzQUMLFEj9k12YwI2UhOwpSTx9p0OgBIbclGgJmlA0n0MMtAIsMeGs05EXL5GBihGD9zqiLmZAkLjMV9ikFVRJzc0Hj68Obol7twRDuAdOILyBqWOYhLuYJpBaeJFhFOIVYEs5RIPIijNifxT3yMxYAN2VtLiJKIlyGkGvdldfJ4WDHDtB1H@Ff/lY67u2H7@wadw1X0Vo8HnoWrqHRovAgJScV9B4n/0NEgl47SUQDf4E1nO9oQtSZZfWzrygRVXdCc4keCZ90aCq@cXb8A "JavaScript (Node.js) – Try It Online")
### How?
Given a polynomial \$P\$ of length \$n\$ (i.e. of degree \$n-1\$) as input, this builds all polynomials \$Q\$ of length \$\lfloor\sqrt{n}\rfloor+1\$ with integer coefficients in \$[-m\dots m]\$ such that \$Q(Q(x))=P(x)\$ for all \$x\in[0\dots n-1]\$. We start with \$m=1\$ and widen the search window until a solution is found.
[Answer]
# [Wolfram Language (Mathematica)](https://www.wolfram.com/wolframscript/), 86 bytes
-10 bytes thanks to @att.
```
(t=a/@(r=Range[0,Tr[1^#]^.5]))/.Solve[(p=t.#^r&)@p@x~CoefficientList~x==#,t,Integers]&
```
[Try it online!](https://tio.run/##fY9Ba4MwFMfv/RSBgCiLmpcap4cMYafBDmPbLShIia3Q2mLDCIz1q7sXW3foYaf3S97//V5yaO3OHFrbb9qpI2oKrWrTKhzVeztsjebsc9TQ0LpJZB1FafJx3H8ZHZ6UTWgzBlF1qtzl@Wi6rt/0ZrCv/dlenFKUWfYyWLM147kOprexH6ymJH4iXXUX15S5uiYBSavVt2Mr4sgDERxBSOSYQMEI3jYZcoZ1jf0Cq8BzMadzDMfXgLg1xHWm9F2shXcCwiOC9NM5AkCJJD2V2bwBU7m47RA5n2X@NTDvAfhXmy1avl60ctFCsWjhTwtcztrMfxalP9Mv "Wolfram Language (Mathematica) – Try It Online")
Non-brute force. Just writes it as a equation in the coefficients, and solves it with `Solve`. Returns all solutions.
[Answer]
# Python/NumPy, 249 bytes
```
from numpy import*
from itertools import*
def f(p):
 l=sqrt(abs(p[0]));P=poly1d(p);t=P-[1,0];o=P-P.o*P[0]//2
 for s,c in product((-l,l),combinations(t.r,int(sqrt(P.o)))):a=poly1d(int_(around((s-(P.o==1))*poly(c))))+[1,0];o=[o,a][a(a)==P]
 return o.c
```
[Attempt This Online!](https://ato.pxeger.com/run?1=hVPLcpswFF10x1fcySZSImzEwyXO8AFexXsPZAiGlhmMqBAt_pZusmln-kntn3TXewHHuEmnHoOQdO455x7B1-_N0XxU9fPzt84Udvjzd6HVAeru0ByhPDRKmxtrWCpNro1SVfuyvM8LKFjD1xZUUftJG5Y-tazZOTHn99uoUdVR7nH_3kRbeyeFE98rfNou1M0WQcula0GhNLQig7KGRqt9lxnG7EpUXGTq8FTWqSlV3TKz0KKsDRtUkIDjb52eJHDnkaVadfWesdYmQBRJzm9on2UEvj3p75RI413KUh5F29gCnZtO16AW2RTAj7E9aI-t9eVjWeUgscGNeIh2-ee0QrGmM4wv2qYqDbteXvOdLePTdH05VRrnmMfQ5yN1eSXdK5Q1-ois2DN1tREF2tfpkW04Fw_cgrzP8sasx210ssCFx7IuFOO4O1kYDT__eqc3tLCGHpaQqbwoyqzMa9OuAbsGJ7YeOjMCBNh4oRe7h1vIBltpjSddm_xDrnHlbQYqwb4EZLFlvajdguu8WeDONREm34ZJ58zmBj3YIMPXQDdA1zI8Mw5I7zWQcN7MX-Ijzu8TDx2EfeLiLCQ3KywledsXEOKA12rmdwC6BPROQFfAjNgemd0JSSCKxhmBw3Cms0eUnKPk3OQdOeqTkNKUffIexwDdrig2edcnAT3c-aSJiJU79uOuHCKmE5BDwtKdzK4wfhwC7ArrBdUKqqObQ0dDwUt33vBAONJ5ROZPVNS2AP9_bv3JreNNboPJrQwnt_LkVjrBIOD8ZdanlwG15PAUorAkswhHfeciTe8cvUcvgn9KlsK3PTrWS8OjkvOv_2UQw6eReJdF59f_pWb89v4A)
Not very golfed yet. But perhaps good as a readable baseline.
Also, not sure this will work for every single edge case.
It works by solving for fixed points of f(f(x)), i.e. roots of f(f(x))-x (the case f(f(x))=x+c is treated separately). The fixed points of f(x) are a subset.
We then simply try all subsets of correct size.
Here is the code again with a few comments. Note that the (deprecated) poly1d class has the very golfing friendly attributes `.c` (coefficients as array) `o` (order) and `.r` (roots). Also, note that indexing a `poly1d` starts from the **end** of the underlying array.
```
from numpy import*
from itertools import*
def f(p):
 # l: leading coefficient of f
 # t: fixed point equation, i.e. P-x ([1,0] is x)
 # o: output; first set at the appropriate value for the
 #    no fixed-point case; will be overwritten otherwise
 l=sqrt(abs(p[0]));P=poly1d(p);t=P-[1,0];o=P-P.o*P[0]//2
 # try: +/- and all subsets (of size sqrt of input order)
 # of the roots of t
 for s,c in product((-l,l),combinations(t.r,int(sqrt(P.o)))):a=poly1d(int_(around((s-(P.o==1))*poly(c))))+[1,0];o=[o,a][a(a)==P]
 return o.c
```
]