## Factor a polynomial Root into Roots of smallest possible degree

15

4

Suppose I have a polynomial Root representing an algebraic number. I want to represent it as a product of several polynomial Roots (if possible) such that the largest degree among the factors is as small as possible. For example, Root[-1 - #1 + 3 #1^3 - #1^4 + #1^5 - 3 #1^6 + 2 #1^7 + #1^9 &, 1] has to be represented as Root[1 + # + #^3 &, 1] Root[1 - # + #^3 &, 1].

The second problem is the same, but replacing "product" with "sum". For example, Root[8 - 4 #1 + 24 #1^2 - 15 #1^3 + 3 #1^5 + 6 #1^6 + #1^9 &, 1] has to be represented as Root[1 + # + #^3 &, 1] + Root[1 - # + #^3 &, 1].

How can I implement solutions to these problems in Mathematica?

5I'm sure you are aware of this but I'll say it anyway: this is not an easy problem. – Daniel Lichtblau – 2016-05-09T14:14:15.037

Would this perhaps be better asked on [Math.SE]? – Mr.Wizard – 2016-07-08T06:31:44.520

Recently we worked on a Metropolis-Hasting algorithm here in permutation-98 space. Consider two polynomials in rational coefficients of degree 3. That's rational-8 space. Is it inconceivable to adapt this algorithm to search this space and converge to a solution? – Dominic – 2019-08-07T20:49:00.127

1@Dominic for Metropolis-Hastings you need to define a reasonably continuous quality function that you can try to maximize or minimize. Do you have any ideas for how to define it? – Roman – 2019-08-07T22:57:37.927

You might get some inspiration from the Fiedler factorization of companion matrices. – yarchik – 2019-08-08T07:30:25.013

2

Just to state the obvious, going in opposite direction can be done with the help of direct products of companion matrices. Here, the problem is to write a companion matrix as a direct product of companion matrices. https://math.stackexchange.com/questions/331017/enlightening-proof-that-the-algebraic-numbers-form-a-field?noredirect=1&lq=1

– yarchik – 2019-08-08T07:35:34.090

@yarchik - you can go in the opposite direction using RootReduce

– Bob Hanlon – 2019-08-09T03:34:43.240

12

# 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}]

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]
]

1. 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

1. 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}

1. 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}

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

1Very nice (+ upvote, of course). I missed a beat on this one-- it can similarly be done using resultants and solving. Basically equivalent to your approach, if you drill deep enough. – Daniel Lichtblau – 2019-08-08T14:54:54.630

Thanks! In general, I would like an approach that would work not only for cubic root factors, but would be able to factor or split into a sum roots of any order (if such factorization exists). – Vladimir Reshetnikov – 2019-08-08T16:22:46.197

@VladimirReshetnikov Nice, can you provide some challenging case? I would be interested to try. – yarchik – 2019-08-08T16:28:22.537

1E.g. Root[-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 &, 1] can be represented as Root[-656 - 150 # - 80 #^2 - 20 #^3 + #^5 &, 1] * Root[-1 + # + #^3 &, 1] and ... – Vladimir Reshetnikov – 2019-08-08T18:38:29.930

1... actually, can be expressed in radicals: $\frac{\left(\sqrt[3]{6} \beta ^2-2\ 3^{2/3}\right) \left(\alpha \left(\sqrt[10]{2} \gamma ^4+\left(\sqrt[5]{2}+2^{7/10}\right) \gamma ^2+2 \gamma +2^{3/10} \left(5+3 \sqrt{2}\right)\right)+\left(-68-35 \sqrt{2}+2 \sqrt{4100-1918 \sqrt{2}}\right) \gamma-2\ 2^{3/10} \left(68+35 \sqrt{2}\right)\right)}{3\ 2^{17/30} \alpha \beta \gamma ^3},,\small\text{where }\alpha =\sqrt{2050-959 \sqrt{2}},,\beta =\sqrt[3]{9+\sqrt{93}},,\gamma =\sqrt[5]{102+72 \sqrt{2}-\sqrt{20714+14647 \sqrt{2}}}$ – Vladimir Reshetnikov – 2019-08-08T18:39:09.717

@VladimirReshetnikov See edits.. – yarchik – 2019-08-08T19:43:34.827

So, if we do not know a form of the solution beforehand, we can iterate through all combinations of lower powers and try them, right? – Vladimir Reshetnikov – 2019-08-10T03:01:21.890

@VladimirReshetnikov That would be my suggestion too. One may also consider >2 factors. One interesting example of this kind is the Ramanujan expression for the singular modulus $k_{210}$. It satisfies a 96-degree equation and can be written as a product of 22 factors. – yarchik – 2019-08-10T14:13:00.787

I was working on simplifying Root[-5 + 20 # - 125 #^2 + 305 #^3 - 275 #^4 + 65 #^5 + 40 #^6 - 25 #^7 + #^8 &, 1] that arises in evaluation of InverseBetaRegularized[3/5, 1/3, 1/3]. I could not directly represent it as a sum or product of roots of smaller degree, but using ad hoc methods found that is can be written as Root[-5 + 190 # + 60 #^2 - 200 #^3 + 16 #^4 &, 1] - Sqrt[Root[3645 + 72900 # + 41040 #^2 - 31680 #^3 + 256 #^4 &, 3]] and, hence, can be represented in radicals. Can your method be adapted to find representations like this? – Vladimir Reshetnikov – 2019-08-11T02:38:41.680

@VladimirReshetnikov I do not know, maybe it is not possible because the functional form is different. But what is the incentive of writing it as a sum of 2 roots? The equation is solvable (with extension $\sqrt{5}$) and the solution can already be written as 1 quartic root. – yarchik – 2019-08-11T13:14:26.330

Is there a systematic way to find extension to factor a polynomial, or $\sqrt5$ was just a guess? – Vladimir Reshetnikov – 2019-08-11T18:34:00.540

How would you approach to expressing Root[4096 + 49152 # - 189696 #^2 + 123136 #^3 - 5808 #^4 - 564 #^5 + #^6 &, 3] in radicals? This root arises in an evaluation of InverseBetaRegularized[4/7, 1/3, 1/2] – Vladimir Reshetnikov – 2019-08-11T19:04:58.297

1

@VladimirReshetnikov Try $\sqrt{21}$ as an extension. Yes, these were my guesses, more systematic way is through the computation of the Galois group of the polynomial. There are some MA packages available (http://library.wolfram.com/infocenter/TechNotes/158/) and some nice posts on math.stackexchange: https://math.stackexchange.com/questions/45893/how-to-find-the-galois-group-of-a-polynomial .

– yarchik – 2019-08-11T19:54:54.877

1

Thanks! By the way, I tried the package Solving the Quintic with Mathematica you linked, and it worked unreliable for me, failing to found existing solutions sometimes. Here is my implementation of a method from Daniel Lazard's paper Solving Quintics by Radicals: https://mathematica.stackexchange.com/a/134551/7288. I corrected a few typos found in the paper to make it work, and so far have not found any bugs in it.

Suppose we have a root $P$ of degree $n$ that has a (yet unknown to us) representation $P = Q + R,S$, where $Q,,R,,S,$ are some roots we want to find, all of degrees smaller than $n$. Suppose that the product $R,S$ is of degree $n$, so we cannot, as a first step, split $P$ into a sum of roots of smaller degrees. Is there a way find $Q,,R,,S,$ in such a case? – Vladimir Reshetnikov – 2019-08-13T20:19:23.030

@VladimirReshetnikov See edits. For illustrative purposes I still do the decomposition sequentially. This is faster. However, provided sufficient computer power and memory, one can try to do it in one step. But in general, you probably need some general mathematical statement whether a multiplicative or an additive decomposition is always possible. This is an interesting math question, why not asking it on the mathoverflow. – yarchik – 2019-08-14T08:06:10.753

@VladimirReshetnikov Thank you too and good luck with your problem solving! That was an interesting project. – yarchik – 2019-08-14T16:09:24.797

9

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.

Nice trick with 2 Resultants (+). I thought, you would be using 1 resultant. – yarchik – 2019-08-09T11:12:10.743

@yarchik It can be done with one but then one of the polynomials needs to have a variable substitution, something like x-->y/x, and I can never remember exactly what it is. – Daniel Lichtblau – 2019-08-09T15:00:20.500