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
- Define a generic polynomial with 1 as the leading coefficient
pX[a_,n_]:=x^n+Sum[a[i]x^i,{i,0,n-1}]
- 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]
]
- 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
- 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}
- 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}
- The coefficients in
rule[1] and rule[2] have been selected as to
match the original equation. However, other choices are possible.
