# Real solutions of a zero-dimensional system
#
# Maple version of the worksheet 03_zero_dimensional_systems.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 03_zero_dimensional_systems.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 system of polynomial equations with rational coefficients that has finitely
# many complex solutions (a *zero-dimensional* system). The output is the list of its real
# solutions, each one given by a box with rational bounds holding exactly one solution, all
# the real solutions being found.
#
# PACE goes through a *Rational Univariate Representation* (RUR): a univariate polynomial
# f(T) and rational functions that give the coordinates of the solutions at the roots of f.
# Solving the system then comes down to the real roots of f (see the worksheet "Real roots of
# a univariate polynomial"), plus the evaluation of the coordinates.

interface(prettyprint = 0, warnlevel = 0):
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 box (a list of INTERVAL) with 10-digit bounds
fbox := b -> evalf[10](bounds(b)):

# ## 1. A first system
#
# The intersection of the cubic y^2 = x^3 - x with the circle x^2 + y^2 = 3:

sys := [y^2 - x^3 + x, x^2 + y^2 - 3]:

# real_solutions gives the real solutions, as boxes of intervals INTERVAL(lo .. hi) (one per
# variable, in the order of the variables: here x, y, sorted by name) of 64 significant bits:

sols := real_solutions(sys, prec = 64):
printf("## 1. %a\n%d real solutions:\n", sys, nops(sols));
for b in sols do printf("  %a\n", fbox(b)) end do:

# Two real solutions, symmetric with respect to the x-axis.

# ## 2. The Rational Univariate Representation
#
# zdim_parameterization returns the RUR itself: the eliminant f, then one numerator g_k per
# variable, all given by their coefficients from the constant term; the solutions are
# x_k = g_k(T)/f'(T) at the roots T of f. The *separating form* sep is the linear form in the
# variables that takes the value T at each solution (here its coefficients on x, y; a name
# sep would be evaluated in the keyword sep = sep, hence sf):

rs := zdim_parameterization(sys, get_separating_element = true, verbose = false):
rur, sf := rs[1], rs[2]:
f := add(rur[1][i] * T^(i - 1), i = 1 .. nops(rur[1])):
printf("\n## 2. The RUR\nseparating form: %a\neliminant: %a\n", sf, f);

# Here T = y. The degree of f, 6, is the number of complex solutions of the system, counted
# with multiplicities: the RUR describes all of them, and only the real roots of f give real
# solutions. Its real roots are +-sqrt(a) with a the real root of the cubic
# a^3 - 6a^2 + 16a - 12, two of them:

printf("real roots of the eliminant:\n");
show(isolate(isolation(rur[1])));

# The coordinates being rational functions of T, a solution is known to any precision by
# refining the root of f only, without solving the system again. To 200 bits, from the same
# RUR:

sols200 := real_solutions(rur, prec = 200, sep = sf):
printf("first solution to 200 bits:\n  %a\n", evalf[62](bounds(sols200[1])));

# ## 3. A check
#
# Evaluated in interval arithmetic on each box, each equation gives an interval holding 0:

printf("\n## 3. Interval evaluation of the equations on the boxes\n");
for b in sols do printf("  %a\n", map(has_zero, evaluate(sys, b, precision = 128))) end do:

# This is only a necessary condition, not a proof: the proof comes from the RUR, whose roots
# are isolated exactly.

# ## 4. Multiple solutions
#
# The line y = 1 is tangent to the parabola y = x^2 - 2x + 2 at (1, 1): a solution of
# multiplicity 2. With multiplicities = true the RUR comes with, for each multiplicity j, the
# polynomial whose roots are the values of T at the solutions of multiplicity j:

tangent := [y - x^2 + 2*x - 2, y - 1]:
rt := zdim_parameterization(tangent, multiplicities = true, verbose = false):
topoly := c -> add(c[i] * T^(i - 1), i = 1 .. nops(c)):
printf("\n## 4. %a\nmultiplicity = polynomial: %a\n", tangent, map(e -> lhs(e) = topoly(rhs(e)), rt[2]));

# One solution of multiplicity 2, and it is real:

printf("real solutions: %a\n", map(fbox, real_solutions(tangent, prec = 64)));

# ## 5. Two equations in two unknowns
#
# For two equations in two unknowns, biv_real_solutions takes another route: it eliminates a
# variable by subresultants, without Groebner basis. The solutions are the same:

printf("\n## 5. biv_real_solutions\n");
for b in biv_real_solutions(op(sys), prec = 64) do printf("  %a\n", fbox(b)) end do:

# On two random curves of degree 16 (256 complex solutions), drawn by Maple:

dense := proc(d) local r;
    r := rand(-100 .. 100);
    add(add(r() * x^i * y^j, j = 0 .. d - i), i = 0 .. d)
end proc:
randomize(3):
a, b := dense(16), dense(16):
t := time[real](): s_biv := biv_real_solutions(a, b, prec = 64): t_biv := time[real]() - t:
t := time[real](): s_rur := real_solutions([a, b], prec = 64): t_rur := time[real]() - t:
printf("biv_real_solutions: %d real solutions in %.2f s\n", nops(s_biv), t_biv);
printf("real_solutions:     %d real solutions in %.2f s\n", nops(s_rur), t_rur);

# Both find the same real solutions among the 256 complex ones, in a fraction of a second
# (the times include the transfer between Maple and Julia).

# ## 6. More variables
#
# The Katsura systems, n + 1 equations in n + 1 unknowns u0, ..., un with 2^n complex
# solutions, are a classical benchmark:

katsura := proc(n)
    local U, eqs, m, l;
    U := i -> `if`(abs(i) <= n, u || (abs(i)), 0);
    eqs := [seq(add(U(l) * U(m - l), l = -n .. n) - U(m), m = 0 .. n - 1)];
    [op(eqs), add(U(l), l = -n .. n) - 1]
end proc:
printf("\n## 6. Katsura systems\nkatsura(2) = %a\n", katsura(2));

# The number of real solutions among the 2^n complex ones, and the time:

for n in [4, 6, 8] do
    t := time[real]();
    s := real_solutions(katsura(n), prec = 64);
    printf("  n = %d: %3d complex, %3d real, %.2f s\n", n, 2^n, nops(s), time[real]() - t);
end do: