# One question, four paradigms: is this polynomial Hurwitz?
#
# A characteristic polynomial with exact coefficients is given. Are all its roots in the
# open left half-plane? (A polynomial with this property is called Hurwitz, or Hurwitz
# stable.) We answer with floating-point arithmetic, exact arithmetic, interval
# arithmetic and certified root isolation, and compare.
#
# The example is built to be hard for floating point: five times the same pair of
# roots, close to the imaginary axis,
#     q(s) = (s^2 + a*s + 1)^5,   a = 2^(-10).
# The answer is *yes*: every root has real part -a/2 = -2^(-11) ~ -4.9e-4.
# Every coefficient of q is exactly representable in Float64 (IEEE double precision,
# 53 significant bits), so no error comes from reading the data.
#
# Maple version of the Julia worksheet 01_one_question_four_paradigms.jl.
# Run in batch with:  maple -q 01_one_question_four_paradigms.mpl

interface(rtablesize = 20, displayprecision = -1):
DynamicSystems:-RouthTable(s + 1, s):   # load it before q is assigned (q is one of its names)
with(LACE):

a := 1/2^10:
q := expand((s^2 + a*s + 1)^5):
n := degree(q, s):
ce := [seq(coeff(q, s, n - k), k = 0 .. n)]:   # highest degree first
cf := Vector(ce, datatype = float[8]):          # the same coefficients as Float64
if map(convert, convert(cf, list), rational, exact) <> ce then error "not exact in Float64" end if:
printf("q(s) = %a\n", sort(q, s));
printf("coefficients exact in Float64: true\n");

# ## 1. Floating point
#
# The usual way: the eigenvalues of the companion matrix, which is what `roots` does in
# Matlab. The companion matrix of s^n + c_1 s^(n-1) + ... + c_n has first row
# (-c_1, ..., -c_n) and ones on its subdiagonal; its eigenvalues are the roots.
# Maple's LinearAlgebra:-Eigenvalues on a float[8] matrix calls LAPACK, in Float64.

MaxRealPart := proc(c::Vector)
    local m, C, i, j;
    m := LinearAlgebra:-Dimension(c) - 1;
    C := Matrix(m, m, datatype = float[8]);
    for j to m do C[1, j] := -c[j + 1]/c[1] end do;
    for i from 2 to m do C[i, i - 1] := 1 end do;
    max(map(Re, convert(LinearAlgebra:-Eigenvalues(C), list)))
end proc:

printf("\n## 1. Floating point\n");
printf("largest real part of the eigenvalues (Float64): %.6e\n", MaxRealPart(cf));

# The largest real part is positive: the computed answer is *unstable*, and it is wrong.
# A root of multiplicity m is perturbed by about eps^(1/m), here
# (2^(-53))^(1/5) ~ 6e-4, more than the distance 2^(-11) to the axis.
#
# Changing the algorithm does not help. The Hurwitz criterion (Hurwitz, 1895): with a
# positive leading coefficient, a polynomial is Hurwitz if and only if the leading
# principal minors of its Hurwitz matrix are all positive. The Hurwitz matrix of
# c_0 s^n + c_1 s^(n-1) + ... + c_n is the n x n matrix of entries c_(2i-j), zero
# outside 0..n; a leading principal minor is the determinant of its upper left k x k
# block. We compute the minors on the same Float64 coefficients, by elimination without
# pivoting: the k-th minor is the product of the first k pivots. The procedure runs
# under evalhf, that is in hardware Float64 arithmetic.

HurwitzMatrix := proc(c::list)
    local m, H, i, j, idx;
    m := nops(c) - 1;
    H := Matrix(m, m);
    for i from 0 to m - 1 do
        for j from 0 to m - 1 do
            idx := 2*i + 1 - j;
            if 0 <= idx and idx <= m then H[i + 1, j + 1] := c[idx + 1] end if
        end do
    end do;
    H
end proc:

# Leading principal minors by elimination without pivoting; H is overwritten, D receives them.
LeadingMinors := proc(H, D, m)
    local k, i, j, f, p;
    p := 1;
    for k to m do
        p := p*H[k, k];
        D[k] := p;
        for i from k + 1 to m do
            f := H[i, k]/H[k, k];
            for j from k to m do H[i, j] := H[i, j] - f*H[k, j] end do
        end do
    end do;
    0
end proc:

Hf := Matrix(HurwitzMatrix(ce), datatype = float[8]):
Df := Vector(n, datatype = float[8]):
evalhf(LeadingMinors(var(Hf), var(Df), n)):
printf("signs of the Hurwitz minors (Float64): %a\n", map(signum, convert(Df, list)));
printf("all positive: %a\n", evalb(min(convert(Df, list)) > 0));

# The answer is *unstable* again, and wrong: the exact minors (section 2) are all
# positive, the last ones of order 1e-43 and 1e-70.
#
# The Routh table is another form of the same criterion: its first column must have
# a constant sign. Computed in Float64 it happens to answer correctly here, with pivots
# of order 1e-27, but nothing tells which floating-point answer to trust.

# First column of the Routh table; c highest degree first, r receives the column.
RouthColumn := proc(c, r, m)
    local A, B, T, k, j, w;
    w := iquo(m, 2) + 2;
    A := Array(1 .. w, datatype = float[8]);
    B := Array(1 .. w, datatype = float[8]);
    T := Array(1 .. w, datatype = float[8]);
    for j to w do
        if 2*j - 1 <= m + 1 then A[j] := c[2*j - 1] end if;
        if 2*j <= m + 1 then B[j] := c[2*j] end if
    end do;
    r[1] := A[1];
    r[2] := B[1];
    for k from 3 to m + 1 do
        for j to w - 1 do T[j] := (B[1]*A[j + 1] - A[1]*B[j + 1])/B[1] end do;
        T[w] := 0;
        for j to w do A[j] := B[j]; B[j] := T[j] end do;
        r[k] := B[1]
    end do;
    0
end proc:

rf := Vector(n + 1, datatype = float[8]):
evalhf(RouthColumn(cf, var(rf), n)):
printf("first column of the Routh table (Float64):\n  %s\n",
       StringTools:-Join(map(x -> sprintf("%.2e", x), convert(rf, list)), " "));
printf("all positive: %a\n", evalb(min(convert(rf, list)) > 0));

# Maple's own commands behave differently, and it is worth knowing why.
# PolynomialTools:-Hurwitz and DynamicSystems:-RouthTable accept floating-point
# coefficients, but they first convert each float to the rational number it denotes,
# then compute exactly. On the Float64 coefficients above, they are exact and answer
# correctly (*stable*), with an exact last Routh pivot of order 1e-27.

qf := add(cf[k + 1]*s^(n - k), k = 0 .. n):
printf("PolynomialTools:-Hurwitz on the Float64 coefficients: %a\n",
       PolynomialTools:-Hurwitz(qf, s));
rt := DynamicSystems:-RouthTable(qf, s):
printf("signs of the first column of RouthTable on them: %a\n",
       map(signum, convert(rt[.., 1], list)));

# The usual Maple way, evalf(q), rounds each coefficient to Digits = 10 *decimal* digits.
# This changes the polynomial, and the rounded polynomial is no longer Hurwitz: both
# commands answer *unstable*, exactly, on the wrong data. Section 5 comes back to this.

qe := evalf(q):
printf("PolynomialTools:-Hurwitz(evalf(q), s) with Digits = 10: %a\n",
       PolynomialTools:-Hurwitz(qe, s));
rt := DynamicSystems:-RouthTable(qe, s):
printf("signs of the first column of RouthTable(evalf(q), s): %a\n",
       map(signum, convert(rt[.., 1], list)));

# ## 2. Exact arithmetic
#
# The same criterion on rational numbers decides without any tolerance:

He := HurwitzMatrix(ce):
De := Vector(n):
LeadingMinors(He, De, n):
printf("\n## 2. Exact arithmetic\n");
printf("exact Hurwitz minors:\n  %s\n",
       StringTools:-Join(map(x -> sprintf("%.2e", x), convert(De, list)), " "));
printf("all positive: %a\n", evalb(min(convert(De, list)) > 0));

# PolynomialTools:-Hurwitz of Maple (a continued fraction expansion, exact) agrees:

printf("PolynomialTools:-Hurwitz(q, s): %a\n", PolynomialTools:-Hurwitz(q, s));

# ## 3. Interval arithmetic with adaptive precision
#
# Each number is replaced by an interval that contains it, and each operation rounds
# outwards: the lower end down, the upper end up. A sign is decided when the interval
# does not contain 0. We compute the leading principal minors of the Hurwitz matrix by
# elimination without pivoting, as above.
#
# The intervals below emulate MPFI with a working precision of prec bits: each end is a
# binary number of prec significant bits; each operation is done exactly on the ends,
# then rounded outwards to prec bits. Once a pivot contains 0, the remaining minors
# (multiples of it) contain 0 too and stay undecided.

# x rounded to prec significant bits, down (dir = -1) or up (dir = 1).
RoundBits := proc(x::rational, prec::posint, dir)
    local e, sc;
    if x = 0 then return 0 end if;
    e := ilog2(abs(numer(x))) - ilog2(denom(x));
    if 2^e > abs(x) then e := e - 1 end if;
    sc := 2^(prec - 1 - e);
    if dir < 0 then floor(x*sc)/sc else ceil(x*sc)/sc end if
end proc:

IvOut := (l, u, prec) -> [RoundBits(l, prec, -1), RoundBits(u, prec, 1)]:
IvSub := (x, y, prec) -> IvOut(x[1] - y[2], x[2] - y[1], prec):
IvMul := proc(x, y, prec)
    local P;
    P := [x[1]*y[1], x[1]*y[2], x[2]*y[1], x[2]*y[2]];
    IvOut(min(P), max(P), prec)
end proc:
IvHasZero := x -> evalb(x[1] <= 0 and 0 <= x[2]):
IvSign := x -> `if`(IvHasZero(x), "?", `if`(x[1] > 0, "+", "-")):

# Quotient by an interval y not containing 0, rounded once at the end.
IvDiv := proc(x, y, prec)
    local P;
    P := [x[1]/y[1], x[1]/y[2], x[2]/y[1], x[2]/y[2]];
    IvOut(min(P), max(P), prec)
end proc:

#  The minors for coefficients given as intervals ci = [[l, u], ...] (section 6), or as
#  numbers c (each the interval [x, x]).
HurwitzMinorsIv := proc(c::list, prec::posint)
    HurwitzMinorsIvData(map(x -> [x, x], c), prec)
end proc:

HurwitzMinorsIvData := proc(ci::list, prec::posint)
    local m, H, D, p, k, i, j, f;
    m := nops(ci) - 1;
    H := map(k -> `if`(k = 0, [0, 0], IvOut(ci[k][1], ci[k][2], prec)),
             HurwitzMatrix([seq(1 .. m + 1)]));
    D := Array(1 .. m, fill = [-infinity, infinity]);
    p := [1, 1];
    for k to m do
        if IvHasZero(H[k, k]) then break end if;
        p := IvMul(p, H[k, k], prec);
        D[k] := p;
        for i from k + 1 to m do
            f := IvDiv(H[i, k], H[k, k], prec);
            for j from k to m do H[i, j] := IvSub(H[i, j], IvMul(f, H[k, j], prec), prec) end do
        end do
    end do;
    convert(D, list)
end proc:

Decide := proc(c::list)
    local prec, D, tr;
    tr := [];
    prec := 53;
    do
        D := HurwitzMinorsIv(c, prec);
        tr := [op(tr), sprintf("%4d bits: %s", prec, cat(op(map(IvSign, D))))];
        if not ormap(IvHasZero, D) then
            return andmap(d -> d[1] > 0, D), tr
        end if;
        prec := 2*prec
    end do
end proc:

printf("\n## 3. Interval arithmetic with adaptive precision\n");
stable, tr := Decide(ce):
for t in tr do printf("%s\n", t) end do;
printf("stable: %a\n", stable);

# At low precision some signs are undecided ("?"): intervals do not give a wrong
# answer, they say they do not know. Doubling the precision until every sign is
# decided gives the right answer, at a cost that adapts to the difficulty of the input.

# ## 4. Certified root isolation
#
# The question can also be answered by locating the roots. Exact arithmetic first
# removes the multiplicities: the square-free part q/gcd(q, q') has the same roots,
# each simple.

printf("\n## 4. Certified root isolation\n");
qs := quo(q, gcd(q, diff(q, s)), s):
printf("square-free part: %a\n", qs);

# Maple has no certified isolator of complex roots in balls, such as PWPoly used by the
# Julia version. We isolate them as real solutions of a system instead: s = x + i y is a
# root of qs if and only if the real part R(x, y) and the imaginary part I(x, y) of
# qs(x + i y) both vanish. LACE:-real_solutions computes a rational univariate
# representation of {R = 0, I = 0} and returns one box [x interval, y interval] per
# solution, with rational ends, each containing exactly one solution.

z := evalc(subs(s = x + I*y, qs)):
Rxy := evalc(Re(z)):
Ixy := evalc(Im(z)):
printf("R = %a,  I = %a\n", Rxy, Ixy);
B := LACE:-real_solutions([Rxy, Ixy], [x, y], prec = 80):
for b in B do
    printf("x in [%.22e, %.22e],  y in [%.22e, %.22e]\n",
           evalf[25](b[1][1]), evalf[25](b[1][2]), evalf[25](b[2][1]), evalf[25](b[2][2]))
end do;
printf("number of roots: %d, all real parts negative: %a\n",
       nops(B), andmap(b -> evalb(b[1][2] < 0), B));

# Each box contains exactly one root and lies in the open left half-plane (the upper end
# of its x interval, a rational number, is negative): the polynomial is Hurwitz, and we
# also know where its roots are, and with which multiplicity.

# ## 5. When the data are not floating-point numbers
#
# Take now delta = 2^(-20) and
#     p(s) = ((s + delta)^2 + 1)^4,
# Hurwitz again (real parts -delta). Its coefficients need up to 161 bits:

printf("\n## 5. When the data are not floating-point numbers\n");
p := expand(((s + 1/2^20)^2 + 1)^4):
cp := [seq(coeff(p, s, 8 - k), k = 0 .. 8)]:
printf("bits of the denominators: %a\n", map(x -> ilog2(denom(x)) + 1, cp));

# Reading them as Float64 changes the polynomial, and the exact test shows that the
# rounded polynomial is *not* Hurwitz:

cr := map(convert, convert(Vector(cp, datatype = float[8]), list), rational, exact):
pr := add(cr[k + 1]*s^(8 - k), k = 0 .. 8):
printf("PolynomialTools:-Hurwitz: p %a, p rounded to Float64 %a\n",
       PolynomialTools:-Hurwitz(p, s), PolynomialTools:-Hurwitz(pr, s));

# Here no floating-point algorithm, however good, can give the right answer: the
# question is lost before any computation. This is why the inputs of the exact and
# certified tools are rational numbers, and why a decimal coefficient such as 0.1
# should be read as 1/10. In Maple, evalf(q) in section 1 was the same trap.

# ## 6. When the data are not known exactly
#
# Coefficients coming from measurements are only known to lie in intervals. The question
# becomes: is *every* polynomial of the family Hurwitz? Take (s+1)^6, whose coefficients
# are 1, 6, 15, 20, 15, 6, 1, and let each of them vary independently by +-12.5 %
# (delta = 1/8):

printf("\n## 6. When the data are not known exactly\n");
c6 := [1, 6, 15, 20, 15, 6, 1]:
clo := map(x -> x*(1 - 1/8), c6):
chi := map(x -> x*(1 + 1/8), c6):

# A first idea is to compute the Hurwitz minors with these interval coefficients, by the
# elimination of section 3. If all of them are positive intervals, every polynomial of
# the family is Hurwitz. Here some signs stay undecided, and raising the precision does
# not help: the width comes from the data, and interval arithmetic overestimates it.

for prec in [53, 1000] do
    printf("%4d bits: %s\n", prec,
           cat(op(map(IvSign, HurwitzMinorsIvData(zip((l, h) -> [l, h], clo, chi), prec)))))
end do;

# Kharitonov's theorem (1978): when the coefficients vary independently and the interval
# of the leading coefficient does not contain 0, every polynomial of the family is
# Hurwitz if and only if four particular polynomials are. Their coefficients are bounds of
# the intervals: from the constant term upwards, lower (l) or upper (u) bounds in the
# patterns l l u u, l u u l, u l l u and u u l l, repeated every four powers. With
# rational bounds, they are tested exactly.

Kharitonov := proc(l::list, u::list)
    local n, pat, i;
    n := nops(l);
    [seq([seq(`if`(pat[irem(n - i, 4) + 1] = 0, l[i], u[i]), i = 1 .. n)],
         pat in [[0, 0, 1, 1], [0, 1, 1, 0], [1, 0, 0, 1], [1, 1, 0, 0]])]
end proc:

HurwitzOf := proc(c)
    local k;
    PolynomialTools:-Hurwitz(add(c[k]*s^(nops(c) - k), k = 1 .. nops(c)), s)
end proc:
printf("the four Kharitonov polynomials Hurwitz: %a\n", map(HurwitzOf, Kharitonov(clo, chi)));

# The four polynomials are Hurwitz: every polynomial of the family is, which the interval
# minors could not decide. (PACE.hurwitz_stable of the Julia worksheet, which does the same
# on interval coefficients, has no Maple counterpart.)
#
# With +-25 % (delta = 1/4), two of the four polynomials are not Hurwitz. They belong to
# the family, so the family is not robustly stable, and this is again a proof:

printf("with delta = 1/4: %a\n", map(HurwitzOf, Kharitonov(map(x -> x*(1 - 1/4), c6),
                                                           map(x -> x*(1 + 1/4), c6))));

# ## Summary
#
# | Paradigm                                         | Answer                | Guarantee        |
# |--------------------------------------------------|-----------------------|------------------|
# | Floating point (eigenvalues, or Hurwitz minors)  | unstable              | none: wrong here |
# | Exact (Hurwitz minors over Q, Hurwitz of Maple)  | stable                | proof            |
# | Intervals, adaptive precision                    | stable                | proof, or        |
# |                                                  |                       | "don't know"     |
# | Exact square-free part + certified isolation     | stable, roots located | proof            |
#
# And when the data are rounded before the computation, even an exact method answers
# another question. When they are only known to lie in intervals, the minors computed
# in interval arithmetic may not decide, and Kharitonov's theorem reduces the question to
# four exact tests.

LACE:-finish():