# Polynomial decomposition algorithm with a hint

If $p,f,g$ are polynomials of degree 2 or higher and $p(x)=f(g(x))$ we have a polynomial decomposition of $p$. When Mathematica (as noted below this was actually on Wolfram Alpha) says decomposition not found'' does it mean it does not exist? Or is it that the algorithm is not conclusive? In the latter case how can I provide $g$ and ask Mathematica to check for existence of $f$?

Edit:

In response to my second question: having a guess for $g$ I suppose I could do repeated long division to recover $f$.

Here is a reference for polynomial decomposition.

• In response to what input do you receive this message? – Oleksandr R. Nov 1 '15 at 16:35
• @OleksandrR. decompose 19456*z^8+98368*z^7+195820*z^6+140530*z^5+60493*z^4+14413*z^3+2149*z^2+196*z+16 .This was was done on Wolfram alpha – Maesumi Nov 1 '15 at 16:38
• Wolfram|Alpha is decidedly not Mathematica, so please do not confuse them. However, Decompose[19456*z^8 + 98368*z^7 + 195820*z^6 + 140530*z^5 + 60493*z^4 + 14413*z^3 + 2149*z^2 + 196*z + 16, z] also does not yield any decomposition: it returns the input as its answer without producing any messages. I would take this as indicating that the decomposition does not exist. – Oleksandr R. Nov 1 '15 at 16:50
• If Decompose does not find a polynomial decomposition of a univariate polynomial with rational (or integer) coefficients then either it does not exist, or there is a bug in Decompose. I do not recall seeing a bug reported in Decompose in a very long time (possibly 20 years). – Daniel Lichtblau Nov 1 '15 at 17:48
• The work by Mark Giesbrecht was just a bit after that of Alagar and Tranh. Jochim von zur Ga"then (Mark's doctoral advisor) published related work though I think a bit later. Kozen and Landau were somewhere in between. I'm not sure what was available when Decompose was actually written. – Daniel Lichtblau Nov 1 '15 at 21:38

I'll give a "proof of no such decomposition" response since it can shed light on the underlying method. Our method of implementation is slightly dated, going back to a paper from 30 years ago by Alagar and Thanh (it was not dated when originally coded, I believe by Dan Grayson). Among other things, this will show that the putative inner polynomial was a close guess to the only one that could work, up to a nonunique normalization.

First note that if $p(x)=f(g(x))$ then by the chain rule we have $p'(x)=f'(g(x))g'(x)$. So we factor $p$ to get candidates for $g'$. In this case there is but one such.

poly = 19456*z^8 + 98368*z^7 + 195820*z^6 + 140530*z^5 + 60493*z^4 +
14413*z^3 + 2149*z^2 + 196*z + 16;
derivpoly = D[poly, z];
Factor[derivpoly]

(* Out[173]= (1 + 8 z) (196 + 2730 z + 21399 z^2 + 70780 z^3 +
136410 z^4 + 83640 z^5 + 19456 z^6) *)


The degree of the second factor rules it out so if there is a nontrivial decomposition then we must have $g'(z)=1+8z$. We integrate to get $g(z)$, setting the undetermined constant term to $0$ (this is fine because decomposition does not determine the constant for the last polynomial).

gpoly = Integrate[1 + 8 z, z]

(* Out[174]= z + 4 z^2 *)


This agrees with the proposed polynomial up to constant term and scale, and as noted that latter is not uniquely determined anyway.

Now we form $f(z)$ with undetermined coefficients and attempt to solve for them. Clearly it must have degree 4.

fdeg = 4;
fpoly = Array[f, fdeg + 1, 0].z^Range[0, fdeg]

(* f[0] + z f[1] + z^2 f[2] + z^3 f[3] + z^4 f[4] *)


Expand with $z$ replaced by $g(z)$ to get the full polynomial in terms of the undetermined coefficients.

newpoly = Expand[fpoly /. z -> gpoly]

(* f[0] + z f[1] + 4 z^2 f[1] + z^2 f[2] + 8 z^3 f[2] + 16 z^4 f[2] +
z^3 f[3] + 12 z^4 f[3] + 48 z^5 f[3] + 64 z^6 f[3] + z^4 f[4] +
16 z^5 f[4] + 96 z^6 f[4] + 256 z^7 f[4] + 256 z^8 f[4] *)


Now we just equate coefficients to those of the input, and solve for the missing coefficients. This is now just a matter of linear algebra.

coeffs = CoefficientList[newpoly - poly, z]

(* Out[208]= {-16 + f[0], -196 + f[1], -2149 + 4 f[1] + f[2], -14413 +
8 f[2] + f[3], -60493 + 16 f[2] + 12 f[3] + f[4], -140530 +
48 f[3] + 16 f[4], -195820 + 64 f[3] + 96 f[4], -98368 +
256 f[4], -19456 + 256 f[4]} *)

soln = Solve[coeffs == 0]

(* Out[209]= {} *)


(This could have simply been done as SolveAlways[newpoly == poly, z] but I wanted to show a more complete set of steps.)

One can work with a subsystem that is not overdetermined to get a polynomial for which the decomposition does in fact work:

soln = Solve[coeffs[[1 ;; 5]] == 0]

newpoly /. soln[[1]]

(*
Out[210]= {{f[0] -> 16, f[1] -> 196, f[2] -> 1365, f[3] -> 3493, f[4] -> -3263}}

Out[211]=
16 + 196 z + 2149 z^2 + 14413 z^3 + 60493 z^4 + 115456 z^5 - 89696 z^6 -
835328 z^7 - 835328 z^8
*)