Related: Calculate \$f^n(x)\$, Polynomialception


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]
  • \$\begingroup\$ Related: codegolf.stackexchange.com/q/157844/9288 \$\endgroup\$
    – alephalpha
    Commented May 7, 2021 at 6:58
  • \$\begingroup\$ Are we allowed to take the degree of the polynomial as an additional input? \$\endgroup\$
    – Delfad0r
    Commented May 7, 2021 at 8:43
  • \$\begingroup\$ I am surprised Wolfram language doesn't have a builtin for this, but I wasn't able to find one. \$\endgroup\$
    – Jonah
    Commented May 7, 2021 at 16:00
  • \$\begingroup\$ @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. \$\endgroup\$
    – soktinpk
    Commented May 7, 2021 at 21:44
  • \$\begingroup\$ @soktinpk, That's interesting. Are you sure it can be used for this problem though? \$\endgroup\$
    – Jonah
    Commented May 7, 2021 at 22:03

7 Answers 7


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*
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):
 for F in(t,-t):
	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!


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\$.

  • 2
    \$\begingroup\$ Welcome to CGCC & what a fantastic first answer! \$\endgroup\$ Commented May 9, 2021 at 16:00
  • \$\begingroup\$ (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...) \$\endgroup\$ Commented May 9, 2021 at 16:02
  • \$\begingroup\$ @DominicvanEssen oh nice thanks! yea it does save 4 bytes lemme update the answer \$\endgroup\$
    – Ariana
    Commented May 9, 2021 at 16:05
  • \$\begingroup\$ ...I think you might be able to lose 13 more bytes by deleting the r=lambda... line & just using iroot(...)[0] in the two places that it was used... \$\endgroup\$ Commented 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[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
    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
def G(f,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))
  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))
 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:
  for y in (t,-t):
   if (f==g(g)):return g


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

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

Haskell, 110 97 bytes

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


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.


Try it online!


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.

  • \$\begingroup\$ 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. \$\endgroup\$
    – Delfad0r
    Commented May 9, 2021 at 14:13
  • 1
    \$\begingroup\$ @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.) \$\endgroup\$
    – Arnauld
    Commented May 9, 2021 at 17:24
  • \$\begingroup\$ 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 :) \$\endgroup\$
    – Delfad0r
    Commented May 9, 2021 at 17:48

Wolfram Language (Mathematica), 86 bytes

-10 bytes thanks to @att.


Try it online!

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

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

Python/NumPy, 249 bytes

from numpy import*
from itertools import*
def f(p):
 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
 # 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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.