12  Roots of a polynomial

In this section we use the following add on packages:

using CalculusWithJuliaSquared
plotly()

CalculusWithJuliaSquared brings the Symbolics computer algebra system, Plots and Roots with it.

---

The roots of a polynomial are the values of \(x\) that when substituted into the expression yield \(0\). For example, the polynomial \(x^2 - x\) has two roots, \(0\) and \(1\). A simple graph verifies this:

f(x) = x^2 - x
plot(f, -2, 2)
plot!(zero, -2, 2)

The graph crosses the \(x\)-axis at both \(0\) and \(1\).

What is known about polynomial roots? Some simple questions might be:

• Will a polynomial always have a root?
• How many roots can there be?
• How large can the roots be?

We look at such questions here.

12.0.1 The factor theorem

We begin with a comment that ties together two concepts related to polynomials. It allows us to speak of roots or factors interchangeably:

The factor theorem relates the roots of a polynomial with its factors: \(r\) is a root of \(p\) if and only if \((x-r)\) is a factor of the polynomial \(p\).

Clearly, if \(p\) is factored as \(a(x-r_1) \cdot (x-r_2) \cdots (x - r_k)\) then each \(r_i\) is a root, as a product involving at least one \(0\) term will be \(0\). The other implication is a consequence of polynomial division.

12.0.2 Polynomial Division

Euclidean division of integers \(a, b\) uniquely writes \(a = b\cdot q + r\) where \(0 \leq r < |b|\). The quotient is \(q\) and the remainder \(r\). There is an analogy for polynomial division, where for two polynomial functions \(f(x)\) and \(g(x)\) it is possible to write

\[ f(x) = g(x) \cdot q(x) + r(x) \]

where the degree of \(r\) is less than the degree of \(g(x)\). The long-division algorithm can be used to find both \(q(x)\) and \(r(x)\).

For the special case of a linear factor where \(g(x) = x - c\), the remainder must be of degree \(0\) (a non-zero constant) or the \(0\) polynomial. The above simplifies to

\[ f(x) = (x-c) \cdot q(x) + r \]

From this, we see that \(f(c) = r\). Hence, when \(c\) is a root of \(f(x)\), then it must be that \(r=0\) and so, \((x-c)\) is a factor.

---

The division algorithm for the case of linear term, \((x-c)\), can be carried out by the synthetic division algorithm. This algorithm produces \(q(x)\) and \(r\), a.k.a \(f(c)\). The Wikipedia page describes the algorithm well.

The following is an example where \(f(x) = x^4 + 2x^2 + 5\) and \(g(x) = x-2\):

2 | 1 0 2  0  5
  |   2 4 12 24
  -------------
    1 2 6 12 29

The polynomial \(f(x)\) is coded in terms of its coefficients (\(a_n\), \(a_{n-1}\), \(\dots\), \(a_1\), \(a_0\)) and is written on the top row. The algorithm then proceeds from left to right. The number just produced on the bottom row is multiplied by \(c\) and placed under the coefficient of \(f(x)\). The values are then added to produce the next number. The sequence produced above is 1 2 6 12 29. The last value (29) is \(r=f(c)\), the others encode the coefficients of q(x), which for this problem is \(q(x)=x^3 + 2x^2 + 6x + 12\). That is, we have written:

\[ x^4 + 2x^2 + 5 = (x-2) \cdot (x^3 + 2x^2 + 6x + 12) + 29. \]

As \(r\) is not \(0\), we can say that \(2\) is not a root of \(f(x)\).

If we were to track down the computation that produced \(f(2) = 29\), we would have

\[ 5 + 2 \cdot (0 + 2 \cdot (2 + 2 \cdot (0 + (2 \cdot 1)))) \]

In terms of \(c\) and the coefficients \(a_0, a_1, a_2, a_3\), and \(a_4\) this is

\[ a_0 + c\cdot(a_1 + c\cdot(a_2 + c\cdot(a_3 + c\cdot a_4))), \]

The above pattern provides a means to compute \(f(c)\) and could easily be generalized for higher degree polynomials. This generalization is called Horner’s method. Horner’s method has the advantage of also being faster and more accurate when floating point issues are accounted for.

A simple implementation of Horner’s algorithm would look like this, if indexing were 0-based:

function horner(p, x)
    n = degree(p)
    Σ = p[n]
    for i in (n-1):-1:0
       Σ = Σ * x + p[i]
    end
    return(Σ)
end

horner (generic function with 1 method)

Recording the different values of Σ would recover the polynomial q.

Julia has a built-in method, evalpoly, to compute polynomial evaluations this way. To illustrate:

@variables x;

p = (1, 2, 3, 4, 5) # 1 + 2x + 3x^2 + 4x^3 + 5x^4
evalpoly(x, p)

\[ x \left( x \left( x \left( 5 x + 4 \right) + 3 \right) + 2 \right) + 1 \]

---

Symbolics can carry out polynomial long division.

This naive attempt to divide won’t “just work” though:

(x^4 + 2x^2 + 5) / (x-2)

\[ \frac{x^{4} + 2 x^{2} + 5}{x - 2} \]

Dividing with / builds a fraction and leaves it as it is. As written, there is no compelling reason to change the expression, though in our example we want the division done.

For this task, divrem is available:

quotient, remainder = divrem(x^4 + 2x^2 + 5, x - 2)

\[ \left( x^{3} + 2 x^{2} + 6 x + 12,\ 29 \right) \]

The answer is a tuple containing the quotient and remainder. The quotient alone can be found with div or ÷, and the remainder alone with poly_rem:

(x^4 + 2x^2 + 5) ÷ (x - 2), poly_rem(x^4 + 2x^2 + 5, x - 2)

\[ \left( x^{3} + 2 x^{2} + 6 x + 12,\ 29 \right) \]

NoteWhy not rem?

Julia’s generic function for a remainder is rem, but it will not do here. Symbolics already defines rem for symbolic expressions, and what it returns is not a polynomial remainder but the call itself, unevaluated:

rem(x^4 + 2x^2 + 5, x - 2)

\[ \operatorname{rem}\left( x^{4} + 2 x^{2} + 5, x - 2 \right) \]

That is why the remainder of a polynomial division has a name of its own, poly_rem.

As well, the partial_fractions function could be used for this task (SymPy calls it apart). This function computes the partial fraction decomposition of a ratio of polynomial functions. It needs to be told which symbol is the variable:

pf = partial_fractions((x^4 + 2x^2 + 5) / (x-2), x)

\[ x^{3} + 2 x^{2} + 6 x + 12 + \frac{29}{x - 2} \]

The function combine_fractions combines such terms over a common denominator, as an “inverse” to partial_fractions (SymPy calls it together):

combine_fractions(pf)

\[ \frac{x^{4} + 2 x^{2} + 5}{x - 2} \]

This isn’t so much of interest at the moment, but will be when techniques of integration are looked at.

12.0.3 The rational root theorem

Factoring polynomials to find roots is a task that most all readers here will recognize, and, perhaps, remember not so fondly. One helpful trick to find possible roots by hand is the rational root theorem: if a polynomial has integer coefficients with \(a_0 \neq 0\), then any rational root, \(p/q\), must have \(p\) dividing the constant \(a_0\) and \(q\) dividing the leading term \(a_n\).

To glimpse why, suppose we have a polynomial with a rational root and integer coefficients. With this in mind, a polynomial with identical roots may be written as \((qx -p)(a_{n-1}x^{n-1}+\cdots a_1 x + a_0)\), where each coefficient is an integer. Multiplying through, we get that the polynomial is \(qa_{n-1}x^n + \cdots + pa_0\). So \(q\) is a factor of the leading coefficient and \(p\) is a factor of the constant.

An immediate consequence is that if the polynomial with integer coefficients is monic, then any rational root must be an integer.

This gives a finite - though possibly large - set of values that can be checked to exhaust the possibility of a rational root. By hand this process can be tedious, though may be speeded up using synthetic division. This task is one of the mainstays of high school algebra where problems are chosen judiciously to avoid too many possibilities.

However, one of the great triumphs of computer algebra is the ability to factor polynomials with integer (or rational) coefficients over the rational numbers. This is typically done by first factoring over modular numbers (akin to those on a clock face) and has nothing to do with the rational root test.

factored_poly can quickly find such a factorization, even for quite large polynomials with rational or integer coefficients. (It hands the work to the Nemo computer algebra library, which CalculusWithJuliaSquared loads.)

For example, factoring \(p = 2x^4 + x^3 -19x^2 -9x +9\). This has possible rational roots of plus or minus \(1\) or \(2\) divided by \(1\), \(3\), or \(9\) - \(12\) possible answers for this modest question. By hand that can be a bit of work, but factored_poly does it without fuss. Like partial_fractions, it is told which symbol is the variable:

p = 2x^4 + x^3 - 19x^2 - 9x + 9
factored_poly(p, x)

\[ \left( x + 3 \right) \left( x + 1 \right) \left( 2 x - 1 \right) \left( x - 3 \right) \]

To work with the factors themselves, poly_factors returns them as a list:

poly_factors(p, x)

\[ \left[ \begin{array}{c} x + 3 \\ x + 1 \\ 2 x - 1 \\ x - 3 \\ \end{array} \right] \]

12.0.4 The fundamental theorem of algebra

There is a basic fact about the roots of a polynomial of degree \(n\). Before formally stating it, we consider the earlier observation that a polynomial of degree \(n\) for large values of \(x\) has a graph that looks like the leading term. However, except at \(0\), monomials do not cross the \(x\) axis, the roots must be the result of the interaction of lower order terms. Intuitively, since each term can contribute only one basic shape up or down, there can not be arbitrarily many roots. In fact, a consequence of the Fundamental Theorem of Algebra (Gauss) is:

A polynomial of degree \(n\) with real or complex coefficients has at most \(n\) real roots.

This statement can be proved with the factor theorem and the division algorithm.

In fact the fundamental theorem states that there are exactly \(n\) roots, though, in general, one must consider multiple roots and possible complex roots to get all \(n\). (Consider \(x^2\) to see why multiplicity must be accounted for and \(x^2 + 1\) to see why complex values may be necessary.)

WarningWarning

The special case of the \(0\) polynomial having no degree defined eliminates needing to exclude it, as it has infinitely many roots. Otherwise, the language would be qualified to have \(n \geq 0\).

12.1 Finding roots of a polynomial

Knowing that a certain number of roots exist and actually finding those roots are different matters. For the simplest cases (the linear case) with \(a_0 + a_1x\), we know by solving algebraically that the root is \(-a_0/a_1\). (We assume \(a_1\neq 0\).) Of course, when \(a_1 \neq 0\), the graph of the polynomial will be a line with some non-zero slope, so will cross the \(x\)-axis as the line and this axis are not parallel.

For the quadratic case, there is the famous quadratic formula (known since \(2000\) BC) to find the two roots guaranteed by the formula:

\[ \frac{-b \pm \sqrt{b^2 - 4ac}}{2a}. \]

The discriminant is defined as \(b^2 - 4ac\). When this is negative, the square root requires the concept of complex numbers to be defined, and the formula shows the two complex roots are conjugates. When the discriminant is \(0\), then the root has multiplicity two, e.g., the polynomial will factor as \(a_2(x-r)^2\). Finally, when the discriminant is positive, there will be two distinct, real roots. This figure shows the \(3\) cases, that are illustrated by \(x^2 -1\), \(x^2\) and \(x^2 + 1\):

plot(x^2 - 1,  -2, 2, legend=false)  # two roots
plot!(x^2, -2, 2)                           # one (double) root
plot!(x^2 + 1, -2, 2)                       # no real root
plot!(zero, -2, 2)

There are similar formulas for the cubic and quartic cases. (The cubic formula was known to Cardano in \(1545\), though through Tartaglia, and the quartic was solved by Ferrari, Cardano’s roommate.)

In general, there is no such formula using radicals for \(5\)th degree polynomials or higher, a proof first given by Ruffini in \(1803\) with improvement by Abel in \(1824\). Even though the fundamental theorem shows that any polynomial can be factored into linear and quadratic terms, there is no general method as to how. (It is the case that some such polynomials may be solvable by radicals, just not all of them.)

The factored_poly function factors over the rational numbers: every factor it finds has rational coefficients, so each linear factor it finds corresponds to a rational root. Irrational roots stay hidden inside a factor it leaves whole. There are alternatives.

Finding roots can also be done through the symbolic_solve function, a function which also has a more general usage, as it can solve other kinds of equations too. (It can also solve several polynomial equations at once, such as where a circle meets a line, but that needs one more package, Groebner, which these notes do not load yet.) Here we illustrate that symbolic_solve can easily handle quadratic expressions:

symbolic_solve(x^2 + 2x - 3 ~ 0, x)

\[ \left[ \begin{array}{c} -3 \\ 1 \\ \end{array} \right] \]

The answer is a vector of values that when substituted in for the free variable x produce \(0.\)

We use the ~ notation to define an equation to pass to symbolic_solve. This convention is not necessary here, as symbolic_solve will assume an expression passed to it is an equation set to 0, but is pedagogically useful. Equations do not have an equals sign, which is reserved for assignment. To solve a more complicated expression of the type \(f(x) = g(x),\) one can solve \(f(x) - g(x) = 0,\) or use f ~ g.

When the expression to solve has more than one free variable, the variable to solve for should be explicitly stated with a second argument. (The specification above is unnecessary.) For example, here we show that symbolic_solve is aware of the quadratic formula:

@variables a b c
symbolic_solve(a*x^2 + b*x + c ~ 0, x)

\[ \left[ \begin{array}{c} \frac{-b + \sqrt{b^{2} - 4 a c}}{2 a} \\ \frac{-b - \sqrt{b^{2} - 4 a c}}{2 a} \\ \end{array} \right] \]

12.1.1 Real variables, complex roots

Some computer algebra systems let a variable carry assumptions: SymPy, for one, can be told that b is real or that c is positive, and its solve then discards any solution that breaks the assumption. symbolic_solve works differently. The variables made by @variables are real by default, yet symbolic_solve does not use that to restrict its answers. Here b is real, and \(b^2 + 1 = 0\) has no real solution, but the two complex solutions come back anyway:

symbolic_solve(b^2 + 1 ~ 0, b)

\[ \left[ \begin{array}{c} i \\ -i \\ \end{array} \right] \]

There is also no way to declare a variable positive. Solving \(c + 1 = 0\) simply gives the negative answer:

symbolic_solve(c + 1 ~ 0, c)

\[ \left[ \begin{array}{c} -1 \\ \end{array} \right] \]

So symbolic_solve finds all the roots, real and complex, and keeping only the ones wanted is a separate step, taken afterwards. For the real roots, the numeric_roots function (described below) has an option to do this. The cubic \(x^3 - x^2 + x - 1\) has one real root and two complex ones:

symbolic_solve(x^3 - x^2 + x - 1 ~ 0, x)

\[ \left[ \begin{array}{c} 1 \\ i \\ -i \\ \end{array} \right] \]

and real_only = true keeps just the real one:

numeric_roots(x^3 - x^2 + x - 1, x; real_only = true)

1-element Vector{Float64}:
 1.0

12.1.2 Factoring with symbolic_solve

Previously, it was mentioned that factored_poly only factors over the rational numbers. However, symbolic_solve can be used to factor. Here is an example:

factored_poly(x^2 - 2, x)

\[ x^{2} - 2 \]

Nothing is found, as the roots are \(\pm \sqrt{2}\), irrational numbers.

rts = symbolic_solve(x^2 - 2 ~ 0, x)
prod(x-r for r in rts)

\[ \left( x - \sqrt{2} \right) \left( x + \sqrt{2} \right) \]

Solving cubics and quartics can be done exactly using radicals. For example, here we see the solutions to a quartic equation can be quite involved, yet still explicit. (As above, complex-valued solutions will be found too: \(x^4 - 2x - 1\) has two real roots and two complex ones.)

symbolic_solve(x^4 - 2x - 1 ~ 0, x)

\[ \left[ \begin{array}{c} \frac{\sqrt{-2 \sqrt[3]{\frac{1}{4} + \frac{\sqrt{\frac{43}{3}}}{12}} - 2 \sqrt[3]{\frac{1}{4} - \frac{\sqrt{\frac{43}{3}}}{12}} - \frac{4}{\sqrt{2 \sqrt[3]{\frac{1}{4} + \frac{\sqrt{\frac{43}{3}}}{12}} + 2 \sqrt[3]{\frac{1}{4} - \frac{\sqrt{\frac{43}{3}}}{12}}}}} - \sqrt{2 \sqrt[3]{\frac{1}{4} + \frac{\sqrt{\frac{43}{3}}}{12}} + 2 \sqrt[3]{\frac{1}{4} - \frac{\sqrt{\frac{43}{3}}}{12}}}}{2} \\ \frac{-\sqrt{-2 \sqrt[3]{\frac{1}{4} + \frac{\sqrt{\frac{43}{3}}}{12}} - 2 \sqrt[3]{\frac{1}{4} - \frac{\sqrt{\frac{43}{3}}}{12}} - \frac{4}{\sqrt{2 \sqrt[3]{\frac{1}{4} + \frac{\sqrt{\frac{43}{3}}}{12}} + 2 \sqrt[3]{\frac{1}{4} - \frac{\sqrt{\frac{43}{3}}}{12}}}}} - \sqrt{2 \sqrt[3]{\frac{1}{4} + \frac{\sqrt{\frac{43}{3}}}{12}} + 2 \sqrt[3]{\frac{1}{4} - \frac{\sqrt{\frac{43}{3}}}{12}}}}{2} \\ \frac{\sqrt{-2 \sqrt[3]{\frac{1}{4} + \frac{\sqrt{\frac{43}{3}}}{12}} - 2 \sqrt[3]{\frac{1}{4} - \frac{\sqrt{\frac{43}{3}}}{12}} + \frac{4}{\sqrt{2 \sqrt[3]{\frac{1}{4} + \frac{\sqrt{\frac{43}{3}}}{12}} + 2 \sqrt[3]{\frac{1}{4} - \frac{\sqrt{\frac{43}{3}}}{12}}}}} + \sqrt{2 \sqrt[3]{\frac{1}{4} + \frac{\sqrt{\frac{43}{3}}}{12}} + 2 \sqrt[3]{\frac{1}{4} - \frac{\sqrt{\frac{43}{3}}}{12}}}}{2} \\ \frac{\sqrt{2 \sqrt[3]{\frac{1}{4} + \frac{\sqrt{\frac{43}{3}}}{12}} + 2 \sqrt[3]{\frac{1}{4} - \frac{\sqrt{\frac{43}{3}}}{12}}} - \sqrt{-2 \sqrt[3]{\frac{1}{4} + \frac{\sqrt{\frac{43}{3}}}{12}} - 2 \sqrt[3]{\frac{1}{4} - \frac{\sqrt{\frac{43}{3}}}{12}} + \frac{4}{\sqrt{2 \sqrt[3]{\frac{1}{4} + \frac{\sqrt{\frac{43}{3}}}{12}} + 2 \sqrt[3]{\frac{1}{4} - \frac{\sqrt{\frac{43}{3}}}{12}}}}}}{2} \\ \end{array} \right] \]

Third- and fourth-degree polynomials can be solved in general, with increasingly more complicated answers. Here symbolic_solve needs help: it stops with an error on the general cubic \(a_3x^3 + a_2x^2 + a_1x + a_0\) with all four coefficients symbolic, and also on the monic cubic \(x^3 + bx^2 + cx + d\). (Dividing by the leading coefficient does not change the roots, so the monic cubic is as general.) The classical route, the one Cardano took, is to depress the cubic first: to remove its \(x^2\) term with the substitution \(x = t - b/3\).

@variables b c d t
cubic = x^3 + b*x^2 + c*x + d
depressed = expand(substitute(cubic, x => t - b/3))

\[ t^{3} - \frac{b^{2} t}{3} + c t + \frac{2 b^{3}}{27} - \frac{b c}{3} + d \]

The \(t^2\) terms cancel, leaving a cubic of the form \(t^3 + pt + q\), where \(p = c - b^2/3\) and \(q = 2b^3/27 - bc/3 + d\). This symbolic_solve can solve in general. The following finds one of the answers:

@variables p q
rts = symbolic_solve(t^3 + p*t + q ~ 0, t)
rts[1]   # there are three roots

\[ \sqrt[3]{-\frac{q}{2} + \sqrt{\frac{p^{3}}{27} + \frac{q^{2}}{4}}} + \sqrt[3]{-\frac{q}{2} - \sqrt{\frac{p^{3}}{27} + \frac{q^{2}}{4}}} \]

This is Cardano’s formula. The other two roots are written in the same cube roots, combined with the imaginary unit \(i\). Each root \(t\) of the depressed cubic gives the root \(x = t - b/3\) of the original.

The general formula is also the way to get exact roots of a cubic whose coefficients are numbers. Called directly on one, symbolic_solve writes the pieces of Cardano’s formula as decimals, so its answer is no longer exact. For \(t^3 + t + 1\), the \(31/108\) inside the square root comes out as a long decimal, and its last digits show the rounding:

symbolic_solve(t^3 + t + 1 ~ 0, t)[1]

\[ \sqrt[3]{-0.5 + \sqrt{0.2870370370370370370370370370370370370370370370370370370370370370370370370370362}} + \sqrt[3]{-0.5 - \sqrt{0.2870370370370370370370370370370370370370370370370370370370370370370370370370362}} \]

(This happens for a cubic with number coefficients that does not factor over the rational numbers. A cubic that does factor, and a cubic with symbolic coefficients like the one above, come out exact.) Substituting the numbers into the general formula instead keeps every piece exact. Here \(p = 1\) and \(q = 1\):

substitute(rts[1], Dict(p => 1, q => 1))

\[ \sqrt[3]{-\frac{1}{2} + \sqrt{\frac{31}{108}}} + \sqrt[3]{-\frac{1}{2} - \sqrt{\frac{31}{108}}} \]

Both answers are the same number, about \(-0.6823\); only the second is exact. Sometimes the substituted formula keeps a piece that could be simplified further by hand: for \(x^3 - 2\), where \(p = 0\) and \(q = -2\), it leaves \(\sqrt{1}\) inside the cube roots.

Some fifth degree polynomials are solvable in terms of radicals, however, symbolic_solve will not seem to have luck with this particular fifth degree polynomial:

symbolic_solve(x^5 - x + 1 ~ 0, x)

\[ \left[ \begin{array}{c} \operatorname{roots\_of}\left( x^{5} - x + 1, x \right) \\ \end{array} \right] \]

What comes back is not a root but a placeholder, roots_of, naming the polynomial whose roots it could not express. (Though there is no formula involving only radicals like the quadratic equation, there is a formula for the roots in terms of a function called the Bring radical.)

12.1.3 Repeated roots

A root that is repeated counts once in the answer from symbolic_solve. Here \(1\) and \(2\) are each double roots, but each is listed once:

symbolic_solve((x-1)^2 * (x-2)^2 ~ 0, x)

\[ \left[ \begin{array}{c} 1 \\ 2 \\ \end{array} \right] \]

To see the multiplicities, factor instead. factored_poly shows them as powers:

factored_poly(expand((x-1)^2 * (x-2)^2), x)

\[ \left( x - 1 \right)^{2} \left( x - 2 \right)^{2} \]

and poly_factors lists each factor as often as it is repeated:

poly_factors(expand((x-1)^2 * (x-2)^2), x)

\[ \left[ \begin{array}{c} x - 1 \\ x - 1 \\ x - 2 \\ x - 2 \\ \end{array} \right] \]

For a polynomial with symbolic coefficients, the difference between the symbol and the coefficients must be identified. The second argument of symbolic_solve does that: it names the variable, and every other symbol is treated as a coefficient. Here x is the variable and a, b and c are coefficients:

symbolic_solve(a*x^2 + b*x + c, x)

\[ \left[ \begin{array}{c} \frac{-b + \sqrt{b^{2} - 4 a c}}{2 a} \\ \frac{-b - \sqrt{b^{2} - 4 a c}}{2 a} \\ \end{array} \right] \]

12.1.4 Numerically finding roots

The numeric_roots function gets numeric approximations to all the roots, real and complex:

numeric_roots(x^5 - x + 1, x)

5-element Vector{ComplexF64}:
  -1.1673039782614187 + 0.0im
 -0.18123244446987538 - 1.0839541013177107im
 -0.18123244446987538 + 1.0839541013177107im
   0.7648844336005848 - 0.35247154603172626im
   0.7648844336005848 + 0.35247154603172626im

As roots are complex numbers in general, all five come back as complex numbers, with an imaginary part of 0.0 for the real ones. This polynomial has \(1\) real root, which real_only = true returns on its own, as a real number:

numeric_roots(x^5 - x + 1, x; real_only = true)

1-element Vector{Float64}:
 -1.1673039782614187

Here we see another example:

ex = x^7 -3x^6 +  2x^5 -1x^3 +  2x^2 + 1x^1  - 2
symbolic_solve(ex ~ 0, x)

\[ \left[ \begin{array}{c} 1 \\ 2 \\ \operatorname{roots\_of}\left( x^{5} - x - 1, x \right) \\ \end{array} \right] \]

This finds two of the seven possible roots. The polynomial factors as \((x-1)(x-2)(x^5 - x - 1)\), and the other five roots, those of the fifth-degree factor, are left inside a roots_of. The remainder of the real roots can be found numerically:

numeric_roots(ex, x; real_only = true)

3-element Vector{Float64}:
 1.0
 1.1673039782614187
 2.0

An exact answer can also be correct and still hard to read. The polynomial \(8x^4 - 8x^2 + 1\) has four real roots, and symbolic_solve finds them all, exactly, in nested square roots:

p = 8x^4 - 8x^2  + 1
symbolic_solve(p ~ 0, x)

\[ \left[ \begin{array}{c} \frac{\sqrt{1 - \frac{\sqrt{2}}{2}} - \sqrt{1 + \frac{\sqrt{2}}{2}}}{2} \\ \frac{-\sqrt{1 + \frac{\sqrt{2}}{2}} - \sqrt{1 - \frac{\sqrt{2}}{2}}}{2} \\ \frac{\sqrt{1 + \frac{\sqrt{2}}{2}} + \sqrt{1 - \frac{\sqrt{2}}{2}}}{2} \\ \frac{\sqrt{1 + \frac{\sqrt{2}}{2}} - \sqrt{1 - \frac{\sqrt{2}}{2}}}{2} \\ \end{array} \right] \]

To get numeric approximations, we can use numeric_roots. These roots are all real, so the imaginary parts are all 0.0:

numeric_roots(p, x)

4-element Vector{ComplexF64}:
 -0.9238795325112867 + 0.0im
 -0.3826834323650898 + 0.0im
  0.3826834323650898 + 0.0im
  0.9238795325112867 + 0.0im

(These numbers may look familiar: \(0.9238\dots\) is \(\cos(\pi/8)\) and \(0.3826\dots\) is \(\cos(3\pi/8)\). The polynomial is the Chebyshev polynomial \(T_4\), and the questions at the end of this section look at that family.)

12.2 Do numeric methods matter when you can just graph?

It may seem that certain practices related to roots of polynomials are unnecessary as we could just graph the equation and look for the roots. This feeling is perhaps motivated by the examples given in textbooks to be worked by hand, which necessarily focus on smallish solutions. But, in general, without some sense of where the roots are, an informative graph itself can be hard to produce. That is, technology doesn’t displace thinking—it only supplements it.

For another example, consider the polynomial \((x-20)^5 - (x-20) + 1\). In this form we might think the roots are near \(20\). However, were we presented with this polynomial in expanded form: \(x^5 - 100x^4 + 4000x^3 - 80000x^2 + 799999x - 3199979\), we might be tempted to just graph it to find roots. A naive graph might be to plot over \([-10, 10]\):

p = x^5 - 100x^4 + 4000x^3 - 80000x^2 + 799999x - 3199979
plot(p, -10, 10)

This seems to indicate a root near \(10\). But look at the scale of the \(y\) axis. The value at \(-10\) is around \(-25,000,000\) so it is really hard to tell if \(f\) is near \(0\) when \(x=10\), as the range is too large.

A graph over \([10,20]\) is still unclear:

plot(p, 10,20)

We see that what looked like a zero near \(10\), was actually a number around \(-100,000\).

Continuing, a plot over \([15, 20]\) still isn’t that useful. It isn’t until we get close to \(18\) that the large values of the polynomial allow a clear vision of the values near \(0\). That being said, plotting anything bigger than \(22\) quickly makes the large values hide those near \(0\), and might make us think where the function dips back down there is a second or third zero, when only \(1\) is the case. (We know that, as this is the same \(x^5 - x + 1\) shifted to the right by \(20\) units.)

plot(p, 18, 22)

Not that it can’t be done, but graphically solving for a root here can require some judicious choice of viewing window. Even worse is the case where something might graphically look like a root, but in fact not be a root. Something like \((x-100)^2 + 0.1\) will demonstrate.

For another example, the following polynomial when plotted over \([-5,7]\) appears to have two real roots:

h = x^7 - 16129x^2 + 254x - 1
plot(h, -5, 7)

in fact there are three, two are very close together:

numeric_roots(h, x; real_only = true)

3-element Vector{Float64}:
 0.007874015406930342
 0.007874016089132754
 6.939437409621392

NoteNote

The difference of the two roots is around 7e-10. For the graph over the interval of \([-5,7]\) there are about \(800\) “pixels” used, so each pixel represents a size of about 1.5e-2. So the cluster of roots would safely be hidden under a single “pixel.”

Two floating-point numbers this close together could also be a single double root, split in two by rounding error. root_enclosures settles the question. It works from the exact polynomial, and each real root comes back inside a small interval that is guaranteed to contain it, printed as [midpoint +/- radius]:

enclosures = root_enclosures(h, x)

3-element Vector{Nemo.ArbFieldElem}:
 [0.0078740154069303411575550030281616333765 +/- 5.65e-41]
 [0.0078740160891327544036087278987797271342 +/- 6.32e-41]
 [6.9394374096213921244367134924476102722 +/- 3.83e-38]

The first two intervals do not overlap, so they contain two different roots:

Nemo.overlaps(enclosures[1], enclosures[2])

false

The point of this is to say, that it is useful to know where to look for roots, even if graphing calculators or graphing programs make drawing graphs relatively painless. A better way in this case would be to find the real roots first, and then incorporate that information into the choice of plot window.

12.3 Some facts about the real roots of a polynomial

A polynomial with real coefficients may or may not have real roots. The following discusses some simple checks on the number of real roots and bounds on how big they can be. This can be roughly used to narrow viewing windows when graphing polynomials.

12.3.1 Descartes’ rule of signs

The study of polynomial roots is an old one. In \(1637\) Descartes published a simple method to determine an upper bound on the number of positive real roots of a polynomial.

Descartes’ rule of signs: if \(p=a_n x^n + a_{n-1}x^{n-1} + \cdots + a_1x + a_0\) then the number of positive real roots is either equal to the number of sign differences between consecutive nonzero coefficients, or is less than it by an even number. Repeated roots are counted separately.

One method of proof (sketched at the end of this section) first shows that in synthetic division by \((x-c)\) with \(c > 0\), we must have that any sign change in \(q\) is related to a sign change in \(p\) and there must be at least one more in \(p\). This is then used to show that there can be only as many positive roots as sign changes. That the difference comes in pairs is related to complex roots of real polynomials always coming in pairs.

An immediate consequence, is that a polynomial whose coefficients are all non-negative will have no positive real roots.

Applying this to the polynomial \(x^5 -x + 1\) we get that the coefficients have signs: + 0 0 0 - + which collapses to the sign pattern +, -, +. This pattern has two changes of sign. The number of positive real roots is either \(2\) or \(0\). In fact there are \(0\) for this case.

What about negative roots? Clearly, any negative root of \(p\) is a positive root of \(q(x) = p(-x)\), as the graph of \(q\) is just that of \(p\) flipped through the \(y\) axis. But the coefficients of \(q\) are the same as \(p\), except for the odd-indexed coefficients (\(a_1, a_3, \dots\)) have a changed sign. Continuing with our example, for \(q(x) = -x^5 + x + 1\) we get the new sign pattern -, +, + which yields one sign change. That is, there must be a negative real root, and indeed there is, \(x \approx -1.1673\).

With this knowledge, we could have known that in an earlier example the graph of p = x^7 - 16129x^2 + 254x - 1 – which indicated two positive real roots – was misleading, as there must be \(1\) or \(3\) by a count of the sign changes.

For another example, if we looked at \(f(x) = x^5 - 100x^4 + 4000x^3 - 80000x^2 + 799999x - 3199979\) again, we see that there could be \(1\), \(3\), or \(5\) positive roots. However, changing the signs of the odd powers leaves all “-” signs, so there are \(0\) negative roots. From the graph, we saw just \(1\) real root, not \(3\) or \(5\). We can verify numerically with:

j =  x^5 - 100x^4 + 4000x^3 - 80000x^2 + 799999x - 3199979
numeric_roots(j, x)

5-element Vector{ComplexF64}:
  18.83269602173858 + 0.0im
 19.818767555530126 - 1.0839541013177107im
 19.818767555530126 + 1.0839541013177107im
 20.764884433600585 - 0.35247154603172626im
 20.764884433600585 + 0.35247154603172626im

12.3.2 Cauchy’s bound on the magnitude of the real roots.

Descartes’ rule gives a bound on how many real roots there may be. Cauchy provided a bound on how large they can be. Assume our polynomial is monic (if not, divide by \(a_n\) to make it so, as this won’t affect the roots). Then any real root is no larger in absolute value than \(h = 1 + |a_0| + |a_1| + |a_2| + \cdots + |a_{n-1}|\), (this is expressed in different ways.)

To see precisely why this bound works, suppose \(x\) is a root with \(|x| > 1\) and let \(h\) be the bound. Then since \(x\) is a root, we can solve \(a_0 + a_1x + \cdots + 1 \cdot x^n = 0\) for \(x^n\) as:

\[ x^n = -(a_0 + a_1 x + \cdots a_{n-1}x^{n-1}) \]

Which after taking absolute values of both sides, yields by the triangle inequality:

\[ |x^n| \leq |a_0| + |a_1||x| + |a_2||x^2| + \cdots |a_{n-1}| |x^{n-1}| \leq (h-1) (1 + |x| + |x^2| + \cdots |x^{n-1}|). \]

The last sum can be computed using a formula for geometric sums, \((|x^n| - 1)/(|x|-1)\). Rearranging, gives the inequality:

\[ |x| - 1 \leq (h-1) \cdot (1 - \frac{1}{|x^n|} ) \leq (h-1) \]

from which it follows that \(|x| \leq h\), as desired.

For our polynomial \(x^5 -x + 1\) we have the sum above is \(3\). The lone real root is approximately \(-1.1673\) which satisfies \(|-1.1673| \leq 3\).

12.4 Questions

Question

What is the remainder of dividing \(x^4 - x^3 - x^2 + 2\) by \(x-2\)?

Select an item

\(6\)

\(x^3 + x^2 + x + 2\)

\(x-2\)

\(0\)

Question

What is the remainder of dividing \(x^4 - x^3 - x^2 + 2\) by \(x^3 - 2x\)?

Select an item

\(x^2 - 2x + 2\)

\(x - 1\)

\(2\)

Question

We have that \(x^5 - x + 1 = (x^3 + x^2 - 1) \cdot (x^2 - x + 1) + (-2x + 2)\).

What is the remainder of dividing \(x^5 - x + 1\) by \(x^2 - x + 1\)?

Select an item

\(x^3 + x^2 - 1\)

\(-2x + 2\)

\(x^2 - x + 1\)

Question

Consider this output from synthetic division

2 | 1 0 0 0 -1  1
  |   2 4 8 16 30
  ---------------
    1 2 4 8 15 31

representing \(p(x) = q(x)\cdot(x-c) + r\).

What is \(p(x)\)?

Select an item

\(x^4 +2x^3 + 4x^2 + 8x + 15\)

\(2x^4 + 4x^3 + 8x^2 + 16x + 30\)

\(x^5 + 2x^4 + 4x^3 + 8x^2 + 15x + 31\)

\(31\)

\(x^5 - x + 1\)

What is \(q(x)\)?

Select an item

\(x^4 +2x^3 + 4x^2 + 8x + 15\)

\(31\)

\(x^5 + 2x^4 + 4x^3 + 8x^2 + 15x + 31\)

\(x^5 - x + 1\)

\(2x^4 + 4x^3 + 8x^2 + 16x + 30\)

What is \(r\)?

Select an item

\(x^5 + 2x^4 + 4x^3 + 8x^2 + 15x + 31\)

\(x^5 - x + 1\)

\(2x^4 + 4x^3 + 8x^2 + 16x + 30\)

\(31\)

\(x^4 +2x^3 + 4x^2 + 8x + 15\)

Question

Let \(p=x^4 -9x^3 +30x^2 -44x + 24\)

Factor \(p\). What are the factors?

Select an item

\((x+2)\) and \((x+3)\)

\((x-2)\) and \((x-3)\)

\(2\) and \(3\)

Question

Does the expression \(x^4 - 5\) factor over the rational numbers?

Select an item

Yes

No

Using numeric_roots, how many real roots does \(x^4 - 5\) have? (Then compare what symbolic_solve returns: every root it gives is written with an \(i\). An exact form can look complex and still be real.)

Question

The Soviet historian I. Y. Depman claimed that in \(1486\), Spanish mathematician Valmes was burned at the stake for claiming to have solved the quartic equation. Here we don’t face such consequences.

Find the largest real root of \(x^4 - 10x^3 + 32x^2 - 38x + 15\).

Question

What are the numeric values of the real roots of \(f(x) = x^6 - 5x^5 + x^4 - 3x^3 + x^2 - x + 1\)?

Select an item

[-0.434235+0.613836im, -0.434235-0.613836im]

[0.578696, 4.91368]

[-0.434235, -0.434235, 0.188049, 0.188049]

[-0.434235, -0.434235, 0.188049, 0.188049, 0.578696, 4.91368]

Question

Odd polynomials must have at least one real root.

Consider the polynomial \(x^5 - 3x + 1\). Does it have more than one real root?

Select an item

Yes

No

Consider the polynomial \(x^5 - 1.5x + 1\). Does it have more than one real root?

Select an item

Yes

No

Question

What is the maximum number of positive, real roots that Descartes’ bound says \(p=x^5 + x^4 - x^3 + x^2 + x + 1\) can have?

How many positive, real roots does it actually have?

What is the maximum number of negative, real roots that Descartes’ bound says \(p=x^5 + x^4 - x^3 + x^2 + x + 1\) can have?

How many negative, real roots does it actually have?

Question

Let \(f(x) = x^5 - 4x^4 + x^3 - 2x^2 + x\). What does Cauchy’s bound say is the largest possible magnitude of a root?

What is the largest magnitude of a real root?

Question

As \(1 + 2 + 3 + 4\) is \(10\), Cauchy’s bound says that the magnitude of the largest real root of \(x^3 - ax^2 + bx - c\) is at most \(10\) where \(a,b,c\) is one of \(2,3,4\). By considering all 6 such possible polynomials (such as \(x^3 - 3x^2 + 2x - 4\)) what is the largest magnitude of a root?

Question

The roots of the Chebyshev polynomials are helpful for some numeric algorithms. These are a family of polynomials related by \(T_{n+1}(x) = 2xT_n(x) - T_{n-1}(x)\) (a recurrence relation in the manner of the Fibonacci sequence). The first two are \(T_0(x) = 1\) and \(T_1(x) =x\).

• Based on the relation, figure out \(T_2(x)\). It is

Select an item

\(x\)

\(2x\)

\(2x^2\)

\(2x^2 - 1\)

• True or false, the \(degree\) of \(T_n(x)\) is \(n\): (Look at the defining relation and reason this out).

Select an item

Yes

No

• The fifth one is \(T_5(x) = 16x^5 - 20x^3 + 5x\). Cauchy’s bound says that the largest root has absolute value

1 + 20/16 + 5/16

2.5625

The Chebyshev polynomials have the property that in fact all \(n\) roots are real, distinct, and in \([-1, 1]\). Using numeric_roots, find the magnitude of the largest root:

• Plotting p over the interval \([-2,2]\) does not help graphically identify the roots:

plot(16x^5 - 20x^3 + 5x, -2, 2)

Does graphing over \([-1,1]\) show clearly the \(5\) roots?

Select an item

Yes

No

12.5 Appendix: Proof of Descartes’ rule of signs

Proof modified from this post.

First, we can assume \(p\) is monic (\(p_n=1\) and positive), and \(p_0\) is non zero. The latter, as we can easily deflate the polynomial by dividing by \(x\) if \(p_0\) is zero.

Let var(p) be the number of sign changes and pos(p) the number of positive real roots of p.

First: For a monic \(p\) if \(p_0 < 0\) then var(p) is odd and if \(p_0 > 0\) then var(p) is even.

This is true for degree \(n=1\) the two sign patterns under the assumption are +- (\(p_0 < 0\)) or ++ (\(p_0 > 0\)). If it is true for degree \(n-1\), then the we can consider the sign pattern of such an \(n\) degree polynomial having one of these patterns: +...+- or +...-- (if \(p_0 < 0\)) or +...++ or +...-+ if (\(p_0>0\)). An induction step applied to all but the last sign for these four patterns leads to even, odd, even, odd as the number of sign changes. Incorporating the last sign leads to odd, odd, even, even as the number of sign changes.

Second: For a monic \(p\) if \(p_0 < 0\) then pos(p) is odd, if \(p_0 > 0\) then pos(p) is even.

This is clearly true for monic degree \(1\) polynomials: if \(c\) is positive \(p = x - c\) has one real root (an odd number) and \(p = x + c\) has \(0\) real roots (an even number). Now, suppose \(p\) has degree \(n\) and is monic. Then as \(x\) goes to \(\infty\), it must be \(p\) goes to \(\infty\).

If \(p_0 < 0\) then there must be a positive real root, say \(r\), (Bolzano’s intermediate value theorem). Dividing \(p\) by \((x-r)\) to produce \(q\) requires \(q_0\) to be positive and of lower degree. By induction \(q\) will have an even number of roots. Add in the root \(r\) to see that \(p\) will have an odd number of roots.

Now consider the case \(p_0 > 0\). There are two possibilities either pos(p) is zero or positive. If pos(p) is \(0\) then there are an even number of roots. If pos(p) is positive, then call \(r\) one of the real positive roots. Again divide by \(x-r\) to produce \(p = (x-r) \cdot q\). Then \(q_0\) must be negative for \(p_0\) to be positive. By induction \(q\) must have an odd number of roots, meaning \(p\) must have an even numbers.

So there is parity between var(p) and pos(p): if \(p\) is monic and \(p_0 < 0\) then both var(p) and pos(p) are both odd; and if \(p_0 > 0\) both var(p) and pos(p) are both even.

Descartes’ rule of signs will be established if it can be shown that var(p) is at least as big as pos(p). Suppose \(r\) is a positive real root of \(p\) with \(p = (x-r)q\). We show that var(p) > var(q) which can be repeatedly applied to show that if \(p=(x-r_1)\cdot(x-r_2)\cdot \cdots \cdot (x-r_l) q\), where the \(r_i\)s are the positive real roots, then var(p) >= l + var(q) >= l = pos(p).

As \(p = (x-c)q\) we must have the leading term is \(p_nx^n = x \cdot q_{n-1} x^{n-1}\) so \(q_{n-1}\) will also be + under our monic assumption. Looking at a possible pattern for the signs of \(q\), we might see the following unfinished synthetic division table for a specific \(q\):

  + ? ? ? ? ? ? ? ?
+   ? ? ? ? ? ? ? ?
  -----------------
  + - - - + - + + 0

But actually, we can fill in more, as the second row is formed by multiplying a positive \(c\):

  + ? ? ? ? ? ? ? ?
+   + - - - + - + +
  -----------------
  + - - - + - + + 0

What’s more, using the fact that to get 0 the two summands must differ in sign and to have a ? plus + yield a -, the ? must be - (and reverse), the following must be the case for the signs of p:

  + - ? ? + - + ? -
+   + - - - + - + +
  -----------------
  + - - - + - + + 0

If the bottom row represents \(q_7, q_6, \dots, q_0\) and the top row \(p_8, p_7, \dots, p_0\), then the sign changes in \(q\) from + to - are matched by sign changes in \(p\). The ones in \(q\) from \(-\) to \(+\) are also matched regardless of the sign of the first two question marks (though \(p\) could possibly have more). The last sign change in \(p\) between \(p_2\) and \(p_0\) has no counterpart in \(q\), so there is at least one more sign change in \(p\) than \(q\).

As such, the var(p) \(\geq 1 +\) var(q).