# # Real solutions of a zero-dimensional system
#
# 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.

#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

# ## 1. A first system
#
# The intersection of the cubic $y^2 = x^3 - x$ with 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]

# `real_solutions` gives the real solutions, as boxes of `BigInterval`s (one interval per
# variable, in the order of the variables) of 64 significant bits:

sols = real_solutions(sys; prec = 64)

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

rur, sep = zdim_parameterization(sys; get_separating_element = true, verbose = false)
sep

# Here $T = y$. The eliminant:

S, T = polynomial_ring(QQ, "T")
f = S(rur[1])

# Its degree, 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 $\pm\sqrt{a}$ with $a$ the real root of the cubic
# $a^3 - 6a^2 + 16a - 12$, two of them:

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 = sep)
sols200[1]

# ## 3. A check
#
# Evaluated in interval arithmetic on each box, each equation gives an interval holding 0
# (`PACE.evaluate`: AbstractAlgebra exports a function of the same name):

[has_zero.(PACE.evaluate(sys, b; precision = 128)) for b in sols]

# 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 + 2x - 2, y - 1]
rur_t, mult = zdim_parameterization(tangent; multiplicities = true, verbose = false)
mult

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

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 Gröbner basis. The solutions are the
# same:

biv_real_solutions(sys...; prec = 64)

# On two random curves of degree 16 (256 complex solutions):

using Random
Random.seed!(3)
dense(d) = sum(rand(-100:100) * x^i * y^j for i in 0:d for j in 0:d-i)
a, b = dense(16), dense(16)
t_biv = @elapsed s_biv = biv_real_solutions(a, b; prec = 64)
t_rur = @elapsed s_rur = real_solutions([a, b]; prec = 64)
(length(s_biv), t_biv), (length(s_rur), t_rur)

# Both find the same 6 real solutions among the 256 complex ones, in a few hundredths of
# a second.

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

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
katsura(2)

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

function count_real(n)
    t = @elapsed s = real_solutions(katsura(n); prec = 64)
    return (n = n, complex = 2^n, real = length(s), seconds = t)
end
[count_real(n) for n in (4, 6, 8)]