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

• – alephalpha May 7 at 6:58
• Are we allowed to take the degree of the polynomial as an additional input? – Delfad0r May 7 at 8:43
• I am surprised Wolfram language doesn't have a builtin for this, but I wasn't able to find one. – Jonah May 7 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. – soktinpk May 7 at 21:44
• @soktinpk, That's interesting. Are you sure it can be used for this problem though? – Jonah May 7 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! – Dominic van Essen May 9 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...) – Dominic van Essen May 9 at 16:02
• @DominicvanEssen oh nice thanks! yea it does save 4 bytes lemme update the answer – Ariana May 9 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... – Dominic van Essen May 9 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) – caird coinheringaahing May 8 at 22:40
• @cairdcoinheringaahing Nice but I think it's because the np.round and the built-in round are conflicting – soktinpk May 8 at 22:43

# 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. – Delfad0r May 9 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.) – Arnauld May 9 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 :) – Delfad0r May 9 at 17:48

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

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