Decompose Polynomials

Question

Given an integral polynomial of degree strictly greater than one, completely decompose it into a composition of integral polynomials of degree strictly greater than one.

Details

• An integral polynomial is a polynomial with only integers as coefficients.
• Given two polynomials p and q the composition is defined by (p∘q)(x):=p(q(x)).
• The decomposition of an integral polynomial p is a finite ordered sequence of integral polynomials q1,q2,...,qn where deg qi > 1 for all 1 ≤ i ≤ n and p(x) = q1(q2(...qn(x)...)), and all qi are not further decomposable. The decomposition is not necessarily unique.
• You can use e.g. lists of coefficients or built in polynomial types as input and output.
• Note that many builtins for this task actually decompose the polynomials over a given field and not necessarily integers, while this challenge requires a decomposition integer polynomials. (Some integer polynomials might admit decomposition into integer polynomials as well as decomposition that contain rational polynomials.)

Examples

x^2 + 1
[x^2 + 1] (all polynomials of degree 2 or less are not decomposable)
x^6 - 6x^5 + 15x^4 - 20x^3 + 15x^2 - 6 x - 1
[x^3 - 2, x^2 - 2x + 1]
x^4 - 8x^3 + 18x^2 - 8x + 2 
[x^2 + 1, x^2 - 4x + 1]
x^6 + x^2 + 1
[x^3 + x + 1, x^2]
x^6
[x^2, x^3]
x^8 + 4x^6 + 6x^4 + 4x^2 + 4 = (x^2 + 1)^4 + 3
[x^2 + 3, x^2, x^2 + 1]
x^6 + 6x^4 + x^3 + 9x^2 + 3x - 5
[x^2 + x - 5, x^3 + 3*x], [x^2 + 5*x + 1, x^3 + 3*x - 2]

Use Maxima for generating examples: Try it online!

Some decomposition algorithms can be found here and here.

Answer 1

Pari/GP, 84 bytes

f(p)=[if(q'',[f(q),r],p)|r<-x*divisors(p\x),r''&&p==subst(q=substpol(p,r,x),x,r)][1]

Based on the algorithm described here.

Try it online!

• flawr: Do you check (or filter out) whether you actually get a decomposition into integral polynomials? (I'm asking because the algorithms in the linked paper describe the factorization over some field, and I don't know any Pari/GP.)
• alephalpha: @flawr I'm using the second algorithm in the paper, which always returns integral polynomials when the input is integral. In fact, the divisors function in Pari/GP always returns primitive polynomials when it takes an integral polynomial. It can be proved that if p=q∘r, where p and r are integral, and r is primitive with r(0)=0, then q must also be integral. Here p, q, r correspond to f, g, h in the paper.

Answer 2

Wolfram Language (Mathematica), 29 bytes

Decompose[#/.x->x+a,x]/.a->0&

Try it online!

I have the example set up here to compose a random polynomial from random quadratics (or less), expand it out, and then try to decompose it.

It's necessary to complicate the polynomial with dummy variable (a) since the built-in will not attempt to decompose a monomial.

I notice that the answer often has much larger coefficients than in the original composition, but they are indeed always integers.

• flawr: Where did you find the information that Decompose[] will always return integral polynomials (if fed with integer polynomials)? When discussing in chat recently we could not find anything about that.
• Kelly Lowder: Do Options@Decompose and it will tell you {Modulus->0}. Now look up Modulus and you'll see "The setting Modulus->0 specifies the full ring [DoubleStruckCapitalZ] of integers. "
• flawr: Ah that is nice, thanks for elaborating!