Real roots of a univariate polynomial

Download: Jupyter notebook, Julia script, Maple script (PACE called from Maple through UDX, see Using PACE.jl in Maple). To install PACE and run them, see Running the worksheets.

The input is a polynomial $f \in \mathbb{Q}[x]$, nothing else. The output is a list of disjoint intervals with rational bounds, each holding exactly one real root of $f$, all the real roots being found. The result is exact: the intervals can then be shrunk to any width, the signs of other polynomials at the roots can be decided, and no rounding error can make a root disappear or appear.

We go through the steps one by one: square-free part, isolation, refinement, then two situations where floating-point arithmetic gives a wrong answer and certified isolation does not.

using PACE, AbstractAlgebra, LinearAlgebra

1. The input

A polynomial is built with AbstractAlgebra (or Nemo). This one has a double root $1/2$, a triple pair $\pm\sqrt 2$ and the real root of $x^3 - x - 1$ (the plastic number), the other two roots of this cubic being complex:

R, x = polynomial_ring(QQ, "x")
f = (x - 1//2)^2 * (x^2 - 2)^3 * (x^3 - x - 1)

The solvers of PACE read a polynomial as the vector of its coefficients, from the constant term to the leading one:

c = Rational{BigInt}.(collect(coefficients(f)))

12-element Vector{Rational{BigInt}}:
    2
   -6
   -3
   15
  19//2
 -43//2
 -49//4
  69//4
    6
 -27//4
   -1
    1

2. Square-free part

A multiple root is a root of $f$ and of $f'$. Dividing $f$ by $\gcd(f, f')$ keeps each root once, which is what the isolation works on:

squarefree_part(f)

The isolation below takes it itself; it is shown here only to see what is isolated.

3. Isolation

isolation isolates the real roots and keeps the result, so that it can be refined later without starting over:

iso = isolation(c)

Isolation(4 roots, degree 6, not refined)

isolate gives its roots, in increasing order:

roots = isolate(iso)

4-element Vector{LACE.Root}:
 Root(-1.5 .. -1.0)
 Root(0.0 .. 1.0)
 Root(1.25 .. 1.375)
 Root(1.375 .. 1.5)

Each root is a Root: either an open interval $]l, r[$ with rational bounds holding exactly one root, or a single rational, equal to the root, when the isolation met the root exactly. Four real roots, as expected; the intervals are as wide as possible: the isolation stops as soon as the roots are separated.

4. Refinement

PACE.refine! shrinks the intervals of the whole isolation to a width below $2^{-\text{prec}}$:

PACE.refine!(iso, 64)
roots = isolate(iso)

4-element Vector{LACE.Root}:
 Root(-1.4142135623730951 .. -1.4142135623730951)
 Root(1//2)
 Root(1.324717957244746 .. 1.324717957244746)
 Root(1.4142135623730951 .. 1.4142135623730951)

The double root $1/2$ is now found exactly. A single root can also be refined alone, to any precision: here $\sqrt 2$ to 1000 bits.

r = roots[end]
PACE.refine!(r, 1000)
Float64(PACE.width(r))

5.832897615645118e-303

The bounds are exact rationals; compared with $\sqrt 2$ computed by MPFR:

s = setprecision(() -> sqrt(BigFloat(2)), BigFloat, 1100)
PACE.lo(r) < s < PACE.hi(r)

true

5. Signs of other polynomials at the roots

Which roots of $f$ are roots of $x^2 - 2$, which ones lie above $1/2$? signs_at gives the sign of each polynomial at each root, a zero being decided exactly (by a gcd), not by an approximate value. The polynomials are given by their integer coefficients (PACE.zx_integral clears the denominators):

zx(p) = PACE.zx_integral(Rational{BigInt}.(collect(coefficients(p))))
boxes, signs = signs_at(iso.sf, [zx(x^2 - 2), zx(2x - 1)], [(PACE.lo(r), PACE.hi(r)) for r in roots])
signs

4-element Vector{Vector{Int64}}:
 [0, -1]
 [-1, 0]
 [-1, 1]
 [0, 1]

Row $i$ gives the signs of $x^2-2$ and $2x-1$ at the $i$-th root: $x^2 - 2$ vanishes at the first and last roots, $2x - 1$ at the second one.

6. Two close roots: where floating point fails

Mignotte's polynomial $x^{20} - 2(100x - 1)^2$ has two real roots very close to $1/100$: they differ by about $10^{-22}$.

m = x^20 - 2 * (100x - 1)^2
cm = Rational{BigInt}.(collect(coefficients(m)));

The usual floating-point method computes the eigenvalues of the companion matrix (this is what roots does in Matlab):

n = degree(m)
C = zeros(n, n)
for i in 2:n
    C[i, i-1] = 1
end
C[:, n] = -Float64.(cm[1:n]) ./ Float64(cm[n+1])
filter(z -> abs(imag(z)) < 1e-6, eigvals(C))

4-element Vector{ComplexF64}:
   -1.734696440260734 + 0.0im
 0.009999999855836147 + 0.0im
  0.01000000014416385 + 0.0im
   1.7324741845653981 + 0.0im

Near $1/100$ it finds a pair of complex conjugate roots, with an imaginary part of about $10^{-10}$: it cannot tell two close real roots from two complex ones, and here it picks the wrong answer. The certified isolation finds four real roots:

isom = isolation(cm)
isolate(isom)

4-element Vector{LACE.Root}:
 Root(-4.0 .. 0.0)
 Root(0.01 .. 0.01)
 Root(0.01 .. 0.01)
 Root(1.0 .. 2.0)

and separates the two close ones once refined:

PACE.refine!(isom, 100)
rm = isolate(isom)
Float64(PACE.lo(rm[3]) - PACE.hi(rm[2]))

1.4142135619248986e-22

The gap between the two intervals is a proof that the roots are distinct.

7. Higher degrees

The Chebyshev polynomial $T_n$ has $n$ real roots in $]-1, 1[$, crowded near $\pm 1$. With $n = 300$:

function chebyshev(n)
    a, b = one(R), x
    for _ in 2:n
        a, b = b, 2x * b - a
    end
    n == 0 ? a : b
end
ct = Rational{BigInt}.(collect(coefficients(chebyshev(300))))
t = @elapsed isot = isolation(ct)
length(isot), t

(300, 0.147261993)

All 300 roots are found, in a fraction of a second. Their refinement to 256 bits:

@elapsed PACE.refine!(isot, 256)

0.362037481

The extreme root is $\cos(\pi/600)$:

r1 = isolate(isot)[end]
c1 = setprecision(() -> cos(BigFloat(pi) / 600), BigFloat, 300)
PACE.lo(r1) < c1 < PACE.hi(r1)

true