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

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.

• Total@Flatten@Last@RootIntervals[poly] will give you the number of real roots. – Virgil May 13 '15 at 16:40
• How about using Roots or NRoots? – Stefan R May 13 '15 at 16:40
• It can be done using quantors as well, but something strange happens then. See my question. – Fred Simons May 14 '15 at 6:28

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

• In your last line, you could use Exponent[poly, x] instead. – Chip Hurst Jun 12 '15 at 21:59
• Thanks! I changed the answer. I didn't know about this function. – David Zwicker Jun 13 '15 at 0:43

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


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

• FindInstance will not tell you if you have a repeated root. – george2079 Oct 21 '15 at 17:46

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

• Using ToRules@Reduce[poly == 0, x, Reals]  is faster. – Virgil May 13 '15 at 16:48

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

• The other thing that may be useful is that Root's ordering appears to put the real roots in order of increasing magnitude, and the complex roots in order of increasing magnitude of the real part. I seem to recall seeing this in the documentation somewhere, but I can't find it today. – John Doty May 14 '15 at 4:25

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