Can Mathematica tell me if a polynomial has all real roots?

Question

I am trying on this polynomial,

poly=331776*x^36 - 11943936*x^34 + 195747840*x^32 - 1932263424*x^30 +
 12809871360*x^28 - 60216016896*x^26 + 206610186240*x^24 -
 524928024576*x^22 + 991718940672*x^20 - 1386996203520*x^18 +
 1415900528640*x^16 - 1026732589056*x^14 + 505483296768*x^12 -
 158084628480*x^10 + 28283240448*x^8 - 2378170368*x^6;

---

• I was trying "Reduce[ (the polynomial) == 0, x, Reals] == Reduce [ (the same polynomial) == 0, x]" but this doesn't work

• Does Mathematica have its own implementation of the Sturm sequence algorithm to test for existence of real roots in a (finite?) interval?

• And if a polynomial has all real roots then I just want to know the value of the root with the highest magnitude.

Answer 1

There is a specific function called CountRoots, which does exactly what you want:

CountRoots[poly, x]
(* 12 *)

You can then compare this number to the number of roots given by the length of the coefficient list:

Exponent[poly, x]
(* 36 *)

Answer 2

There is probably a better way to do this, but one can find the degree of a polynomial via the length of the CoefficientList, and one can get the number of real roots by adding up the multiplicities given by RootIntervals:

poly = 331776*x^36 - 11943936*x^34 + 195747840*x^32 - 1932263424*x^30 
     + 12809871360*x^28 - 60216016896*x^26 + 206610186240*x^24
     - 524928024576*x^22 + 991718940672*x^20 - 1386996203520*x^18 
     + 1415900528640*x^16 - 1026732589056*x^14 + 505483296768*x^12 
     - 158084628480*x^10 + 28283240448*x^8 - 2378170368*x^6;

degree = Length@Rest@CoefficientList[poly, x]
realroots = Total@Flatten@Last@RootIntervals[poly]
degree == realroots

36
12
False

So, in this case there are complex roots.

To get the position of the maximum real root, you can use Reduce (or Roots, or Solve):

Max[x /. List@ToRules@Reduce[poly == 0, x, Reals]]

2

Answer 3

You can also use FindInstance, specifying the domain over all reals

pol = 331776*x^36 - 11943936*x^34 + 195747840*x^32 - 
  1932263424*x^30 + 12809871360*x^28 - 60216016896*x^26 + 
  206610186240*x^24 - 524928024576*x^22 + 991718940672*x^20 - 
  1386996203520*x^18 + 1415900528640*x^16 - 1026732589056*x^14 + 
  505483296768*x^12 - 158084628480*x^10 + 28283240448*x^8 - 
  2378170368*x^6

FindInstance[pol == 0, x, Reals, Exponent[pol, x]]

(*{{x -> 0},{x -> -2},{x -> 2}}*)

or within a particular subset of the reals

FindInstance[pol == 0 && 0 <= x < 3, x, Reals, Exponent[pol, x]]

(*{{x -> 0},{x -> 2}}*)

The maximum root value can be determined straightforwardly by

Max@FindInstance[pol == 0, x, Reals, Exponent[pol, x]][[;; , ;; , 2]]

(* 2 *)

Answer 4

I think that you can use RootIntervals to do what you want:

poly = 331776*x^36 - 11943936*x^34 + 195747840*x^32 - 1932263424*x^30 + 
  12809871360*x^28 - 60216016896*x^26 + 206610186240*x^24 - 
  524928024576*x^22 + 991718940672*x^20 - 1386996203520*x^18 + 
  1415900528640*x^16 - 1026732589056*x^14 + 505483296768*x^12 - 
  158084628480*x^10 + 28283240448*x^8 - 2378170368*x^6;
foo = RootIntervals[poly]

(* {{{-(143/64), -(1/512)}, {0, 0}, {1/512, 143/64}}, {{1, 1, 1}, {1, 1, 
   1, 1, 1, 1}, {1, 1, 1}}} *)

Solve[{poly == 0, Element[x, Interval[Last[foo[[1]]]]]}, x]

(* {{x -> {2}}, {x -> {2}}, {x -> {2}}} *)

Answer 5

I'd use Root. First, put your polynomial in pure functional form:

fpoly = Function[Evaluate[poly /. x -> #]]

Now, note that Root rigorously orders roots, putting real roots before complex ones. This is a 36th order polynomial, so look at the 36th root:

Im[Root[fpoly, 36]] == 0 (* False *)

Answer 6

First, let's simplify the problem slightly by removing the six zeroes at the origin, which we know to be real anyway:

p = -2378170368 + 28283240448 x^2 - 158084628480 x^4 + 505483296768 x^6 -
    1026732589056 x^8 + 1415900528640 x^10 - 1386996203520 x^12 +
    991718940672 x^14 - 524928024576 x^16 + 206610186240 x^18 -
    60216016896 x^20 + 12809871360 x^22 - 1932263424 x^24 +
    195747840 x^26 - 11943936 x^28 + 331776 x^30;

Here is a method that combines features from polynomial and linear algebra. The method, based on the work of Miroslav Fiedler and Gerhard Schmeisser, constructs (or attempts to construct) a tridiagonal companion matrix from a polynomial (i.e., a tridiagonal matrix whose characteristic polynomial is the given polynomial), using a modified Euclidean algorithm. Even if the eigenvalues of this matrix are not needed, the matrix still has diagnostic value for checking if a given polynomial has all its roots real. Here's how to build the diagonal and one of the off-diagonals:

n = Exponent[p, x];
p0 = p/Coefficient[p, x, n]; p1 = D[p0/n, x];
{d, e} = MapAt[Most, Transpose[
               Table[{q, r} = PolynomialQuotientRemainder[p0, p1, x];
                     s = If[TrueQ[r == 0], 0,
                            Coefficient[-r, x, Exponent[r, x]]];
                     If[k < n, p0 = p1;
                        p1 = If[TrueQ[r == 0],
                                (#/Coefficient[#, x, Exponent[#, x]]) &[D[p1, x]],
                                 -r/s]];
                     {-Coefficient[q, x, 0], s}, {k, n}]], 2]

(I've omitted the output, as the result has fractions with a lot of digits.)

Checking if all the roots are real is as easy as

And @@ NonNegative[e]
   False

(This previous answer also used this method to check for real roots.)