# # Real roots of a univariate polynomial
#
# 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.

#nb # Notebook only: Nemo and AbstractAlgebra objects are shown as plain text
#nb # (their HTML and LaTeX output throws on purpose, which IJulia 1.34 no longer catches).
#nb filter!(m -> !(m == MIME("text/latex") || (m isa AbstractVector && MIME("text/html") in m)), IJulia.ijulia_mime_types)

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

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

# `isolate` gives its roots, in increasing order:

roots = isolate(iso)

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

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

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

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

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

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

# and separates the two close ones once refined:

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

# 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

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

@elapsed PACE.refine!(isot, 256)

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