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

Question

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?)

Answer 1

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)

Answer 2

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)

Answer 3

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.

Answer 4

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)

Answer 5

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

Answer 6

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.