Given a polynomial \$p\$ with real coefficients of order \$1\$ and degree \$n\$, find another polynomial \$q\$ of degree at most \$n\$ such that \$(p∘q)(X) = p(q(X)) \equiv X \mod X^{n+1}\$, or in other words such that \$p(q(X)) = X + h(X)\$ where \$h\$ is an arbitrary polynomial with \$ord(h) \geqslant n+1\$. The polynomial \$q\$ is uniquely determined by \$p\$.

For a polynomial \$p(X) = a_nX^n + a_{n+1}X^{n+1} + ... + a_m X^m\$ where \$n \leqslant m\$ and \$a_n ≠ 0\$,\$a_m ≠ 0\$, we say \$n\$ is the order of \$p\$ and \$m\$ is the degree of \$p\$.

Simplification: You can assume that \$p\$ has integer coefficients, and \$a_1 = 1\$ (so \$p(X) = X + \text{[some integral polynomial of order 2]}\$). In this case \$q\$ has integral coefficients too.

The purpose of this simplification is to avoid the issues with floating point numbers. There is however a non-integral example for illustration purposes.


  • Consider the Taylor series of \$\exp(x)-1 = x + x^2/2 + x^3/6 + x^4/24 + ...\$ and \$\ln(x+1) = x - x^2/2 + x^3/3 - x^4/4 + ... \$ then obviously \$\ln(\exp(x)-1+1)= x\$. If we just consider the Taylor polynomials of degree 4 of those two functions we get with the notation from below (see testcases) \$p = [-1/4,1/3,-1/2,1,0]\$ and \$q = [1/24, 1/6, 1/2, 1,0]\$ and \$(p∘q)(X) \equiv X \mod X^5\$
  • Consider the polynomial \$p(X) = X + X^2 + X^3 + X^4\$. Then for \$q(X) = X - X^2 + X^3 - X^4\$ we get

$$(p \circ q)(X) = p(q(X)) = X - 2X^5 + 3X^6 - 10X^7 +...+ X^{16} \equiv X \mod X^5$$


Here the input and output polynomials are written as lists of coefficients (with the coefficient of the highest degree monomial first, the constant term last):

p = [4,3,2,0];  q=[0.3125,-.375,0.5,0]

Integral Testcases:

p = [1,0]; q = [1,0]

p = [9,8,7,6,5,4,3,2,1,0]; q = [4862,-1430,429,-132,42,-14,5,-2,1,0]

p = [-1,3,-3,1,0]; q = [91,15,3,1,0]

4 Answers 4


Python 2 + sympy, 128 bytes

We locally invert the polynomial by assuming that q(x) = x, composing it with p, checking the coefficient for x2, and subtracting it from q. Let's say the coefficient was 4, then the new polynomial becomes q(x) = x - 4x2. We then again compose this with p, but look up the coefficient for x3. Etc...

from sympy import*
for _ in i[2:]:q-=compose(p,q).nth(n)*x**n;n+=1
print q.all_coeffs()

Mathematica, 45 bytes


Yeah, Mathematica has a builtin for that....

Unnamed function taking as input a polynomial in the variable x, such as -x^4+3x^3-3x^2+x for the last test case, and returning a polynomial with similar syntax, such as x+3x^2+15x^3+91x^4 for the last test case.

#+O@x^(#~Exponent~x+1) turns the input # into a power series object, truncated at the degree of #; InverseSeries does what it says; and Normal turns the resulting truncated power series back into a polynomial. (We could save those initial 7 bytes if an answer in the form x+3x^2+15x^3+91x^4+O[x]^5 were acceptable. Indeed, if that were an acceptable format for both input and output, then InverseSeries alone would be a 13-byte solution.)


JavaScript (ES6), 138 bytes


Port of @orlp's answer. I/O is in the form of arrays of coefficients in reverse order i.e. the first two coefficients are always 0 and 1.


Pari/GP, 29 bytes


Try it online!


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.