# Complex solutions of a zero-dimensional system
#
# Maple version of the worksheet 07_complex_solutions.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 07_complex_solutions.mpl
#
# with UDX_MPL_LIB set, and UDX_SRC, PACE_SRC, JULIA, JULIA_PROJECT, UDX_LIB set or replaced
# below by their paths.
#
# Experimental: PWPoly is an experimental add-on of PACE, not installed with it. It must be
# added once to the Julia environment of UDX (import Pkg; Pkg.add("PWPoly")); load("PWPoly")
# below loads it. It is not available on Windows: this worksheet runs on Linux and macOS.
#
# The worksheet "Real solutions of a zero-dimensional system" computed the real solutions
# from a Rational Univariate Representation (RUR): an eliminant f(T) and coordinates
# x_k = g_k(T)/f'(T) at its roots. The same RUR describes *all* the complex solutions: we
# isolate the complex roots of f in certified balls with PWPoly, then evaluate the coordinates
# on these balls in ball arithmetic. Each solution is given by a box of the complex space,
# with a proof that it holds it.

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):
load("PWPoly"):

# The complex balls (of Arb, in Nemo) stay in Julia, as JuliaObject; the operations on them
# are those of Julia, called by callpkg. Their real and imaginary parts come to Maple as
# intervals INTERVAL(lo .. hi) holding them.
cdiv := (a, b) -> callpkg("Base", "/", a, b):
re := z -> callpkg("Nemo", "real", z):
im := z -> callpkg("Nemo", "imag", z):
rad := J -> (rhs(op(1, J)) - lhs(op(1, J))) / 2:
# does the ball p(s) hold 0? (p a polynomial in the variables vars, s a list of balls)
holds_zero := (p, s, vars) -> callpkg("Nemo", "contains_zero",
                                      callpkg("Nemo", "evaluate", p, s, ring = vars)):
# the precision of the balls of Nemo (a global of Julia): Balls is the Julia type
Balls := julia_global("Nemo", "Balls"):

# ## 1. The complex roots of the eliminant
#
# The system of the worksheet on real solutions, the cubic y^2 = x^3 - x and the circle
# x^2 + y^2 = 3:

sys := [y^2 - x^3 + x, x^2 + y^2 - 3]:
rur := zdim_parameterization(sys, verbose = false):
f := rur[1]:
printf("## 1. %a\neliminant: %a\n", sys, add(f[i] * T^(i - 1), i = 1 .. nops(f)));

# The eliminant T^6 - 6T^4 + 16T^2 - 12 has degree 6: the system has 6 complex solutions.
# isolate_complex_roots_sqrfree_poly isolates the complex roots of a square-free polynomial
# with integer coefficients (the denominators are cleared for it). With acb = true the roots are balls of
# Arb, each holding exactly one root; log2eps sets their relative precision. It returns the
# roots, which ones are real, and their number:

res := isolate_complex_roots_sqrfree_poly(f, log2eps = -64, acb = true):
Tr, is_real, nb_real := op(res):
show(Tr);
printf("real: %a, %d real roots\n", is_real, nb_real);

# Two real roots, +-sqrt(a) with a ~ 1.15, and two pairs of complex conjugate roots.

# ## 2. The coordinates
#
# The coordinates are x_k = g_k(T)/f'(T). evaluate_at_complex_points evaluates a polynomial
# with integer coefficients on balls; the numerators and f' are multiplied by the same
# integer L, which does not change their quotient:

solutions_from_rur := proc(rur, Tr, bits)
    local f, gs, df, L, den, nums, k, i;
    f, gs := rur[1], rur[2 .. -1];
    df := [seq(k * f[k + 1], k = 1 .. nops(f) - 1)];
    L := ilcm(seq(seq(denom(q), q in g), g in [df, op(gs)]));
    den := evaluate_at_complex_points(L * df, Tr, precision = bits);
    nums := [seq(evaluate_at_complex_points(L * g, Tr, precision = bits), g in gs)];
    [seq([seq(cdiv(nums[k][i], den[i]), k = 1 .. nops(gs))], i = 1 .. nops(Tr))]
end proc:

callpkg("Nemo", "set_precision!", Balls, 64):
sols := solutions_from_rur(rur, Tr, 64):
printf("\n## 2. The six solutions (x, y)\n");
show(sols);

# Six solutions (x, y): the two real ones, and two pairs of conjugate ones.

# ## 3. Checks
#
# Each equation, evaluated on each box in ball arithmetic, gives a ball holding 0:

printf("\n## 3. Checks\nthe equations hold 0 on the boxes: %a\n",
       [seq([seq(holds_zero(p, s, [x, y]), p in sys)], s in sols)]);

# The real solutions are those of real_solutions: their boxes meet (the real parts of the
# balls, as intervals, compared in Maple).

meets := (J, K) -> evalb(lhs(op(1, J)) <= rhs(op(1, K)) and lhs(op(1, K)) <= rhs(op(1, J))):
reals := real_solutions(sys, prec = 64):
printf("the real ones meet the boxes of real_solutions: %a\n",
       [seq(`if`(is_real[i], ormap(r -> andmap(k -> meets(re(sols[i][k]), r[k]), [1, 2]), reals),
                 NULL), i = 1 .. nops(sols))]);

# ## 4. Precision
#
# The coefficients of a RUR are large, and their evaluation loses about as many bits as they
# have. To get the solutions with b bits, the roots of f and the evaluation need more: the
# precision is doubled until the boxes are narrow enough. The balls of Nemo work at the
# precision set by set_precision!(Balls, ...).

bwidth := sols -> max(seq(seq(max(rad(re(c)), rad(im(c))), c in s), s in sols)):

complex_solutions := proc(sys, {target := 64})
    local rur, bits, res, sols;
    rur := zdim_parameterization(sys, ':-verbose' = false);
    bits := 2 * target;
    do
        callpkg("Nemo", "set_precision!", Balls, bits);
        res := isolate_complex_roots_sqrfree_poly(rur[1], ':-log2eps' = -bits, ':-acb' = true);
        sols := solutions_from_rur(rur, res[1], bits);
        if bwidth(sols) <= 2^(-target) then return [sols, res[2], bits] end if;
        bits := 2 * bits;
    end do;
end proc:

# The Katsura-6 system (worksheet on real solutions), 64 complex solutions:

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## 4. Katsura-6: [target, bits of work, width of the boxes]\n");
for tg in [64, 200] do
    r := complex_solutions(katsura(6), target = tg);
    printf("  %a\n", [tg, r[3], evalf[3](bwidth(r[1]))]);
end do:

# 64 bits on the solutions need 128 bits of work; 200 bits need 800.

# ## 5. Larger systems
#
# Katsura-4 and Katsura-6, and the cyclic-5 system sum_i prod_{j<k} z_{i+j} = 0
# (k = 1, ..., 4), z_1 ... z_5 = 1:

cyclic := proc(n)
    local z, eqs, k, i, j;
    z := i -> zz || (modp(i - 1, n) + 1);
    eqs := [seq(add(mul(z(i + j), j = 0 .. k - 1), i = 1 .. n), k = 1 .. n - 1)];
    [op(eqs), mul(z(i), i = 1 .. n) - 1]
end proc:

report := proc(name, sys, vars)
    local t, r, ok;
    t := time[real]();
    r := complex_solutions(sys);
    t := time[real]() - t;
    ok := andmap(s -> andmap(p -> holds_zero(p, s, vars), sys), r[1]);
    [name, nops(r[1]), numboccur(r[2], true), nops(real_solutions(sys, prec = 32)), ok, r[3],
     evalf[3](t)]
end proc:
printf("\n## 5. [system, complex, real, real_solutions, checked, bits, seconds]\n");
printf("  %a\n", report("Katsura-4", katsura(4), [seq(u || i, i = 0 .. 4)]));
printf("  %a\n", report("Katsura-6", katsura(6), [seq(u || i, i = 0 .. 6)]));
printf("  %a\n", report("cyclic-5", cyclic(5), [seq(zz || i, i = 1 .. 5)]));

# The numbers of real solutions agree with real_solutions, and every box holds a zero of its
# system. The times include the transfers between Maple and Julia: here every division and
# every check is one call.