Real solutions of a zero-dimensional system

Download: Jupyter notebook, Julia script, Maple script (PACE called from Maple through UDX, see Using PACE.jl in Maple). To install PACE and run them, see Running the worksheets.

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.

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]

2-element Vector{AbstractAlgebra.Generic.MPoly{Rational{BigInt}}}:
 -x^3 + x + y^2
 x^2 + y^2 - 3

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

sols = real_solutions(sys; prec = 64)

2-element Vector{Vector{BigInterval}}:
 [[1.359304085971776420730660392788800573278, 1.359304085971776420730660392812266330097], [-1.07344883523176522996096257952156660135, -1.073448835231765229960962579519779626699]]
 [[1.359304085971776420730660392788800573278, 1.359304085971776420730660392812266330097], [1.073448835231765229960962579519779626699, 1.07344883523176522996096257952156660135]]

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

2-element Vector{Int64}:
 0
 1

Here $T = y$. The eliminant:

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

T^6 - 6*T^4 + 16*T^2 - 12

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

2-element Vector{LACE.Root}:
 Root(-2.0 .. -1.0)
 Root(1.0 .. 2.0)

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]

2-element Vector{BigInterval}:
   [1.359304085971776420730660392800532555346481278067672283797141251583875588944174, 1.359304085971776420730660392800532555346481278067672283797141251583875588945072]
 [-1.073448835231765229960962579520667916192623147582044134487669184447604210904965, -1.073448835231765229960962579520667916192623147582044134487669184447604210904895]

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]

2-element Vector{BitVector}:
 [1, 1]
 [1, 1]

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

1-element Vector{Pair{Int64, Vector{Rational{BigInt}}}}:
 2 => [-1, 1]

One solution of multiplicity 2, and it is real:

real_solutions(tangent; prec = 64)

1-element Vector{Vector{BigInterval}}:
 [[1.0, 1.0], [1.0, 1.0]]

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)

2-element Vector{Vector{Any}}:
 [[1.359304085971776420730660392793470116071, 1.359304085971776420730660392807679904994], [-1.073448835231765229960962579528413282514, -1.073448835231765229960962579512944864983]]
 [[1.359304085971776420730660392793470116071, 1.359304085971776420730660392807679904994], [1.073448835231765229960962579512944864983, 1.073448835231765229960962579528413282514]]

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)

((6, 0.038581247), (6, 0.074505041))

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)

3-element Vector{AbstractAlgebra.Generic.MPoly{Rational{BigInt}}}:
 u0^2 - u0 + 2*u1^2 + 2*u2^2
 2*u0*u1 + 2*u1*u2 - u1
 u0 + 2*u1 + 2*u2 - 1

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

3-element Vector{@NamedTuple{n::Int64, complex::Int64, real::Int64, seconds::Float64}}:
 (n = 4, complex = 16, real = 12, seconds = 0.006309466)
 (n = 6, complex = 64, real = 32, seconds = 0.029267179)
 (n = 8, complex = 256, real = 84, seconds = 0.554012029)