# Recover polynomial $f(x)$ from $f^2(x)$

## 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 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]

• May 7, 2021 at 6:58
• Are we allowed to take the degree of the polynomial as an additional input? May 7, 2021 at 8:43
• I am surprised Wolfram language doesn't have a builtin for this, but I wasn't able to find one. May 7, 2021 at 16:00
• @Jonah It actually sort of does: reference.wolfram.com/language/ref/Decompose.html. If I had more time I would figure out how to turn it into a solution to the problem (e.g. see en.wikipedia.org/wiki/Polynomial_decomposition), but I'm leaving this comment in case anyone else wants to use it. May 7, 2021 at 21:44
• @soktinpk, That's interesting. Are you sure it can be used for this problem though? May 7, 2021 at 22:03

# Python 2, 638634608603572560508 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),N)
for F in(t,-t):
f=[F]
for i in R(1,N):f+=;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*F**~-n)/n/F**n
if n<2:f=g/-~F
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)&=\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 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 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_{ 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\$$.

• Welcome to CGCC & what a fantastic first answer! May 9, 2021 at 16:00
• (by the way, you can move all the test-cases into the 'footer' bit in TIO, so that it automatically only counts the bytes in your code; and I think you can also save 4 bytes by importing everything (*) from gmpy2...) May 9, 2021 at 16:02
• @DominicvanEssen oh nice thanks! yea it does save 4 bytes lemme update the answer May 9, 2021 at 16:05
• ...I think you might be able to lose 13 more bytes by deleting the r=lambda... line & just using iroot(...) in the two places that it was used... May 9, 2021 at 20:21

# 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 + ff_derivs**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((, 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:
# 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+s**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((,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:
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. 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 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, 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. 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.

• 684 bytes (doesn't work on TIO, I'm guessing because TIO is outdated) May 8, 2021 at 22:40
• @cairdcoinheringaahing Nice but I think it's because the np.round and the built-in round are conflicting May 8, 2021 at 22:43

• -3 bytes + a provably correct solution by shamelessly stealing the idea of trying increasing values of b from Arnauld's answer.
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! 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.

# J, 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!

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 ps 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.

# 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!

### 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.

• You mention in the explanation that you test all the values of $x$ between $0$ and $n-1$, where $n$ is the degree of $P$. However, I believe you should test at least $n+1$ different values of $x$ in order to be sure that $Q\circ Q=P$. As an example, consider the polynomials $x^2-x$ and $2x^2-2x$. If you evaluate them at $x=0$ and $x=1$ you get the same value, but they are not the same polynomial. May 9, 2021 at 14:13
• @Delfad0r The explanation was misleading: $n$ is the length of the input list, so the degree is $n-1$. (I've added your test cases.) May 9, 2021 at 17:24
• Ok, thanks! Very nice solution by the way, of all the "golfed" ones (including mine, sadly), this is the only one which, as far as I can tell, is provably right :) May 9, 2021 at 17:48

# Wolfram Language (Mathematica), 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!

Non-brute force. Just writes it as a equation in the coefficients, and solves it with Solve. Returns all solutions.

• Is there a chance it won't find one of the solutions? Jan 1 at 2:44
• @user974006 If that happens it would be a bug of Mathematica. Jan 1 at 8:39
• 88 bytes with coefficient list I/O
– att
Jan 2 at 23:29
• 86
– att
Jan 2 at 23:37

# Python/NumPy, 249 bytes

from numpy import*
from itertools import*
def f(p):
l=sqrt(abs(p));P=poly1d(p);t=P-[1,0];o=P-P.o*P//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!

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));P=poly1d(p);t=P-[1,0];o=P-P.o*P//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