[Question]
  [
### Challenge
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.
### Examples
* 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\$
![](https://i.stack.imgur.com/3BIF3.png)
* 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$$
### Testcases
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]
```
  
[Answer]
# 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*
i=input()
p=Poly(i,var('x'));q=p*0+x
n=2
for _ in i[2:]:q-=compose(p,q).nth(n)*x**n;n+=1
print q.all_coeffs()
```
[Answer]
# Mathematica, 45 bytes
```
Normal@InverseSeries[#+O@x^(#~Exponent~x+1)]&
```
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.)
[Answer]
## JavaScript (ES6), 138 bytes
```
a=>a.reduce((r,_,i)=>[...r,i<2?i:a.map(l=>c=p.map((m,j)=>(r.map((n,k)=>p[k+=j]=m*n+(p[k]||0)),m*l+(c[j]||0)),p=[]),c=[],p=[1])&&-c[i]],[])
```
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.
[Answer]
# [Pari/GP](http://pari.math.u-bordeaux.fr/), 29 bytes
```
p->Pol(serreverse(p+O(x^#p)))
```
[Try it online!](https://tio.run/##LYrLDoJQDER/pcFNG9tExBcL@QXdE0xYXAwJwUkxRL/@CsJq5pwZ1N7aE7GhK0VYcX91PAT3MAYfAmN7489jAxGJNdB9GWQFwdv@zfMXQslsEmpWFlEqy4NmutddNfV0iVwvetaTHnXZVm3pRJb9sZL4Aw "Pari/GP – Try It Online")
]