# Real roots of a univariate polynomial
#
# Maple version of the worksheet 01_univariate_real_roots.jl: the functions of PACE.jl are
# called from Maple through UDX (see "Using PACE.jl in Maple with UDX.jl"). Run in batch with
#
#     maple -q 01_univariate_real_roots.mpl
#
# with UDX_MPL_LIB set, and UDX_SRC, PACE_SRC, JULIA, JULIA_PROJECT, UDX_LIB set or replaced
# below by their paths.
#
# The input is a polynomial f in 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.

interface(prettyprint = 0, warnlevel = 0, rtablesize = 20):
read cat(getenv("UDX_SRC"), "/Maple/udx.mpl"):
read cat(getenv("PACE_SRC"), "/maple/PACE.mpl"):
UDX:-udx_set_julia(getenv("JULIA"), env = ["JULIA_PROJECT" = getenv("JULIA_PROJECT"),
                                           "UDX_LIB" = getenv("UDX_LIB")]):
with(PACE):

# a root of PACE (a JuliaObject): its bounds, as rationals or 10-digit floats
bnds := r -> [lo(r), hi(r)]:
fbnds := r -> evalf[10](bnds(r)):
# is the exact real number c inside the interval of r? (decided with d digits)
inside := (r, c, d) -> evalb(evalf[d](lo(r) - c) < 0 and evalf[d](hi(r) - c) > 0):

# ## 1. The input
#
# This polynomial has a double root 1/2, a triple pair +-sqrt(2) and the real root of
# x^3 - x - 1 (the plastic number), the other two roots of this cubic being complex. A
# univariate polynomial is sent to Julia as the vector of its coefficients, from the constant
# term to the leading one.

f := expand((x - 1/2)^2 * (x^2 - 2)^3 * (x^3 - x - 1)):
printf("## 1. The input\nf = %a\n", factor(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. The result comes back as coefficients.

sf := squarefree_part(f):
printf("\n## 2. Square-free part\ncoefficients: %a\npolynomial:   %a\n", sf,
       factor(add(sf[i] * x^(i - 1), i = 1 .. nops(sf))));

# ## 3. Isolation
#
# isolation isolates the real roots and keeps the result in Julia (a JuliaObject), so that
# it can be refined later without starting over; isolate gives its roots, in increasing
# order. Each root is 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.

iso := isolation(f):
rts := isolate(iso):
printf("\n## 3. Isolation\n%d real roots:\n", nops(rts));
show(rts);

# Four real roots, as expected; the intervals are as wide as possible: the isolation stops
# as soon as the roots are separated.

# ## 4. Refinement
#
# refine shrinks the intervals of the whole isolation to a width below 2^(-prec):

refine(iso, 64):
rts := isolate(iso):
printf("\n## 4. Refinement\nto 64 bits:\n");
for r in rts do printf("  %a\n", fbnds(r)) end do:

# 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. The bounds are exact rationals, compared with
# sqrt(2) evaluated by Maple with 400 digits.

r := rts[-1]:
refine(r, 1000):
printf("width of the interval of sqrt(2): %.3e\n", evalf(width(r)));
printf("lo < sqrt(2) < hi: %a\n", inside(r, sqrt(2), 400));

# ## 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. Its arguments: the square-free part, the polynomials, the intervals of
# the roots.

res := signs_at(sf, [x^2 - 2, 2*x - 1], map(bnds, rts)):
printf("\n## 5. Signs of x^2 - 2 and 2x - 1 at the roots\n");
for i to nops(rts) do printf("  root %d: %a\n", i, res[2][i]) end do:

# 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 (100 x - 1)^2 has two real roots very close to 1/100: they
# differ by about 10^(-22). The usual floating-point method computes the eigenvalues of the
# companion matrix (this is what roots does in Matlab), here in hardware Float64 arithmetic.

m := expand(x^20 - 2 * (100*x - 1)^2):
n := degree(m, x):
C := Matrix(n, n, datatype = float[8]):
for i from 2 to n do C[i, i - 1] := 1 end do:
for i to n do C[i, n] := -coeff(m, x, i - 1) / lcoeff(m, x) end do:
ev := LinearAlgebra:-Eigenvalues(C):
printf("\n## 6. Mignotte's polynomial\neigenvalues of the companion matrix near the real axis:\n");
for z in ev do if abs(Im(z)) < 1e-6 then printf("  %a\n", z) end if end do:

# 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(m):
printf("certified isolation: %d real roots\n", nops(isolate(isom)));

# and separates the two close ones once refined:

refine(isom, 100):
rm := isolate(isom):
printf("gap between the 2nd and 3rd intervals: %.3e\n", evalf(lo(rm[3]) - 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 +-1. With n = 300:

ct := expand(orthopoly[T](300, x)):
t := time[real]():
isot := isolation(ct):
nb := op(2, select(e -> lhs(e) = nb, fields(isot))[1]):
printf("\n## 7. Chebyshev T_300\n%d real roots, isolated in %.2f s (with the transfer)\n",
       nb, time[real]() - t);

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

t := time[real]():
refine(isot, 256):
printf("refined to 256 bits in %.2f s\n", time[real]() - t);

# The extreme root is cos(pi/600):

r1 := isolate(isot)[-1]:
printf("lo < cos(pi/600) < hi: %a\n", inside(r1, cos(Pi/600), 120));