diff --git a/src/AlgebraicSolving.jl b/src/AlgebraicSolving.jl index 94f4046..7e313a8 100644 --- a/src/AlgebraicSolving.jl +++ b/src/AlgebraicSolving.jl @@ -22,6 +22,9 @@ include("algorithms/param-ideal/groebner-bases.jl") include("algorithms/param-ideal/multiplication-matrices.jl") include("algorithms/param-ideal/hermite-matrices.jl") include("algorithms/param-ideal/parametrizations.jl") +include("algorithms/param-ideal/sign-determination.jl") +include("algorithms/param-ideal/semialgebraic-set.jl") +include("algorithms/param-ideal/real-root-classification.jl") #= siggb =# include("siggb/siggb.jl") #= progress =# diff --git a/src/algorithms/param-ideal/real-root-classification.jl b/src/algorithms/param-ideal/real-root-classification.jl new file mode 100644 index 0000000..4b40d11 --- /dev/null +++ b/src/algorithms/param-ideal/real-root-classification.jl @@ -0,0 +1,292 @@ +@doc Markdown.doc""" + __pick_points(points::Vector{Vector{Vector{QQFieldElem}}}) -> Vector{Vector{QQFieldElem}} + +Returns a list of sample points, one from each isolating box in `points`. + +**Note**: This is an internal function. +""" +function _pick_points(points::Vector{Vector{Vector{QQFieldElem}}})::Vector{Vector{QQFieldElem}} + return map(point -> map(x -> (x[1] + x[2]) / 2, point), points) +end + +@doc Markdown.doc""" + __unique_by_sign(points::Vector{Vector{QQFieldElem}}, f::Vector{QQMPolyRingElem}) -> Vector{Vector{QQFieldElem}} + +Returns a list of points from `points` such that the sign vectors of `f` at these points are unique. + +**Note**: This is an internal function. +""" +function _unique_by_sign(points::Vector{Vector{QQFieldElem}}, f::Vector{QQMPolyRingElem})::Vector{Vector{QQFieldElem}} + signs = Set{Vector{QQFieldElem}}() + res = Vector{Vector{QQFieldElem}}() + for point in points + s = [sign(evaluate(q, point)) for q in f] + if !(s in signs) + push!(signs, s) + push!(res, point) + end + end + res +end + +@doc Markdown.doc""" + _identity_matrix(m::Int) -> Matrix{Int} + +Returns the identity matrix of size `m x m`. + +**Note**: This is an internal function. +""" +function _identity_matrix(m::Int)::Matrix{Int} + I = zeros(Int, m, m) + for i in 1:m + I[i, i] = 1 + end + return I +end + +@doc Markdown.doc""" + _random_matrix(m::Int, n::Int) -> Matrix{Int} + +Returns a random integer matrix of size `m x n` with entries in the range [1, 99]. + +**Note**: This is an internal function. +""" +function _random_matrix(m::Int, n::Int)::Matrix{Int} + map(a -> 1 + abs(a) % 99, rand(Int, m, n)) +end + +@doc Markdown.doc""" + real_root_classification(I::ParametricIdeal{K}, g::Vector{MPoly{K}}; ) -> SemialgebraicSet + +Computes the real root classification of the polynomials `g` with respect to the parametric ideal `I`. The result is a semi-algebraic set that describes the regions in the parameter space where the number of real roots such that `g` is positive is constant. Each region is accompanied by the constant number of real roots in that region, and a witness point in that region. + +**Note**: The regions described by the semi-algebraic set may differ from the actual regions due to the fact that the Hermite matrices may not specialize well at certain parameter values. However, their symmetric difference is guaranteed to be contained in a proper algebraic set. + +# Arguments +- `I::ParametricIdeal{K}`: input parametric ideal. +- `g::Vector{MPoly{K}}`: input list of polynomials. +- `nr_thrds::Int=1`: the number of threads to use for parallel computations. +- `worker_pool::AbstractWorkerPool=default_worker_pool()`: the worker pool to use for parallel computations. +- `info_level::Int=0`: info level printout: off (`0`, default), summary (`1`), detailed (`2`). +- `show_progress::Bool=info_level >= 1`: whether to show progress bars during computations. +- `output_form::Symbol=:sign`: the output form of the real root classification, can be either polynomials accompanied by their signs (`:sign`), Hermite matrices accompanied by their signatures (`:signature`), or Hermite matrices and multiplication matrices (`:matrix`). +- `ignore_no_real_roots::Bool=false`: whether to ignore sign conditions with no real roots. +""" +function real_root_classification( + I::ParametricIdeal{K} where {K<:FracFieldElem}, + g::Vector{MPoly{K}} where {K<:FracFieldElem}; + nr_thrds::Int=1, + worker_pool::AbstractWorkerPool=default_worker_pool(), + info_level::Int=0, + show_progress::Bool=info_level >= 1, + output_form::Symbol=:sign, + ignore_no_real_roots::Bool=false, +)::SemialgebraicSet + R = base_ring(I.base_frac_field) + # We first handle the case where a certain Hermite matrix is singular. + # This can be handle by ideal saturation/radicalization. + α = fill(0, length(g)) + H_1 = hermite_matrix(g, α, I; nr_thrds=nr_thrds, worker_pool=worker_pool, show_progress=show_progress) + if size(H_1, 1) == 0 + return Empty() + end + h_1 = det(H_1) + if is_zero(h_1) + @warn "The Hermite matrix associated to one has zero determinant." + @warn "Trying again with the radical ideal." + I.radicalize = true + return real_root_classification(I, g; nr_thrds=nr_thrds, worker_pool=worker_pool, show_progress=show_progress) + end + M = Vector{MatSpaceElem{FracFieldElem{QQMPolyRingElem}}}() + m = Vector{FracFieldElem{QQMPolyRingElem}}() + for i in 1:length(g) + M_α = multiplication_matrix(I, g[i]; nr_thrds=nr_thrds, worker_pool=worker_pool, show_progress=show_progress) + if M_α == zero_matrix(QQ, size(M_α, 1), size(M_α, 2)) && i == 1 + @warn "The multiplication matrix associated to the first polynomial is zero." + @warn "This indicates an additional rank deficiency." + @warn "Trying again by removing the first polynomial." + return real_root_classification(I, g[2:end]; nr_thrds=nr_thrds, worker_pool=worker_pool, show_progress=show_progress) + end + m_α = det(M_α) + if is_zero(m_α) + @warn "The multiplication matrix associated to one polynomial has zero determinant." + @warn "Trying again with the saturation ideal." + push!(I.sats, g) + return real_root_classification(I, g; nr_thrds=nr_thrds, worker_pool=worker_pool, show_progress=show_progress) + end + push!(M, M_α) + push!(m, m_α) + end + if output_form == :matrix + return MatEnum(hm=H_1, mm=M) + end + # We now compute a polynomial w whose vanishing defines a proper algebraic set of parameters, + # where at least one of the Hermite matrices is singular or does not specialize well. + # While this does not guarantee that the Hermite matrices specialize well outside of V(w), + # the problematic parameters are contained in a proper algebraic set, + # and thus the real root classification is generically correct. + w = Vector{QQMPolyRingElem}() + if is_empty(I.gens_alt) + push!(w, _change_ring(denominator(h_1), R)) + push!(w, _change_ring(numerator(h_1), R)) + for m_α in m + push!(w, _change_ring(denominator(m_α), R)) + push!(w, _change_ring(numerator(m_α), R)) + end + else + if info_level >= 2 + println("Computing GCD of determinants of Hermite matrices associated to the original and alternative generators.") + println("Old factors: $(factor(numerator(h_1)))") + end + I.gens, I.gens_alt = I.gens_alt, I.gens + H_1_alt = hermite_matrix(g, α, I; nr_thrds=nr_thrds, worker_pool=worker_pool, show_progress=show_progress) + I.gens, I.gens_alt = I.gens_alt, I.gens + h_1_alt = det(H_1_alt) + push!(w, _change_ring(denominator(h_1), R)) + push!(w, _change_ring(denominator(h_1_alt), R)) + push!(w, _change_ring(gcd(numerator(h_1), numerator(h_1_alt)), R)) + if info_level >= 2 + println("New factors: $(factor(numerator(h_1_alt)))") + end + if info_level >= 1 + println("GCD of determinants of Hermite matrices: $(factor(w[end]))") + end + for m_α in m + push!(w, _change_ring(denominator(m_α), R)) + end + end + # We now compute at least one point in each connected component of the semi-algebraic set + # defined by the non-vanishing of the polynomials in w. + # This allows us to determine the sign conditions of g at each connected component, + # and thus to determine the adapted family of Hermite matrices that we need to consider. + p = _pick_points(points_per_components(QQMPolyRingElem[], QQMPolyRingElem[], w; nr_thrds=nr_thrds, worker_pool=worker_pool, info_level=info_level)) + if info_level >= 1 + println("Number of sample points: ", length(p)) + end + prog = Progress.ProgressBar(total=length(p); desc="Sign conditions", enabled=show_progress) + Progress.update!(prog, 0) + Σ = Vector{Vector{Vector{Int}}}([]) + c = Vector{Vector{Int}}([]) + for pᵢ in p + Σᵢ, cᵢ = sign_determination(I, g, pᵢ) + push!(Σ, Σᵢ) + push!(c, cᵢ) + Progress.next!(prog) + end + Progress.finish!(prog) + if info_level >= 1 + println("Sign conditions at sample points: ", Set(Σ)) + end + σ = fill(1, length(g)) + if !(σ in vcat(Σ...)) + if info_level >= 1 + println("Sign $σ not found at sample points; returning empty set.") + end + return Empty() + end + if info_level >= 1 + printstyled("[ Good news: ", bold=true, color=:light_magenta) + println("sign $σ found at sample points!") + end + A = Vector{Vector{Int}}([]) + for i in 1:length(Σ) + for j in i:length(Σ) + append!(A, adapted_family(sort_signs(unique([Σ[i]; Σ[j]])))) + end + end + A = sort_indices(unique(A)) + if info_level >= 1 + println("Adapted family for real root classification: ", A) + end + H = Vector{MatSpaceElem{FracFieldElem{QQMPolyRingElem}}}() + for α in A + H_α = hermite_matrix(g, α, I; nr_thrds=nr_thrds, worker_pool=worker_pool, show_progress=show_progress) + push!(H, H_α) + end + if output_form == :signature + sigs = Vector{Vector{Int}}() + counts = Vector{Int}() + witnesses = Vector{Vector{QQFieldElem}}() + for (pᵢ, Σᵢ, cᵢ) in zip(p, Σ, c) + j = findfirst(==(σ), Σᵢ) + cᵢⱼ = isnothing(j) ? 0 : cᵢ[j] + if ignore_no_real_roots && cᵢⱼ == 0 + continue + end + sig = [tarski_query(g, α, I, pᵢ; nr_thrds=nr_thrds, worker_pool=worker_pool, show_progress=show_progress) for α in A] + if (!(sig in sigs)) + push!(sigs, sig) + push!(counts, cᵢⱼ) + push!(witnesses, pᵢ) + end + end + return SigEnum(hms=H, sigs=sigs, counts=counts, witnesses=witnesses) + end + # We now compute the real root classification defined with semi-algebraic formulas. + # The polynomials appearing in the semi-algebraic descriptions are leading principal minors + # of the previously-identified Hermite matrices. + # To determine their possible signs, another round of sample points is needed. + Q = parent(H_1)(_identity_matrix(size(H_1, 1))) + h = Vector{QQMPolyRingElem}() + success = false + while !success + empty!(h) + success = true + for i in eachindex(A) + α = A[i] + if info_level >= 1 + println("Computing leading principal minors of Hermite matrix for: ", α) + end + H_α = transpose(Q) * H[i] * Q + h_α = [det(H_α[1:j, 1:j]) for j in 1:size(H_α, 1)] + if any(is_zero, h_α) + @warn "A leading principal minor of the Hermite matrix for α = $(α) is zero." + @warn "Applying a random change of basis to all Hermite matrices and trying again." + # We could also try to generate Q deterministically, + # but a random change of basis is simpler and works well in practice. + Q = parent(H_1)(_random_matrix(size(H_1, 1), size(H_1, 1))) + success = false + break + end + append!(h, _change_ring(map(p -> numerator(p) * denominator(p), h_α), R)) + end + end + h = filter(!is_constant, h) + p = _unique_by_sign(_pick_points(points_per_components(QQMPolyRingElem[], QQMPolyRingElem[], [w; h]; nr_thrds=nr_thrds, worker_pool=worker_pool, info_level=info_level)), h) + if info_level >= 1 + println("Number of refined sample points: ", length(p)) + end + prog = Progress.ProgressBar(total=length(p); desc="Refined sign conditions", enabled=show_progress) + Progress.update!(prog, 0) + signs = Vector{Vector{Int}}() + counts = Vector{Int}() + witnesses = Vector{Vector{QQFieldElem}}() + for pᵢ in p + Σᵢ, cᵢ = sign_determination(I, g, pᵢ) + j = findfirst(==(σ), Σᵢ) + cᵢⱼ = isnothing(j) ? 0 : cᵢ[j] + if ignore_no_real_roots && cᵢⱼ == 0 + Progress.next!(prog) + continue + end + new_sign = map(p -> Int(numerator(sign(evaluate(p, pᵢ)))), h) + if !(new_sign in signs) + push!(signs, new_sign) + push!(counts, cᵢⱼ) + push!(witnesses, pᵢ) + end + Progress.next!(prog) + end + Progress.finish!(prog) + if info_level >= 1 + println("Signs at refined sample points: ", signs) + end + w = unique([q for wᵢ in w for (q, _) in factor(wᵢ)]) + # Finally, the semi-algebraic description of the real root classification is given by + # the semi-algebraic set defined by the signs of the leading principal minors of the Hermite matrices, + # restricted to the semi-algebraic set defined by the non-vanishing of the polynomials in w. + Intersect([ + BasicSemialgebraicSet(eqs=[], ineqs=w, pos=[], nonneg=[]), + SignEnum(polys=h, signs=signs, counts=counts, witnesses=witnesses) + ]) +end diff --git a/src/algorithms/param-ideal/semialgebraic-set.jl b/src/algorithms/param-ideal/semialgebraic-set.jl new file mode 100644 index 0000000..559c117 --- /dev/null +++ b/src/algorithms/param-ideal/semialgebraic-set.jl @@ -0,0 +1,351 @@ +abstract type SemialgebraicSet end + +@kwdef struct Empty <: SemialgebraicSet end + +@kwdef struct Universe <: SemialgebraicSet end + +@kwdef struct BasicSemialgebraicSet <: SemialgebraicSet + eqs::Vector{QQMPolyRingElem} + ineqs::Vector{QQMPolyRingElem} + pos::Vector{QQMPolyRingElem} + nonneg::Vector{QQMPolyRingElem} +end + +@kwdef struct SignEnum <: SemialgebraicSet + polys::Vector{QQMPolyRingElem} + signs::Vector{Vector{Int}} + counts::Vector{Int} + witnesses::Vector{Vector{QQFieldElem}} +end + +@kwdef struct SigEnum <: SemialgebraicSet + hms::Vector{MatSpaceElem{FracFieldElem{QQMPolyRingElem}}} + sigs::Vector{Vector{Int}} + counts::Vector{Int} + witnesses::Vector{Vector{QQFieldElem}} +end + +@kwdef struct MatEnum <: SemialgebraicSet + hm::MatSpaceElem{FracFieldElem{QQMPolyRingElem}} + mm::Vector{MatSpaceElem{FracFieldElem{QQMPolyRingElem}}} +end + +struct Complement <: SemialgebraicSet + S::SemialgebraicSet +end + +struct Intersect <: SemialgebraicSet + S::Vector{SemialgebraicSet} +end + +struct UnionSet <: SemialgebraicSet + S::Vector{SemialgebraicSet} +end + +# print methods + +function Base.show(io::IO, _::Empty) + print(io, "Empty()") +end + +function Base.show(io::IO, _::Universe) + print(io, "Universe()") +end + +function Base.show(io::IO, S::BasicSemialgebraicSet) + print(io, "BasicSemialgebraicSet(") + print(io, "eqs=[") + for (i, eq) in enumerate(S.eqs) + if i > 1 + print(io, ", ") + end + print(io, eq) + end + print(io, "], ineqs=[") + for (i, ineq) in enumerate(S.ineqs) + if i > 1 + print(io, ", ") + end + print(io, ineq) + end + print(io, "], pos=[") + for (i, p) in enumerate(S.pos) + if i > 1 + print(io, ", ") + end + print(io, p) + end + print(io, "], nonneg=[") + for (i, p) in enumerate(S.nonneg) + if i > 1 + print(io, ", ") + end + print(io, p) + end + print(io, "])") +end + +function Base.show(io::IO, S::SignEnum) + print(io, "SignEnum(") + print(io, "polys=[") + for (i, p) in enumerate(S.polys) + if i > 1 + print(io, ", ") + end + print(io, p) + end + print(io, "], signs=") + print(io, S.signs) + print(io, ", counts=") + print(io, S.counts) + print(io, ", witnesses=") + print(io, S.witnesses) + print(io, ")") +end + +function Base.show(io::IO, S::SigEnum) + print(io, "SigEnum(") + print(io, "hms=[") + for (i, hm) in enumerate(S.hms) + if i > 1 + print(io, ", ") + end + print(io, hm) + end + print(io, "], sigs=") + print(io, S.sigs) + print(io, ", counts=") + print(io, S.counts) + print(io, ", witnesses=") + print(io, S.witnesses) + print(io, ")") +end + +function Base.show(io::IO, S::MatEnum) + print(io, "MatEnum(") + print(io, "hm=") + print(io, S.hm) + print(io, ", mm=[") + for (i, mm) in enumerate(S.mm) + if i > 1 + print(io, ", ") + end + print(io, mm) + end + print(io, "])") +end + +function Base.show(io::IO, S::Complement) + print(io, "Complement(") + show(io, S.S) + print(io, ")") +end + +function Base.show(io::IO, S::Intersect) + print(io, "Intersect([") + for (i, Sᵢ) in enumerate(S.S) + if i > 1 + print(io, ", ") + end + show(io, Sᵢ) + end + print(io, "])") +end + +function Base.show(io::IO, S::UnionSet) + print(io, "UnionSet([") + for (i, Sᵢ) in enumerate(S.S) + if i > 1 + print(io, ", ") + end + show(io, Sᵢ) + end + print(io, "])") +end + +function Base.show(io::IO, S::SemialgebraicSet) + if S isa Empty + show(io, S::Empty) + elseif S isa Universe + show(io, S::Universe) + elseif S isa BasicSemialgebraicSet + show(io, S::BasicSemialgebraicSet) + elseif S isa SignEnum + show(io, S::SignEnum) + elseif S isa SigEnum + show(io, S::SigEnum) + elseif S isa MatEnum + show(io, S::MatEnum) + elseif S isa Complement + show(io, S::Complement) + elseif S isa Intersect + show(io, S::Intersect) + elseif S isa UnionSet + show(io, S::UnionSet) + else + error("Unknown SemialgebraicSet type") + end +end + +# equality methods + +function Base.:(==)(_::Empty, _::Empty) + return true +end + +function Base.:(==)(_::Universe, _::Universe) + return true +end + +function Base.:(==)(S1::BasicSemialgebraicSet, S2::BasicSemialgebraicSet) + return S1.eqs == S2.eqs && S1.ineqs == S2.ineqs && S1.pos == S2.pos && S1.nonneg == S2.nonneg +end + +function Base.:(==)(S1::SignEnum, S2::SignEnum) + return S1.polys == S2.polys && S1.signs == S2.signs && S1.counts == S2.counts +end + +function Base.:(==)(S1::SigEnum, S2::SigEnum) + return S1.hms == S2.hms && S1.sigs == S2.sigs && S1.counts == S2.counts +end + +function Base.:(==)(S1::MatEnum, S2::MatEnum) + return S1.hm == S2.hm && S1.mm == S2.mm +end + +function Base.:(==)(S1::Complement, S2::Complement) + return S1.S == S2.S +end + +function Base.:(==)(S1::Intersect, S2::Intersect) + return S1.S == S2.S +end + +function Base.:(==)(S1::UnionSet, S2::UnionSet) + return S1.S == S2.S +end + +function Base.:(==)(S1::SemialgebraicSet, S2::SemialgebraicSet) + if S1 isa BasicSemialgebraicSet && S2 isa BasicSemialgebraicSet + return S1 == S2 + elseif S1 isa SignEnum && S2 isa SignEnum + return S1 == S2 + elseif S1 isa SigEnum && S2 isa SigEnum + return S1 == S2 + elseif S1 isa MatEnum && S2 isa MatEnum + return S1 == S2 + elseif S1 isa Complement && S2 isa Complement + return S1 == S2 + elseif S1 isa Intersect && S2 isa Intersect + return S1 == S2 + elseif S1 isa UnionSet && S2 isa UnionSet + return S1 == S2 + else + return false + end +end + +# evaluation methods + +function evaluate(_::Empty, _::Vector{QQFieldElem})::Bool + return false +end + +function evaluate(_::Universe, _::Vector{QQFieldElem})::Bool + return true +end + +function evaluate(S::BasicSemialgebraicSet, vals::Vector{QQFieldElem})::Bool + for eq in S.eqs + if evaluate(eq, vals) != 0 + return false + end + end + for ineq in S.ineqs + if evaluate(ineq, vals) == 0 + return false + end + end + for p in S.pos + if evaluate(p, vals) <= 0 + return false + end + end + for p in S.nonneg + if evaluate(p, vals) < 0 + return false + end + end + return true +end + +function evaluate(S::SignEnum, vals::Vector{QQFieldElem})::Bool + signs = [sign(evaluate(p, vals)) for p in S.polys] + i = findfirst(==(signs), S.signs) + return isnothing(i) ? false : S.counts[i] > 0 +end + +function evaluate(S::SigEnum, vals::Vector{QQFieldElem})::Bool + hm_specs = map(hm -> map(a -> evaluate(numerator(a), vals) // evaluate(denominator(a), vals), hm), S.hms) + sig = map(hm -> signature(hm), hm_specs) + i = findfirst(==(sig), S.sigs) + return isnothing(i) ? false : S.counts[i] > 0 +end + +function evaluate(S::MatEnum, vals::Vector{QQFieldElem})::Bool + function _sum_of_signatures(_S::MatEnum, _vals::Vector{QQFieldElem})::Int + if isempty(_S.mm) + return signature(map(a -> evaluate(numerator(a), _vals) // evaluate(denominator(a), _vals), _S.hm)) + end + total = _sum_of_signatures(MatEnum(hm=_S.hm, mm=_S.mm[2:end]), _vals) + total += _sum_of_signatures(MatEnum(hm=_S.hm * _S.mm[1], mm=_S.mm[2:end]), _vals) + return total + end + return _sum_of_signatures(S, vals) > 0 +end + +function evaluate(S::Complement, vals::Vector{QQFieldElem})::Bool + return !evaluate(S.S, vals) +end + +function evaluate(S::Intersect, vals::Vector{QQFieldElem})::Bool + for Sᵢ in S.S + if !evaluate(Sᵢ, vals) + return false + end + end + return true +end + +function evaluate(S::UnionSet, vals::Vector{QQFieldElem})::Bool + for Sᵢ in S.S + if evaluate(Sᵢ, vals) + return true + end + end + return false +end + +function evaluate(S::SemialgebraicSet, vals::Vector{QQFieldElem})::Bool + if S isa Empty + return evaluate(S::Empty, vals) + elseif S isa Universe + return evaluate(S::Universe, vals) + elseif S isa BasicSemialgebraicSet + return evaluate(S::BasicSemialgebraicSet, vals) + elseif S isa SignEnum + return evaluate(S::SignEnum, vals) + elseif S isa SigEnum + return evaluate(S::SigEnum, vals) + elseif S isa MatEnum + return evaluate(S::MatEnum, vals) + elseif S isa Complement + return evaluate(S::Complement, vals) + elseif S isa Intersect + return evaluate(S::Intersect, vals) + elseif S isa UnionSet + return evaluate(S::UnionSet, vals) + else + error("Unknown SemialgebraicSet type") + end +end diff --git a/src/algorithms/param-ideal/sign-determination.jl b/src/algorithms/param-ideal/sign-determination.jl new file mode 100644 index 0000000..de1116e --- /dev/null +++ b/src/algorithms/param-ideal/sign-determination.jl @@ -0,0 +1,243 @@ +# Notation 2.43 (Sign variations) +@doc Markdown.doc""" + sign_variations(s::Vector{T}) -> Int + +Computes the number of sign variations in the sequence `s`, ignoring zeros. +""" +function sign_variations(s::Vector{T}) where {T} + s′ = filter(x -> !is_zero(x), s) + if isempty(s′) + return 0 + end + count = 0 + for i in 1:(length(s′)-1) + if s′[i] * s′[i+1] < 0 + count += 1 + end + end + return count +end + +@doc Markdown.doc""" + alternate_signs(s::Vector{T}) -> Vector{T} + +Returns a new vector where the signs of the elements in `s` are alternated, starting with the first element unchanged. +""" +function alternate_signs(s::Vector{T})::Vector{T} where {T} + s′ = copy(s) + for i in (length(s)%2+1):2:length(s) + s′[i] = -s′[i] + end + return s′ +end + +@doc Markdown.doc""" + signature(M::QQMatrix) -> Int + +Computes the signature of the real symmetric matrix `M`, i.e., the number of positive eigenvalues minus the number of negative eigenvalues. +""" +function signature(M::QQMatrix)::Int + @assert is_symmetric(M) + l = collect(coefficients(charpoly(M))) + sign_variations(l) - sign_variations(alternate_signs(l)) +end + +@doc Markdown.doc""" + matrix_of_signs(A::Vector{Vector{Int}}, Σ::Vector{Vector{Int}}) -> Array{Int,2} + +Computes the matrix of signs of `A` on `Σ`. +""" +function matrix_of_signs(A::Vector{Vector{Int}}, Σ::Vector{Vector{Int}})::Array{Int,2} + m = length(A) + n = length(Σ) + M = ones(Int, m, n) + for i in 1:m + for j in 1:n + for (αᵢ, σⱼ) in zip(A[i], Σ[j]) + if αᵢ == 0 + continue + elseif σⱼ == 0 + M[i, j] = 0 + elseif αᵢ == 1 + M[i, j] *= σⱼ + elseif αᵢ == -1 + M[i, j] /= σⱼ + end + end + end + end + M +end + +@doc Markdown.doc""" + sort_indices(A::Vector{Vector{Int}}) -> Vector{Vector{Int}} + +Sorts the list of multi-indices `A`. +""" +function sort_indices(A::Vector{Vector{Int}})::Vector{Vector{Int}} + sort(A, by=α -> sum(α[i] * 3^(i - 1) for i in eachindex(α))) +end + +@doc Markdown.doc""" + sort_signs(Σ::Vector{Vector{Int}}) -> Vector{Vector{Int}} + +Sorts the list of sign vectors `Σ`. +""" +function sort_signs(Σ::Vector{Vector{Int}})::Vector{Vector{Int}} + sort(Σ, by=σ -> map(σᵢ -> σᵢ == -1 ? 2 : σᵢ, σ)) +end + +@doc Markdown.doc""" + extend_signs(Σ::Vector{Vector{Int}}, T::Vector{Int}) -> Vector{Vector{Int}} + +Extends the list of sign vectors `Σ` by appending each sign in `T` to each vector in `Σ`. +""" +function extend_signs(Σ::Vector{Vector{Int}}, T::Vector{Int})::Vector{Vector{Int}} + Σ′ = Vector{Vector{Int}}() + for σ in Σ + for τ in T + push!(Σ′, [σ; τ]) + end + end + Σ′ +end + +@doc Markdown.doc""" + extend_indices(A::Vector{Vector{Int}}, B::Vector{Int}) -> Vector{Vector{Int}} + +Extends the list of multi-indices `A` by appending each element in `B` to each vector in `A`. +""" +function extend_indices(A::Vector{Vector{Int}}, B::Vector{Int})::Vector{Vector{Int}} + A′ = Vector{Vector{Int}}() + for β in B + for α in A + push!(A′, [α; β]) + end + end + A′ +end + +@doc Markdown.doc""" + modified_tensor_product(M::Matrix{Int}, M′::Matrix{Int}) -> Matrix{Int} + +Computes the modified tensor product of the matrices `M` and `M′`. +""" +function modified_tensor_product(M::Matrix{Int}, M′::Matrix{Int}) + n, m = size(M) + n′, m′ = size(M′) + M′′ = zeros(Int, n * n′, m * m′) + for i in 1:n + for j in 1:m + for i′ in 1:n′ + for j′ in 1:m′ + M′′[n*(i′-1)+i, m′*(j-1)+j′] = M[i, j] * M′[i′, j′] + end + end + end + end + M′′ +end + +@doc Markdown.doc""" + adapted_family(Σ::Vector{Vector{Int}}) -> Vector{Vector{Int}} + +Computes a family of multi-indices adapted to sign determination on the list of sign vectors `Σ`. +""" +function adapted_family(Σ::Vector{Vector{Int}})::Vector{Vector{Int}} + if isempty(Σ) + return [] + end + s = length(Σ[1]) + if s == 0 + return [[]] + end + Ξ⁽⁰⁾ = Vector{Vector{Int}}([]) + Ξ⁽¹⁾ = Vector{Vector{Int}}([]) + Ξ⁽²⁾ = Vector{Vector{Int}}([]) + for σ in Σ + τ = σ[1:(s-1)] + if τ in Ξ⁽⁰⁾ && τ in Ξ⁽¹⁾ + Ξ⁽²⁾ = push!(Ξ⁽²⁾, τ) + elseif τ in Ξ⁽⁰⁾ + Ξ⁽¹⁾ = push!(Ξ⁽¹⁾, τ) + else + Ξ⁽⁰⁾ = push!(Ξ⁽⁰⁾, τ) + end + end + [ + extend_indices(adapted_family(Ξ⁽⁰⁾), [0]); + extend_indices(adapted_family(Ξ⁽¹⁾), [1]); + extend_indices(adapted_family(Ξ⁽²⁾), [2]) + ] +end + +@doc Markdown.doc""" + tarski_query(g::Vector{MPoly{K}}, α::Vector{Int}, I::ParametricIdeal{K}, vals::Vector{QQFieldElem}=QQFieldElem[]; nr_thrds::Int=1, worker_pool::AbstractWorkerPool=default_worker_pool(), show_progress::Bool=false) -> Int + +Computes the Tarski query of `g^α` with respect to the parametric ideal `I` at the parameter values `vals`. +""" +function tarski_query( + g::Vector{MPoly{K}}, + α::Vector{Int}, + I::ParametricIdeal{K}, + vals::Vector{QQFieldElem}=QQFieldElem[]; + nr_thrds::Int=1, + worker_pool::AbstractWorkerPool=default_worker_pool(), + show_progress::Bool=false +)::Int where {K<:FracFieldElem} + @assert length(vals) == I.num_params + return signature(hermite_matrix(g, α, I, vals; nr_thrds=nr_thrds, worker_pool=worker_pool, show_progress=show_progress)) +end + +@doc Markdown.doc""" + tarski_query(g::Vector{MPoly{K}}, A::Vector{Vector{Int}}, I::ParametricIdeal{K}, vals::Vector{QQFieldElem}=QQFieldElem[]; nr_thrds::Int=1, worker_pool::AbstractWorkerPool=default_worker_pool(), show_progress::Bool=false) -> Vector{Vector{Int}} + +Computes the Tarski query of `g^α` with respect to the parametric ideal `I` at the parameter values `vals` for each multi-index `α` in `A`. +""" +function tarski_query( + g::Vector{MPoly{K}}, + A::Vector{Vector{Int}}, + I::ParametricIdeal{K}, + vals::Vector{QQFieldElem}=QQFieldElem[]; + nr_thrds::Int=1, + worker_pool::AbstractWorkerPool=default_worker_pool(), + show_progress::Bool=false +)::Vector{Vector{Int}} where {K<:FracFieldElem} + [[tarski_query(g, α, I, vals; nr_thrds=nr_thrds, worker_pool=worker_pool, show_progress=show_progress)] for α in A] +end + +@doc Markdown.doc""" + sign_determination(I::ParametricIdeal{K}, g::Vector{MPoly{K}}, vals::Vector{QQFieldElem}=QQFieldElem[]; nr_thrds::Int=1, worker_pool::AbstractWorkerPool=default_worker_pool(), show_progress::Bool=false) -> Tuple{Vector{Vector{Int}}, Vector{Int}} + +Computes the realizable sign conditions of the polynomials `g` with respect to the parametric ideal `I` at the parameter values `vals`. The result is a tuple containing a list of sign vectors and a list of number of real roots corresponding to each sign vector. +""" +function sign_determination( + I::ParametricIdeal{K}, + g::Vector{MPoly{K}}, + vals::Vector{QQFieldElem}=QQFieldElem[]; + nr_thrds::Int=1, + worker_pool::AbstractWorkerPool=default_worker_pool(), + show_progress::Bool=false +)::Tuple{Vector{Vector{Int}},Vector{Int}} where {K<:FracFieldElem} + s = length(g) + r = tarski_query(g, fill(0, s), I, vals; nr_thrds=nr_thrds, worker_pool=worker_pool, show_progress=show_progress) + if r == 0 + return (Vector{Vector{Int}}([]), Vector{Int}()) + end + Σ = [Vector{Int}()] + c = [r] + A = [Vector{Int}()] + for i in 1:length(g) + Σ′ = extend_signs(Σ, [0, 1, -1]) + A′ = extend_indices(A, [0, 1, 2]) + M = matrix(QQ, matrix_of_signs(A′, Σ′)) + T = matrix(QQ, tarski_query(g[1:i], A′, I, vals; nr_thrds=nr_thrds, worker_pool=worker_pool, show_progress=show_progress)) + c′ = inv(M) * T + @assert all(denominator(c′[j, 1]) == 1 for j in 1:length(Σ′)) + Σ = [Σ′[j] for j in 1:length(Σ′) if !is_zero(c′[j, 1])] + c = [Int(numerator(c′[j, 1])) for j in 1:length(Σ′) if !is_zero(c′[j, 1])] + @assert sum(c) == r + A = adapted_family(Σ) + end + (Σ, c) +end diff --git a/src/exports.jl b/src/exports.jl index 726feae..2d3fb04 100644 --- a/src/exports.jl +++ b/src/exports.jl @@ -4,6 +4,7 @@ export polynomial_ring, MPolyRing, GFElem, MPolyRingElem, QQMPolyRingElem, base_ring, coefficient_ring, evaluate, prime_field, sig_groebner_basis, cyclic, leading_coefficient, points_per_components, + multiplication_matrix, hermite_matrix, real_root_classification, equidimensional_decomposition, homogenize, dimension, FqMPolyRingElem, hilbert_series, hilbert_dimension, hilbert_degree, hilbert_polynomial, rational_curve_parametrization diff --git a/src/imports.jl b/src/imports.jl index 07f2043..76dabaf 100644 --- a/src/imports.jl +++ b/src/imports.jl @@ -27,6 +27,7 @@ import Nemo: fraction_field, GF, height, + inv, is_prime, is_probable_prime, is_square, diff --git a/test/algorithms/param-ideal/hermite-matrices.jl b/test/algorithms/param-ideal/hermite-matrices.jl index 6598733..ec5f3c6 100644 --- a/test/algorithms/param-ideal/hermite-matrices.jl +++ b/test/algorithms/param-ideal/hermite-matrices.jl @@ -5,7 +5,7 @@ import Nemo: zero_matrix f = [R(1)] g = AlgebraicSolving.groebner_basis(AlgebraicSolving.Ideal(f)) b = AlgebraicSolving._monomial_basis(g) - H₁ = AlgebraicSolving.hermite_matrix(g, b, one(R)) + H₁ = hermite_matrix(g, b, one(R)) @test H₁ == zero_matrix(QQ, 0, 0) end @@ -14,10 +14,10 @@ end f = [x^2 - 2 * x + 1] g = AlgebraicSolving.groebner_basis(AlgebraicSolving.Ideal(f)) b = AlgebraicSolving._monomial_basis(g) - H₁ = AlgebraicSolving.hermite_matrix(g, b, one(R)) + H₁ = hermite_matrix(g, b, one(R)) @test H₁ == matrix(QQ, [2 2; 2 2]) - Hₓ = AlgebraicSolving.hermite_matrix(g, b, x) - Mₓ = AlgebraicSolving.multiplication_matrix(g, b, x) + Hₓ = hermite_matrix(g, b, x) + Mₓ = multiplication_matrix(g, b, x) @test Hₓ == H₁ * Mₓ end @@ -26,7 +26,7 @@ end R′, (x,) = polynomial_ring(fraction_field(R), [:x], internal_ordering=:degrevlex) f = ParametricIdeal([a * x^2 + b * x + c]) g = one(R′) - H = AlgebraicSolving.hermite_matrix(f, g) + H = hermite_matrix(f, g) @test H == matrix(fraction_field(R), [2 -b//a; -b//a (b^2-2*a*c)//a^2]) end @@ -35,6 +35,6 @@ end R′, (x,) = polynomial_ring(fraction_field(R), [:x], internal_ordering=:degrevlex) f = ParametricIdeal([x^2 - 2 * x + 1]) g = one(R′) - H₁ = AlgebraicSolving.hermite_matrix(f, g) + H₁ = hermite_matrix(f, g) @test H₁ == matrix(QQ, [2 2; 2 2]) end diff --git a/test/algorithms/param-ideal/multiplication-matrices.jl b/test/algorithms/param-ideal/multiplication-matrices.jl index 1c2bb10..64389b1 100644 --- a/test/algorithms/param-ideal/multiplication-matrices.jl +++ b/test/algorithms/param-ideal/multiplication-matrices.jl @@ -5,14 +5,14 @@ import Nemo: identity_matrix, matrix f = [x^2 - 2 * x + 1] g = AlgebraicSolving.groebner_basis(AlgebraicSolving.Ideal(f)) b = AlgebraicSolving._monomial_basis(g) - M₁ = AlgebraicSolving.multiplication_matrix(g, b, one(R)) + M₁ = multiplication_matrix(g, b, one(R)) @test M₁ == identity_matrix(QQ, 2) - Mₓ = AlgebraicSolving.multiplication_matrix(g, b, x) + Mₓ = multiplication_matrix(g, b, x) @test Mₓ == matrix(QQ, [0 -1; 1 2]) end @testset "Algorithms -> Parametric multiplication matrix" begin - M = AlgebraicSolving.multiplication_matrix + M = multiplication_matrix R, (a, b, c) = polynomial_ring(QQ, [:a, :b, :c]) R′, (x,) = polynomial_ring(fraction_field(R), [:x], internal_ordering=:degrevlex) f = ParametricIdeal([a * x^2 + b * x + c]) diff --git a/test/algorithms/param-ideal/real-root-classification.jl b/test/algorithms/param-ideal/real-root-classification.jl new file mode 100644 index 0000000..1855d5b --- /dev/null +++ b/test/algorithms/param-ideal/real-root-classification.jl @@ -0,0 +1,28 @@ +@testset "Algorithms -> Real root classification" begin + R, (a, b, c) = polynomial_ring(QQ, [:a, :b, :c]) + R′, (x,) = polynomial_ring(fraction_field(R), [:x], internal_ordering=:degrevlex) + f = ParametricIdeal([a * x^2 + b * x + c]) + # Mathematically, the set of (a,b,c) in R^3 such that + # a*x^2+b*x+c has a real root satisfying a*x+c > 0 + # is given by (a != 0) and (b^2-4ac != 0) and (2c-b > 0 or c(b-a-c) > 0) + rrc = real_root_classification(f, [a * x + c]) + ineqs = rrc.S[1].ineqs + polys = unique([q for p in rrc.S[2].polys for (q, _) in factor(p)]) + # We verify that the polynomials in the semi-algebraic description + # appear in the expected set of polynomials + @test a in ineqs + @test 4*a*c-b^2 in ineqs + @test b-2*c in polys + @test a-b+c in polys +end + +@testset "Algorithms -> Real root classification 2" begin + R, (a,) = polynomial_ring(QQ, [:a]) + R′, (x,) = polynomial_ring(fraction_field(R), [:x], internal_ordering=:degrevlex) + f = ParametricIdeal([x^2 - x, x^2 - (a + 1) * x + a]) + rrc = real_root_classification(f, [1-2*x]) + # When a = 1, the ideal is generated by x^2-x instead of x-1, so x = 0 is such that 1-2*x > 0. + # In other words, evaluate(rrc, [QQ(1)]) should be true, which is not the case. + # This is why the real root classification is only guaranteed to be correct generically. + @test !evaluate(rrc, [QQ(1)]) +end diff --git a/test/algorithms/param-ideal/semialgebraic-set.jl b/test/algorithms/param-ideal/semialgebraic-set.jl new file mode 100644 index 0000000..b4b5b9d --- /dev/null +++ b/test/algorithms/param-ideal/semialgebraic-set.jl @@ -0,0 +1,9 @@ +@testset "Algorithms -> Semialgebraic sets" begin + R, (x, y) = polynomial_ring(QQ, ["x", "y"]) + S = AlgebraicSolving.Intersect([ + AlgebraicSolving.BasicSemialgebraicSet(eqs=[x^2 + y^2 - 1], ineqs=[x - y], pos=[], nonneg=[]), + AlgebraicSolving.Complement(AlgebraicSolving.BasicSemialgebraicSet(eqs=[], ineqs=[x + y - 1], pos=[], nonneg=[])) + ]) + v = evaluate(S, [QQ(1 // 2), QQ(1 // 2)]) + @test !v +end diff --git a/test/algorithms/param-ideal/sign-determination.jl b/test/algorithms/param-ideal/sign-determination.jl new file mode 100644 index 0000000..56bac7e --- /dev/null +++ b/test/algorithms/param-ideal/sign-determination.jl @@ -0,0 +1,115 @@ +Var = AlgebraicSolving.sign_variations +Sign = AlgebraicSolving.signature +Mat = AlgebraicSolving.matrix_of_signs +Ada = AlgebraicSolving.adapted_family +∧ = AlgebraicSolving.extend_signs +× = AlgebraicSolving.extend_indices +⊗′ = AlgebraicSolving.modified_tensor_product + +@testset "Algorithms -> Sign variations" begin + @test Var([]) == 0 + @test Var([1, -1, 2, 0, 0, 3, 4, -5, -2, 0, 3]) == 4 +end + +@testset "Algorithms -> Signature of quadratic forms" begin + @test Sign(matrix(QQ, [1;;])) == 1 + @test Sign(matrix(QQ, [1 2; 2 3])) == 0 + @test Sign(matrix(QQ, [2 0; 0 0])) == 1 +end + +@testset "Algorithms -> Matrix of signs" begin + Σ = [Vector{Int}()] ∧ [0, 1, -1] + A = [Vector{Int}()] × [0, 1, 2] + @test Σ == [[0], [1], [-1]] + @test A == [[0], [1], [2]] + M = M′ = Mat(A, Σ) + @test M == [ + 1 1 1 + 0 1 -1 + 0 1 1 + ] + @test M⊗′M′ == [ + 1 1 1 1 1 1 1 1 1; + 0 0 0 1 1 1 -1 -1 -1; + 0 0 0 1 1 1 1 1 1; + 0 1 -1 0 1 -1 0 1 -1; + 0 0 0 0 1 -1 0 -1 1; + 0 0 0 0 1 -1 0 1 -1; + 0 1 1 0 1 1 0 1 1; + 0 0 0 0 1 1 0 -1 -1; + 0 0 0 0 1 1 0 1 1 + ] + Σ = Σ ∧ [0, 1, -1] + A = A × [0, 1, 2] + @test Mat(A, Σ) == M⊗′M′ +end + +# Example 10.74 +@testset "Algorithms -> Adapted family 1" begin + Σ = [[0], [1], [-1]] + @test Ada(Σ) == [[0], [1], [2]] + @test Mat(Ada(Σ), Σ) == [1 1 1; 0 1 -1; 0 1 1] + Σ = [[1], [-1]] + @test Ada(Σ) == [[0], [1]] + @test Mat(Ada(Σ), Σ) == [1 1; 1 -1] + Σ = [[0], [1]] + @test Ada(Σ) == [[0], [1]] + @test Mat(Ada(Σ), Σ) == [1 1; 0 1] + Σ = [[0], [-1]] + @test Ada(Σ) == [[0], [1]] + @test Mat(Ada(Σ), Σ) == [1 1; 0 -1] + Σ = [[0]] + @test Ada(Σ) == [[0]] + @test Mat(Ada(Σ), Σ) == [1;;] + Σ = [[1]] + @test Ada(Σ) == [[0]] + @test Mat(Ada(Σ), Σ) == [1;;] + Σ = [[-1]] + @test Ada(Σ) == [[0]] + @test Mat(Ada(Σ), Σ) == [1;;] +end + +# Example 10.81 +@testset "Algorithms -> Adapted family 2" begin + Σ = [ + [0, 1, 0, 0], [0, 1, 0, 1], [0, 1, 0, -1], [0, 1, 1, -1], + [1, -1, 0, 0], [1, -1, 0, 1], [1, -1, 0, -1], [-1, 0, 0, 0], + [-1, 0, 0, -1], [-1, 0, 1, 1], [-1, 0, 1, -1] + ] + @test Ada(Σ) == [ + [0, 0, 0, 0], [1, 0, 0, 0], [2, 0, 0, 0], [0, 0, 1, 0], + [1, 0, 1, 0], [0, 0, 0, 1], [1, 0, 0, 1], [2, 0, 0, 1], + [0, 0, 1, 1], [0, 0, 0, 2], [1, 0, 0, 2] + ] +end + +@testset "Algorithms -> Adapted family 3" begin + s = 10 + Σ = [Vector{Int}()] + for i in 1:s + Σ = Σ ∧ [1, -1] + end + A = Ada(Σ) + M = matrix(QQ, Mat(A, Σ)) + @test size(M) == (2^s, 2^s) + M = inv(M) + @test M[end, end] == QQ(1 // 2^s) +end + +@testset "Algorithms -> Sign determination 1" begin + R, () = polynomial_ring(QQ, Symbol[]) + R′, (x1, x2) = polynomial_ring(fraction_field(R), [:x1, :x2]) + f = ParametricIdeal([x1^2 + x2^2 - 1, x1^2 - x2]) + g₁ = x1 + x2 + g₂ = x1 - x2 + S = AlgebraicSolving.sign_determination(f, [g₁, g₂]) + @test S == ([[1, 1], [-1, -1]], [1, 1]) +end + +@testset "Algorithms -> Sign determination 2" begin + R, () = polynomial_ring(QQ, Symbol[]) + R′, (x1, x2, x3) = polynomial_ring(fraction_field(R), [:x1, :x2, :x3]) + f = ParametricIdeal([x1^2 - 1, x2^2 - 1, x3^3 - x3]) + S = AlgebraicSolving.sign_determination(f, [x1, x2, x3]) + @test length(S[1]) == 12 +end diff --git a/test/runtests.jl b/test/runtests.jl index eb57df5..53d8844 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -16,6 +16,9 @@ include("algorithms/param-ideal/groebner-bases.jl") include("algorithms/param-ideal/multiplication-matrices.jl") include("algorithms/param-ideal/hermite-matrices.jl") include("algorithms/param-ideal/parametrizations.jl") +include("algorithms/param-ideal/sign-determination.jl") +include("algorithms/param-ideal/semialgebraic-set.jl") +include("algorithms/param-ideal/real-root-classification.jl") include("examples/katsura.jl") include("interp/thiele.jl") include("interp/newton.jl")