Rewrite a real polynomial in real (but only linear and quadratic) factors

According to my calculus book: "Every real polynomial can be factored into a product of real (possibly repeated) linear factors and real (also possibly repeated) quadratic factors having no real zeros."

Question: how to do this in Mathematica?

For example: $f(x) = x^5 - 3 x^4 + x^2 - 11 x + 6$

Factor[x^5 - 3 x^4 + x^2 - 11 x + 6]


This doesn't give the desired result (btw, why not?)

• I'm guessing that it's because it tries to factor the expression over the integers (or perhaps over the rationals) since you have the expression in exact form. Since that's not guaranteed (some of the roots will be irrational), it won't in general factor the polynomial. For inexact coefficients, it might do some factoring, but in general it would have to find the roots first. From the documentation: Factor factors a polynomial over the integers (by default). – march Oct 28 '15 at 23:23
• Note that the Fundamental Theorem of Algebra is merely an existence theorem, which implies that every real polynomial has a factorization into linear and quadratic real factors. However, since Abel & Galois it's been known that there is no analog of the quadratic formula for polynomials of degree > 4. So while it may be possible to explicitly factor some polynomials of degree 5, this is impossible for others. – murray Oct 30 '15 at 16:46

As noted, one sticking point is that polynomial roots in general cannot be represented in Mathematica by anything other than Root[] objects. Nevertheless, it is possible to do a few manipulations to factorize a real polynomial into its linear and quadratic factors.

decompose[poly_, x_] /; PolynomialQ[poly, x] := Module[{gr, rts},
rts = x /. Solve[poly == 0, x];
gr = GatherBy[rts, {Composition[Unitize, Im],
Composition[RootReduce, Re]}];
Coefficient[poly, x, Exponent[poly, x]]
Apply[Times, Factor[Collect[Times @@ #, x, RootReduce]] & /@
(x - Apply[Join, gr])]]


First, an easy case:

decompose[200 - 448 x + 204 x^2 + 35 x^3 - 89 x^4 + 46 x^5 - 9 x^6 + x^7, x]
(-2 + x) (1/2 (1 - Sqrt[5]) + x) (1/2 (1 + Sqrt[5]) + x)
(25 - 6 x + x^2) (4 - 2 x + x^2)


Here's the OP's example:

decompose[6 - 11 x + x^2 - 3 x^4 + x^5, x]
(x + Root[-6 - 11 #1 - #1^2 + 3 #1^4 + #1^5 &, 1])
(x + Root[-6 - 11 #1 - #1^2 + 3 #1^4 + #1^5 &, 2])
(x + Root[-6 - 11 #1 - #1^2 + 3 #1^4 + #1^5 &, 3])
(x^2 + Root[1296 - 216 #1 + 648 #1^2 + 1566 #1^3 - 1457 #1^4 + 270 #1^5 -
244 #1^6 + 80 #1^7 + 8 #1^8 + #1^10 &, 4] +
x Root[-718 - 1115 #1 + 116 #1^2 + 820 #1^3 + 539 #1^4 + 186 #1^5 + 102 #1^6 +
107 #1^7 + 54 #1^8 + 12 #1^9 + #1^10 &, 3])


Complicated, yes. But N[] reveals that the Root[]s are mere algebraic numbers after all:

N[%]
(-3.1838 + x) (-0.552473 + x) (1.53612 + x) (2.2206 - 0.799845 x + x^2)

• This doesn't work with x^10 by the way, why? – GambitSquared Oct 29 '15 at 10:38
• What output were you expecting? – J. M. will be back soon Oct 29 '15 at 11:47

You can find the roots using Solve:

sol = Solve[x^5 - 3 x^4 + x^2 - 11 x + 6 == 0, x] // N

{{x -> -1.53612}, {x -> 0.552473}, {x -> 3.1838},
{x -> 0.399923 - 1.4355 I}, {x -> 0.399923 + 1.4355 I}}


The linear factors are those corresponding to the real roots and the quadratic factors correspond to the imaginary roots. You can get the linear and quadratic factors with a bit more effort. This returns the linear and the quadratic terms separately:

roots = sol[[All, 1, 2]];
realRoots = x - Select[roots, Im[#] == 0 &]
complexRoots = x - Select[roots, Im[#] != 0 &];
Chop@ComplexExpand[Times @@ complexRoots]

{1.53612 + x, -0.552473 + x, -3.1838 + x}
2.2206 - 0.799845 x + x^2


In fact, you can save yourself a bit of work... march's comment is correct. All you need to do is to make one of the coefficients real (rather than integer)

Factor[x^5 - 3. x^4 + x^2 - 11 x + 6]

1. (-3.1838 + 1. x) (-0.552473 + 1. x) (1.53612 + 1. x)
(2.2206 -0.799845 x + 1. x^2)


This gathers roots into complex conjugate pairs after solving the system numerically (since factoring a fifth order polynomial analytically is tough):

In[1]:= poly = x^5 - 3 x^4 + x^2 - 11 x + 6;

In[2]:= Chop[Times @@ Flatten[Expand@*Times @@@ (GatherBy[
NSolve[poly],
{Re[x] /. #, Abs[Im[x]] /. #} &] /. Rule -> Subtract)]]
Out[2]= (-3.1838 + x) (-0.552473 + x) (1.53612 + x) (2.2206 - 0.799845 x + x^2)


A word of warning: this won't quite do want you want in the case of repeated roots, but it's possible to generalize the approach to work there, too.

Times@@Map[x-#&, N[x/.{ToRules[Reduce[x^5-3 x^4+x^2-11 x+6==0, x]]}]]


gives the result in factored form and not just as a list of roots or rules

(-3.1838+x)(-0.552473+x)(-0.399923-1.4355 I+x)(-0.399923+1.4355 I+x)(1.53612+x)


In general, factoring polynomials of order greater than four is tricky, but Mathematica can do it in some cases:

Factor[28.35 + 67.545 x - 109.935 x^2 - 66.11 x^3 + 3.1 x^4 + x^5]


(* 1. (-7.5 + 1. x) (-0.7 + 1. x) (0.3 + 1. x) (2. + 1. x) (9. + 1. x) *)

or

NSolve[x^5 - 3 x^4 + x^2 - 11 x + 6 == 0, x]


{{x -> -1.53612}, {x -> 0.399923 - 1.4355 I}, {x -> 0.399923 + 1.4355 I}, {x -> 0.552473}, {x -> 3.1838}}

Thank you very much for all your ideas. I have combined them in this solution:

poly = x^5 - 3 x^4 + x^2 - 11 x + 6
roots = (x /. Solve[poly == 0, x]) // N
gr = x - Join[GatherBy[Cases[roots, n_ /; Im[n] != 0], {Re[#], Abs@Im[#]} &],DeleteCases[roots, n_ /; Im[n] != 0]] /. n_ /; Im[n] == 0 -> {n}
Times @@ Chop[Expand[# /. List -> Times & /@ gr]]


Output:

(-3.1838 + x) (-0.552473 + x) (1.53612 + x) (2.2206 - 0.799845 x + x^2)


It also works with real repeated roots, such as $f(x)=-8 + 8 x + 2 x^2 - 3 x^4 + x^5$:

poly = -8 + 8 x + 2 x^2 - 3 x^4 + x^5
roots = (x /. Solve[poly == 0, x]) // N
gr = x - Join[GatherBy[Cases[roots, n_ /; Im[n] != 0], {Re[#], Abs@Im[#]} &],DeleteCases[roots, n_ /; Im[n] != 0]] /. n_ /; Im[n] == 0 -> {n}
Times @@ Chop[Expand[# /. List -> Times & /@ gr]]


Output:

(-2. + x)^2 (-1. + x) (2. + 2. x + x^2)


Unfortunately it doesn't work with complex multiplicities.