Factor a polynomial Root into Roots of smallest possible degree

A constructive approach

The problem can be solved if the form of the solution is given.

Define the two factors using a hint (that these should be cubic equations) in the original post

y1 = x^3 - p1  x + q1
y2 = x^3 - p2  x + q2

Build a companion matrix of the polynomial $p(x)$

CompanionMatrix[p_,x_]:=Module[{n,w=CoefficientList[p,x]},w=-w/Last[w];
     n=Length[w]-1; SparseArray[{{i_,n}:>w[[i]],{i_,j_}/;i==j+1->1},{n,n}]]

• The roots of a polynomial equation $p_A(x)=0$ are given by the eigenvalues of its companion matrix $A$.
• If $a$ is a root of $p_A(x)$ (with companion matrix $A$) and $b$ is a root of $p_B(x)$ (with companion matrix $B$) then $a b$ is an eigenvalue of $A\otimes B$ and $a+b$ is an eigenvalue of $A\otimes I+I\otimes B$, where $\otimes$ is the direct product of matrices.

Let us focus on the product case and, therefore, determine the characteristic equation of the direct product of companion matrices

(a1=CompanionMatrix[y1,x])//MatrixForm
(a2=CompanionMatrix[y2,x])//MatrixForm
 y3=CharacteristicPolynomial[KroneckerProduct[a1,a2],x]
 y4=y3 Sign[CoefficientList[y3,x]//Last]

The sum is treated similarly.

We find out that the resulting characteristic polynomial is

 y4=-q1^3 q2^3 + p1 p2 q1^2 q2^2 x + (-p2^3 q1^2 - p1^3 q2^2 + 3 q1^2 q2^2) x^3 + p1 p2 q1 q2 x^4 + p1^2 p2^2 x^5 - 3 q1 q2 x^6 - 2 p1 p2 x^7 + x^9

Notice the relatively simple form in the case of trinomial equations. Now, your polynomial is

z = -1 - #1 + 3 #1^3 - #1^4 + #1^5 - 3 #1^6 + 2 #1^7 + #1^9 &@x

We demand that the two polynomials (y3 and z) are equal for any value of x

r = SolveAlways[z == y4, {x}]

Several solutions are obtained. Take, for instance, the first one

{y1, y2} /. r[[1]]

(*{1/q2 + x/q2^(2/3) + x^3, q2 - q2^(2/3) x + x^3}*)

Setting q2=1 we obtain the OP result.

Comment on the method

MA has a nice RootApproximant function, which operates by virtue of the LLL algorithm. One may be tempted to follow this route and implement some kind of this experimental mathematics approach. In contrast, the presented solution is fully constructive (in fact, algebraic) and does not require arbitrary precision computations.

Answer to the 1st challenge question

y1=x^5+ m x^3+n x^2+p1 x+q1;
y2=x^3+ p2 x+q2;
a1=CompanionMatrix[y1,x];
a2=CompanionMatrix[y2,x];
y3=CharacteristicPolynomial[KroneckerProduct[a1,a2],x];
y4=y3 Sign[CoefficientList[y3,x]//Last];
z=-282300416-64550400 #-34426880 #^2-14185880 #^3+8564800 #^4+4231216 #^5-972800 #^6-367820 #^7+27360 #^8+2600 #^9+1680 #^10+100 #^11-240 #^12+40 #^13+#^15&@x;
r=SolveAlways[y4==z,{x}]
({y1,y2}/.r[[1]])/.q2->1
Out[1]= {656-150 x+80 x^2-20 x^3+x^5,1+x+x^3}

• Notice that in the characteristic equation the highest order term can be negative, whereas I assume that in the given equation z it is always 1. Therefore, in order to use SolveAlways y3 is multiplied by the sign.

• From the shown solutions, one can see that there is some arbitrariness in the results. But this is, of course, expected. As for the shown solution in radicals, one probably needs to guess the field extension.

• One can create a Module as to fully automate the derivation. But is there a pressing need?

Answer to the 2st challenge question

Here I demonstrate that the method can be used to write a root $P$ in terms of three roots of lower-order polynomials as $P=Q+R S$. Additionally, the computation is done in a more structured form

1. Define a generic polynomial with 1 as the leading coefficient

---

pX[a_,n_]:=x^n+Sum[a[i]x^i,{i,0,n-1}] 

2. Define two modules that split a root in terms of a sum or a product

---

pSum[pA_,pB_,y_]:=Module[{mA,mB,mAB,nA,nB,p},
  nA=Exponent[pA,y];
  nB=Exponent[pB,y];
  mA=CompanionMatrix[pA,y];
  mB=CompanionMatrix[pB,y];
  mAB=KroneckerProduct[mA,IdentityMatrix[nB]]+KroneckerProduct[IdentityMatrix[nA],mB];
  p=CharacteristicPolynomial[mAB,y];
  p Sign[CoefficientList[p,y]//Last]
]

pProduct[pA_,pB_,y_]:=Module[{mA,mB,mAB,p},
  mA=CompanionMatrix[pA,y];
  mB=CompanionMatrix[pB,y];
  mAB=KroneckerProduct[mA,mB];
  p=CharacteristicPolynomial[mAB,y];
  p Sign[CoefficientList[p,y]//Last]
] 

3. Construct a working example for polynomials of the degrees 3, 2, 2 respectively.

---

{nQ,nR,nS}={3,2,2};
nT=nR nS;    
pP[0]=(RootReduce[Root[#^3+#+3&,1]+Root[#^2+#+5&,1]Root[#^2+#+7&,1]]//First)@x

Out[1]= 1984051825-172403780 x-281288553 x^2+14148329 x^3+17544721 x^4-310509 x^5-619703 x^6-4623 x^7+13443 x^8+244 x^9-167 x^10-3 x^11+x^12

4. Do the first part, namely, split the root of a 12th order equation into a sum of 3rd and 4th order roots

---

(sol[1]=SolveAlways[pP[0]==pSum[pX[q,nQ],pX[t,nT],x],x])//Transpose//TableForm
rule[1]=q[2]->0;
sol[1,1]=First[sol[1]]/.rule[1];
AppendTo[sol[1,1],rule[1]];
{pX[q,nQ],pX[t,nT]}/.sol[1,1]

Out[2]= {3+x+x^3,1225-35 x-58 x^2-x^3+x^4}

5. Do the second part, split the 4th order root into a product of 2 roots of quadratic equations

---

(sol[2]=SolveAlways[(pX[t,nT]/.sol[1,1])==pProduct[pX[r,nR],pX[s,nS],x],x])//Transpose//TableForm
rule[2]=s[1]->1;
sol[2,1]=First[sol[2]]/.rule[2];
AppendTo[sol[2,1],rule[2]];
{pX[r,nR],pX[s,nS]}/.sol[2,1]

Out[3]= {5+x+x^2,7+x+x^2}

6. The coefficients in rule[1] and rule[2] have been selected as to match the original equation. However, other choices are possible.

I'll show the resultant formulation for the degree 15 example.

The polynomial in question:

poly15 = (-282300416 - 64550400 # - 34426880 #^2 - 14185880 #^3 + 
      8564800 #^4 + 4231216 #^5 - 972800 #^6 - 367820 #^7 + 
      27360 #^8 + 2600 #^9 + 1680 #^10 + 100 #^11 - 240 #^12 + 
      40 #^13 + #^15) &[z];

Define monic polynomials of degree 3 and 5 with symbolic coefficients.

poly3[x_] := Array[a, 3, 0].x^Range[0, 2] + x^3
poly5[x_] := Array[b, 5, 0].x^Range[0, 4] + x^5

Here is the double-resultant way to get a polynomial whose roots are products of the roots of poly3 and poly5.

resxy = Resultant[Resultant[z - x*y, poly3[y], y], poly5[x], x];

Equate coefficients with this and the degree 15 polynomial. Since we can scale the two roots, force a constant term to be unity.

coeffs = CoefficientList[resxy - poly15, z] /. {a[0] -> 1}

Now solve:

solns=Solve[coeffs == 0]

(* Out[285]= {{a[1] -> 1, a[2] -> 0, b[0] -> 656, b[1] -> -150, 
  b[2] -> 80, b[3] -> -20, b[4] -> 0}, {a[1] -> -(-1)^(1/3), 
  a[2] -> 0, b[0] -> 656 (-1 + (-1)^(1/3)), b[1] -> 150 (-1)^(1/3), 
  b[2] -> 80, b[3] -> 20 (1 - (-1)^(1/3)), 
  b[4] -> 0}, {a[1] -> (-1)^(2/3), a[2] -> 0, 
  b[0] -> 656 (-1 - (-1)^(2/3)), b[1] -> -150 (-1)^(2/3), b[2] -> 80, 
  b[3] -> 20 (1 + (-1)^(2/3)), b[4] -> 0}} *)

Recover polynomials that give the root pairs. To get integer coefficients we use the first solution above.

{poly3[x], poly5[x]} /. solns[[1]]

(* Out[289]= {x + x^3 + a[0], 656 - 150 x + 80 x^2 - 20 x^3 + x^5} *)

Not so well explained as the approach by @yarchik (which I'd expect to see get the bounty), but at some level it is equivalent.