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