# # Complex solutions of a zero-dimensional system
#
# 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.
#
# !!! warning "Experimental"
#     PWPoly is an experimental add-on of PACE: it is not installed with PACE. Add it once,
#     `import Pkg; Pkg.add("PWPoly")`, and load it after PACE. It is not available on
#     Windows (the DLL of PWPoly_jll exports none of the functions of PWPoly, only those of
#     FLINT it embeds): this worksheet runs on Linux and macOS.

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

# The balls are those of Arb, in Nemo; PWPoly is loaded last:

using Nemo, PACE, PWPoly

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

R, (x, y) = polynomial_ring(QQ, ["x", "y"])
sys = [y^2 - x^3 + x, x^2 + y^2 - 3]
rur = zdim_parameterization(sys; verbose = false)
f = rur[1]

# 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 (`PACE.zx_integral` clears the
# denominators). With `acb = true` the roots are balls of Arb, each holding exactly one
# root; `log2eps` sets their relative precision:

T, is_real, nb_real = isolate_complex_roots_sqrfree_poly(PACE.zx_integral(f); log2eps = -64, acb = true)
T

#-

is_real, nb_real

# Two real roots, $\pm\sqrt a$ with $a \approx 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:

function solutions_from_rur(rur, T, bits)
    f, gs = rur[1], rur[2:end]
    df = [k * f[k+1] for k in 1:length(f)-1]
    L = lcm([denominator(q) for g in vcat([df], gs) for q in g])
    den = evaluate_at_complex_points(BigInt.(L .* df), T; precision = bits)
    nums = [evaluate_at_complex_points(BigInt.(L .* g), T; precision = bits) for g in gs]
    return [[nums[k][i] / den[i] for k in eachindex(gs)] for i in eachindex(T)]
end

sols = solutions_from_rur(rur, T, 64)

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

holds_zero(p, s) = contains_zero(Nemo.evaluate(map_coefficients(c -> parent(s[1])(c), p), s))
[[holds_zero(p, s) for p in sys] for s in sols]

# The real solutions are those of `real_solutions`: their boxes meet.

meets(ball, iv) = BigFloat(midpoint(real(ball))) - BigFloat(radius(real(ball))) <= right(iv) &&
                  left(iv) <= BigFloat(midpoint(real(ball))) + BigFloat(radius(real(ball)))
reals = real_solutions(sys; prec = 64)
[any(all(meets(s[k], r[k]) for k in 1:2) for r in reals) for (s, real_root) in zip(sols, is_real) if real_root]

# ## 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, ...)`.

function complex_solutions(sys; target = 64)
    rur = zdim_parameterization(sys; verbose = false)
    bits = 2 * target
    while true
        sols, is_real = set_precision!(Balls, bits) do
            T, is_real, _ = isolate_complex_roots_sqrfree_poly(PACE.zx_integral(rur[1]);
                                                               log2eps = -bits, acb = true)
            solutions_from_rur(rur, T, bits), is_real
        end
        width = maximum(max(radius(real(c)), radius(imag(c))) for s in sols for c in s)
        width <= ldexp(1.0, -target) && return (solutions = sols, is_real = is_real, bits = bits)
        bits *= 2
    end
end

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

function katsura(n)
    K, u = polynomial_ring(QQ, ["u$i" for i in 0:n])
    U(i) = abs(i) <= n ? u[abs(i)+1] : zero(K)
    eqs = [sum(U(l) * U(m - l) for l in -n:n) - U(m) for m in 0:n-1]
    push!(eqs, sum(U(l) for l in -n:n) - 1)
    return eqs
end
width(sols) = Float64(maximum(max(radius(real(c)), radius(imag(c))) for s in sols for c in s))
[(target, (r = complex_solutions(katsura(6); target = target); (r.bits, width(r.solutions))))
 for target in (64, 200)]

# 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, \dots, 4$), $z_1 \cdots z_5 = 1$:

function cyclic(n)
    K, z = polynomial_ring(QQ, ["z$i" for i in 1:n])
    eqs = [sum(prod(z[mod1(i + j, n)] for j in 0:k-1) for i in 1:n) for k in 1:n-1]
    push!(eqs, prod(z) - 1)
    return eqs
end

function report(name, sys)
    t = @elapsed r = complex_solutions(sys)
    ok = all(holds_zero(p, s) for s in r.solutions for p in sys)
    return (system = name, complex = length(r.solutions), real = count(r.is_real),
            real_solutions = length(real_solutions(sys; prec = 32)), checked = ok,
            bits = r.bits, seconds = round(t, digits = 2))
end
[report(name, s) for (name, s) in ("Katsura-4" => katsura(4), "Katsura-6" => katsura(6),
                                     "cyclic-5" => cyclic(5))]

# The numbers of real solutions agree with `real_solutions`, and every box holds a zero of
# its system. Katsura-8, with 256 complex solutions (84 real), takes about 30 s.