# Structural stability of a 2-D system: from a polydisc to real solutions
#
# A discrete 2-D system (two independent indices, for instance space and time in image
# processing) with transfer function N/D, D a polynomial in two complex variables z1, z2,
# is structurally stable when D(z1, z2) <> 0 on the closed unit bidisc
# {|z1| <= 1, |z2| <= 1}. This is a question about infinitely many points. We turn it,
# exactly, into the question "does a system of two polynomial equations have a real
# solution?", and answer it with certified tools: exact arithmetic over the rationals, and
# real roots enclosed in intervals that are guaranteed to contain them.
#
# Running example of the article:
#
#   D(z1, z2) = (2 z1^2 + 10 z1 + 12) + (z1^2 + 5 z1 + 6) z2 = (z1 + 2)(z1 + 3)(z2 + 2),
#
# stable (its zeros have |z1| = 2, 3 or |z2| = 2), which lets us check each step. The
# algorithms do not factor D.
#
# Run in batch: maple -q 03_2d_stability.mpl

interface(prettyprint = 0, quiet = true):
with(LACE):

D0 := expand((2*z1^2 + 10*z1 + 12) + (z1^2 + 5*z1 + 6)*z2):
printf("D = %a\n", sort(D0));

# ---------------------------------------------------------------------------------------
# Step 1: two 1-D tests (DeCarlo, Murray, Saeks)
#
# Theorem (DeCarlo, Murray, Saeks): D is stable if and only if D(z1, 1) and D(1, z2) have
# no root in the closed unit disc and D has no zero on the torus
# T^2 = {|z1| = 1, |z2| = 1}.
#
# A univariate p of degree d has all its roots outside the closed disc if and only if
# (1 + s)^d p((1 - s)/(1 + s)) is Hurwitz (all its roots have a negative real part): the
# Moebius map s -> (1 - s)/(1 + s) sends the open left half-plane onto the exterior of the
# disc. The exact Hurwitz test of worksheet 1 then decides; here it is Maple's
# PolynomialTools:-Hurwitz, exact on rational coefficients.

no_root_in_disc := proc(p, s)
    local d, h, j;
    d := degree(p, s);
    h := expand(add(coeff(p, s, j)*(1 - s)^j*(1 + s)^(d - j), j = 0 .. d));
    PolynomialTools:-Hurwitz(h, s);
end proc:

p1 := subs(z1 = s, z2 = 1, D0):      # D(z1, 1)
p2 := subs(z1 = 1, z2 = s, D0):      # D(1, z2)
printf("D(s, 1) = %a, outside the disc: %a\n", sort(p1), no_root_in_disc(p1, s));
printf("D(1, s) = %a, outside the disc: %a\n", sort(p2), no_root_in_disc(p2, s));

# ---------------------------------------------------------------------------------------
# Step 2: from the torus to real points (Moebius)
#
# z = (x - i)/(x + i) maps the real line onto the unit circle minus the point 1.
# Substituting it for z1 and z2 and clearing the denominators (x1 + i)^d1 (x2 + i)^d2,
# d1 and d2 the degrees of D in z1 and z2, the numerator of D is R(x1, x2) + i I(x1, x2)
# with R, I polynomials with rational coefficients, and D has no zero on the torus minus
# {z1 = 1 or z2 = 1} if and only if R = I = 0 has no real solution. (The points with a
# coordinate equal to 1 are covered by step 1.)
#
# The substitution is done term by term, z1^j z2^k becoming
# (x1 - i)^j (x1 + i)^(d1 - j) (x2 - i)^k (x2 + i)^(d2 - k); evalc splits the result into
# real and imaginary parts, x1 and x2 being real.

moebius := proc(P)
    local d1, d2, N, j, k;
    d1 := degree(P, z1);
    d2 := degree(P, z2);
    N := expand(add(add(coeff(coeff(P, z1, j), z2, k)
                        *(x1 - I)^j*(x1 + I)^(d1 - j)*(x2 - I)^k*(x2 + I)^(d2 - k),
                        k = 0 .. d2), j = 0 .. d1));
    [expand(evalc(Re(N))), expand(evalc(Im(N)))];
end proc:

Re_D, Im_D := op(moebius(D0)):
printf("R = %a\nI = %a\n", sort(Re_D), sort(Im_D));

# ---------------------------------------------------------------------------------------
# Step 3: real solutions of a bivariate system
#
# R and I have no common factor, so the system has finitely many complex solutions. For
# two equations in two unknowns, LACE has a dedicated solver, which uses subresultants over
# the integers, computed modulo many primes, and no Groebner basis. It returns a rational
# univariate representation (RUR): a polynomial rext in a new variable T (here
# T = x2 - x1), and two rational functions of T that give x1 and x2 at each of its roots.

# Over the integers Maple returns the common integer content; a constant means no common factor.
printf("gcd(R, I) = %a, of degree %d\n", gcd(Re_D, Im_D), degree(gcd(Re_D, Im_D)));

comps := biv_parameterization([Re_D, Im_D], [x1, x2], verbose = false):
printf("rext = %a, T = %a\n", sort(comps[1][1]),
       comps[1][6][1]*x1 + comps[1][6][2]*x2);
printf("factored: %a\n", factor(comps[1][1]));

# This is 324 T^4 + 369 T^2 + 100 = (36 T^2 + 25)(9 T^2 + 4): no real root. The certified
# isolation of the real solutions finds none.
#
# The bivariate RUR of LACE (no Groebner basis), then the certified real solutions: each
# component [rext, rden, num_X, num_Y, mult, lin] gives x1 = num_Y/rden, x2 = num_X/rden,
# and real_solutions returns one box (a list of intervals) per real solution. A constant
# resultant means no solution: LACE reports it as an error, caught here.

biv_real_solutions := proc(a, b)
    local cs, c;
    try
        cs := biv_parameterization([a, b], [x1, x2], verbose = false);
    catch "LACE (biv_parameterization): %1":
        if lastexception[3] <> "the resultant is a constant" then
            error;
        end if;
        cs := [];
    end try;
    [seq(op(real_solutions(c[1], c[2], [c[4], c[3]], prec = 30)), c = cs)];
end proc:

real_torus_zeros := (a, b) -> nops(biv_real_solutions(a, b)):
printf("real zeros on the torus: %d\n", real_torus_zeros(Re_D, Im_D));

# D is structurally stable. This is a proof: every step is exact or certified.

# ---------------------------------------------------------------------------------------
# Near the boundary: a tolerance cannot decide
#
# Take D_eps = 2 + z1 + z2 + eps. For eps > 0 it is stable (|z1 + z2| <= 2 on the bidisc);
# for eps < 0 it vanishes at z1 = z2 = -1 - eps/2, inside the bidisc. With eps = +-1e-12,
# a floating-point test samples |D_eps| on the torus, at the N^2 points
# (exp(2 i pi j/N), exp(2 i pi k/N)), N = 1000, in hardware floating point (evalhf):

min_on_torus := proc(e, N)
    local m, j, k, a, b, re, im;
    m := HFloat(infinity);
    for j from 0 to N - 1 do
        for k from 0 to N - 1 do
            a := 2*Pi*j/N;
            b := 2*Pi*k/N;
            re := 2.0 + cos(a) + cos(b) + e;
            im := sin(a) + sin(b);
            m := min(m, sqrt(re^2 + im^2));
        end do;
    end do;
    m;
end proc:

for e in [1e-12, -1e-12] do
    printf("eps = %.1e: min |D_eps| on the grid = %.17g\n", e, evalhf(min_on_torus(e, 1000)));
end do:

# Both minima are positive and equal (the grid contains z1 = z2 = -1): a sampling (or any
# test with a tolerance above 1e-12) cannot tell the stable case from the unstable one.
# The exact test does:

for e in [1/10^12, -1/10^12] do
    De := 2 + z1 + z2 + e;
    ab := moebius(De);
    q1 := subs(z1 = s, z2 = 1, De);
    q2 := subs(z1 = 1, z2 = s, De);
    printf("eps = %a: disc tests (%a, %a), real zeros on the torus: %d\n", e,
           no_root_in_disc(q1, s), no_root_in_disc(q2, s), real_torus_zeros(op(ab)));
end do:

# For eps = -1e-12, the two real solutions are the two zeros of D_eps on the torus, at
# distance about 1e-6 from z1 = z2 = -1: the system is not stable. The boxes of (x1, x2)
# and, from their midpoints (a floating-point illustration, not a certificate), |z1 + 1|:

sols := biv_real_solutions(op(moebius(2 + z1 + z2 - 1/10^12))):
for b in sols do
    xm := evalf[20](map(iv -> (iv[1] + iv[2])/2, b));
    printf("x1 in [%.12e, %.12e]\nx2 in [%.12e, %.12e], |z1 + 1| ~ %.3e\n",
           b[1][1], b[1][2], b[2][1], b[2][2], abs((xm[1] - I)/(xm[1] + I) + 1));
end do:

# ---------------------------------------------------------------------------------------
# Summary
#
# | Step | Question                                        | Tool                              |
# |------|-------------------------------------------------|-----------------------------------|
# | 1    | two univariate disc tests                       | exact Hurwitz test after a Moebius|
# |      |                                                 | map (PolynomialTools:-Hurwitz)    |
# | 2    | zeros on the torus <=> real zeros of {R, I}     | exact substitution                |
# | 3    | does {R = I = 0} have a real solution?          | bivariate RUR by subresultants +  |
# |      |                                                 | certified isolation (LACE)        |
#
# Sampling the torus gives no proof, and near the boundary no tolerance separates the two
# answers.