# # 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? 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) = \bigl(s^2 + a\,s + 1\bigr)^5, \qquad a = 2^{-10}.$$
# The answer is *yes*: every root has real part $-a/2 = -2^{-11} \approx -4.9\cdot 10^{-4}$.
# Every coefficient of $q$ is exactly representable in `Float64`, so no error comes
# from reading the data.

#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 LinearAlgebra, Nemo, PACE, PaceControl

a = 1 // 2^10
R, s = polynomial_ring(QQ, "s")
q = (s^2 + a * s + 1)^5
ce = Rational{BigInt}.(reverse([coeff(q, k) for k in 0:degree(q)]))   ## highest degree first
@assert Rational{BigInt}.(Float64.(ce)) == ce                          ## exact in Float64

# ## 1. Floating point
#
# The usual way: the eigenvalues of the companion matrix, which is what `roots` does in
# Matlab.

function max_real_part(coefs)
    n = length(coefs) - 1
    C = zeros(n, n)
    C[1, :] = -Float64.(coefs[2:end]) ./ Float64(coefs[1])
    for i in 2:n
        C[i, i-1] = 1
    end
    maximum(real.(eigvals(C)))
end

max_real_part(ce)

# The largest real part is positive: the computed answer is *unstable*, and it is wrong.
# A root of multiplicity $m$ is perturbed by about $\varepsilon^{1/m}$, here
# $(2^{-53})^{1/5}\approx 6\cdot 10^{-4}$, more than the distance $2^{-11}$ to the axis.
#
# Changing the algorithm does not help. The Hurwitz criterion asks that the leading
# principal minors of the Hurwitz matrix be positive. Computed on the same `Float64`
# coefficients, by elimination without pivoting, some come out negative:

function leading_minors(H)
    H = copy(H)
    m = size(H, 1)
    D = similar(H, m)
    p = one(eltype(H))
    for k in 1:m
        p *= H[k, k]
        D[k] = p
        for i in k+1:m
            f = H[i, k] / H[k, k]
            for j in k:m
                H[i, j] -= f * H[k, j]
            end
        end
    end
    D
end

D = leading_minors(PACE.hurwitz_matrix(Float64.(ce)))
(sign.(D), all(>(0), D))

# The answer is *unstable* again, and wrong: the exact minors are all positive, the last
# ones of order $10^{-43}$ and $10^{-70}$. A Routh table in `Float64` happens to answer
# correctly here, with pivots of order $10^{-27}$, but nothing tells which one to trust.

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

all(>(0), leading_minors(PACE.hurwitz_matrix(ce)))

# `hurwitz_stable` of PaceControl (Liénard-Chipart criterion, exact) agrees:

hurwitz_stable(ce)

# ## 3. Interval arithmetic with adaptive precision
#
# Each number is replaced by an interval that contains it, and each operation rounds
# outwards. 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.

function hurwitz_minors(coefs, prec)
    setprecision(BigFloat, prec) do
        H = PACE.hurwitz_matrix([BigInterval(BigFloat(x, RoundDown), BigFloat(x, RoundUp)) for x in coefs])
        m = size(H, 1)
        D = Vector{BigInterval}(undef, m)
        p = one(BigInterval)
        for k in 1:m
            p *= H[k, k]
            D[k] = p
            for i in k+1:m
                f = H[i, k] / H[k, k]
                for j in k:m
                    H[i, j] -= f * H[k, j]
                end
            end
        end
        D
    end
end

signs(D) = join(has_zero(d) ? "?" : (d > 0 ? "+" : "-") for d in D)

function decide(coefs)
    trace = String[]
    prec = 53
    while true
        D = hurwitz_minors(coefs, prec)
        push!(trace, lpad(prec, 4) * " bits: " * signs(D))
        all(d -> !has_zero(d), D) && return all(d -> d > 0, D), trace
        prec *= 2
    end
end

stable, trace = decide(ce)
trace

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

qs = squarefree_part(q)

# Its roots are then isolated in certified complex balls (PWPoly, through PACE). With
# integer coefficients (not on Windows, where PWPoly is not available):

den = lcm([denominator(coeff(qs, k)) for k in 0:degree(qs)])
ci = [BigInt(numerator(coeff(qs, k) * den)) for k in 0:degree(qs)]
rts, = isolate_complex_roots_sqrfree_poly(ci; log2eps = -80, acb = true)

#-

all(r -> Nemo.is_negative(real(r)), rts)

# Each ball contains exactly one root and lies in the open left half-plane: 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) = \bigl((s+\delta)^2+1\bigr)^4,$$
# Hurwitz again (real parts $-\delta$). Its coefficients need up to 161 bits:

p = ((s + 1 // 2^20)^2 + 1)^4
cp = Rational{BigInt}.(reverse([coeff(p, k) for k in 0:degree(p)]))
[ndigits(denominator(x), base = 2) for x in cp]

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

cr = Rational{BigInt}.(Float64.(cp))
(hurwitz_stable(cp), hurwitz_stable(cr))

# 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$.

# ## 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 $\pm 12.5\%$
# ($\delta = 1/8$):

c6 = Rational{BigInt}[1, 6, 15, 20, 15, 6, 1]
δ = 1 // 8
lo, hi = c6 .* (1 - δ), c6 .* (1 + δ)

# A first idea is to compute the Hurwitz minors with these interval coefficients, by the
# elimination of section 1. 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.

interval_signs(prec) = setprecision(BigFloat, prec) do
    signs(leading_minors(PACE.hurwitz_matrix(
        [BigInterval(BigFloat(l, RoundDown), BigFloat(h, RoundUp)) for (l, h) in zip(lo, hi)])))
end
(interval_signs(53), interval_signs(1000))

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

function kharitonov(lo, hi)
    n = length(lo)
    pats = ((0, 0, 1, 1), (0, 1, 1, 0), (1, 0, 0, 1), (1, 1, 0, 0))
    [[pat[(n - i) % 4 + 1] == 0 ? lo[i] : hi[i] for i in 1:n] for pat in pats]
end

[hurwitz_stable(k) for k in kharitonov(lo, hi)]

# The four polynomials are Hurwitz: every polynomial of the family is, which the interval
# minors could not decide. `PACE.hurwitz_stable` on interval coefficients does the same:

PACE.hurwitz_stable([BigInterval(BigFloat(l), BigFloat(h)) for (l, h) in zip(lo, hi)])

# With $\pm 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:

[hurwitz_stable(k) for k in kharitonov(c6 .* (1 - 1 // 4), c6 .* (1 + 1 // 4))]

# ## Summary
#
# | Paradigm | Answer | Guarantee |
# |---|---|---|
# | Floating point (eigenvalues, or Hurwitz minors) | unstable | none: wrong here |
# | Exact (Hurwitz minors over ℚ) | stable | proof |
# | Intervals, adaptive precision | stable | proof, or "don't know" at low precision |
# | 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.