Complex 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 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.

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]

7-element Vector{Rational{BigInt}}:
 -12
   0
  16
   0
  -6
   0
   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

6-element Vector{Nemo.ComplexFieldElem}:
 [1.0734488352317652300 +/- 4.41e-20]
 [-1.0734488352317652300 +/- 4.41e-20]
 [-1.6809120138921139586 +/- 4.45e-21] + [0.63372817467511970984 +/- 1.52e-21]*im
 [-1.6809120138921139586 +/- 4.45e-21] + [-0.63372817467511970984 +/- 1.52e-21]*im
 [1.6809120138921139586 +/- 4.45e-21] + [0.63372817467511970984 +/- 1.52e-21]*im
 [1.6809120138921139586 +/- 4.45e-21] + [-0.63372817467511970984 +/- 1.52e-21]*im

is_real, nb_real

(Bool[1, 1, 0, 0, 0, 0], 2)

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)

6-element Vector{Vector{Nemo.ComplexFieldElem}}:
 [[1.3593040859718 +/- 3.63e-14] + [+/- 2.67e-39]*im, [1.07344883523177 +/- 6.71e-15] + [+/- 3.34e-39]*im]
 [[1.3593040859718 +/- 3.63e-14] + [+/- 2.67e-39]*im, [-1.07344883523177 +/- 6.71e-15] + [+/- 3.34e-39]*im]
 [[-1.179652042985888210 +/- 4.97e-19] + [-0.903013145857004187 +/- 6.61e-19]*im, [-1.680912013892113958 +/- 7.81e-19] + [0.633728174675119710 +/- 3.41e-19]*im]
 [[-1.179652042985888210 +/- 4.97e-19] + [0.903013145857004187 +/- 6.61e-19]*im, [-1.680912013892113958 +/- 7.81e-19] + [-0.633728174675119710 +/- 3.41e-19]*im]
 [[-1.179652042985888210 +/- 4.97e-19] + [0.903013145857004187 +/- 6.61e-19]*im, [1.680912013892113958 +/- 7.81e-19] + [0.633728174675119710 +/- 3.41e-19]*im]
 [[-1.179652042985888210 +/- 4.97e-19] + [-0.903013145857004187 +/- 6.61e-19]*im, [1.680912013892113958 +/- 7.81e-19] + [-0.633728174675119710 +/- 3.41e-19]*im]

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]

6-element Vector{Vector{Bool}}:
 [1, 1]
 [1, 1]
 [1, 1]
 [1, 1]
 [1, 1]
 [1, 1]

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]

2-element Vector{Bool}:
 1
 1

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

complex_solutions (generic function with 1 method)

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

2-element Vector{Tuple{Int64, Tuple{Int64, Float64}}}:
 (64, (128, 3.0790341769646653e-27))
 (200, (800, 3.4027573530778945e-131))

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

3-element Vector{@NamedTuple{system::String, complex::Int64, real::Int64, real_solutions::Int64, checked::Bool, bits::Int64, seconds::Float64}}:
 (system = "Katsura-4", complex = 16, real = 12, real_solutions = 12, checked = 1, bits = 128, seconds = 0.15)
 (system = "Katsura-6", complex = 64, real = 32, real_solutions = 32, checked = 1, bits = 128, seconds = 2.34)
 (system = "cyclic-5", complex = 70, real = 10, real_solutions = 10, checked = 1, bits = 128, seconds = 0.72)

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.