diff --git a/.JuliaFormatter.toml b/.JuliaFormatter.toml new file mode 100644 index 0000000..f196f10 --- /dev/null +++ b/.JuliaFormatter.toml @@ -0,0 +1,9 @@ +# Configuration file for JuliaFormatter.jl +# For more information, see: https://domluna.github.io/JuliaFormatter.jl/stable/config/ + +always_for_in = true +always_use_return = true +margin = 80 +remove_extra_newlines = true +separate_kwargs_with_semicolon = true +short_to_long_function_def = true diff --git a/.github/workflows/format_check.yml b/.github/workflows/format_check.yml new file mode 100644 index 0000000..232768a --- /dev/null +++ b/.github/workflows/format_check.yml @@ -0,0 +1,32 @@ +ame: format-check +on: + push: + branches: + - master + - release-* + pull_request: + types: [opened, synchronize, reopened] +jobs: + build: + runs-on: ubuntu-latest + steps: + - uses: julia-actions/setup-julia@latest + with: + version: '1' + - uses: actions/checkout@v3 + - name: Format check + shell: julia --color=yes {0} + run: | + using Pkg + # If you update the version, also update the style guide docs. + Pkg.add(PackageSpec(name="JuliaFormatter", version="1")) + using JuliaFormatter + format("src", verbose=true) + format("test", verbose=true) + out = String(read(Cmd(`git diff`))) + if isempty(out) + exit(0) + end + @error "Some files have not been formatted !!!" + write(stdout, out) + exit(1) diff --git a/Project.toml b/Project.toml index 72f5b3c..a73359c 100644 --- a/Project.toml +++ b/Project.toml @@ -13,11 +13,3 @@ Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" MultivariatePolynomials = "0.4.2" MutableArithmetics = "0.3.1, 1" julia = "1.6" - -[extras] -DynamicPolynomials = "7c1d4256-1411-5781-91ec-d7bc3513ac07" -Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" -TypedPolynomials = "afbbf031-7a57-5f58-a1b9-b774a0fad08d" - -[targets] -test = ["Test", "DynamicPolynomials", "TypedPolynomials"] diff --git a/src/SemialgebraicSets.jl b/src/SemialgebraicSets.jl index 7cff6cb..968e8ab 100644 --- a/src/SemialgebraicSets.jl +++ b/src/SemialgebraicSets.jl @@ -10,8 +10,10 @@ const MP = MultivariatePolynomials const APL = AbstractPolynomialLike -export AbstractSemialgebraicSet, AbstractBasicSemialgebraicSet, AbstractAlgebraicSet -export FullSpace, AlgebraicSet, BasicSemialgebraicSet, addequality!, addinequality! +export AbstractSemialgebraicSet, + AbstractBasicSemialgebraicSet, AbstractAlgebraicSet +export FullSpace, + AlgebraicSet, BasicSemialgebraicSet, addequality!, addinequality! # Semialgebraic set described by polynomials with coefficients in T abstract type AbstractSemialgebraicSet end @@ -19,12 +21,13 @@ abstract type AbstractSemialgebraicSet end abstract type AbstractBasicSemialgebraicSet <: AbstractSemialgebraicSet end abstract type AbstractAlgebraicSet <: AbstractBasicSemialgebraicSet end -addinequality!(S::AbstractAlgebraicSet, p) = throw(ArgumentError("Cannot add inequality to an algebraic set")) - -struct FullSpace <: AbstractAlgebraicSet +function addinequality!(S::AbstractAlgebraicSet, p) + throw(ArgumentError("Cannot add inequality to an algebraic set")) end + +struct FullSpace <: AbstractAlgebraicSet end function Base.show(io::IO, ::FullSpace) - print(io, "R^n") + return print(io, "R^n") end nequalities(::FullSpace) = 0 equalities(::FullSpace) = [] @@ -39,7 +42,14 @@ Base.intersect(S::FullSpace, T::AbstractAlgebraicSet) = T Base.intersect(S::FullSpace, T::FullSpace) = S # If `intersect(S, T)` is not implemented, this method will `StackOverflow`. -Base.intersect(S::AbstractSemialgebraicSet, T::AbstractSemialgebraicSet, args...; kws...) = intersect(intersect(S, T), args...; kws...) +function Base.intersect( + S::AbstractSemialgebraicSet, + T::AbstractSemialgebraicSet, + args...; + kws..., +) + return intersect(intersect(S, T), args...; kws...) +end # The keywords are only used when transforming `Element` # into `BasicSemialgebraicSet`. Base.intersect(set::AbstractSemialgebraicSet; kws...) = set diff --git a/src/basic.jl b/src/basic.jl index 6b66829..d34bc5f 100644 --- a/src/basic.jl +++ b/src/basic.jl @@ -1,33 +1,61 @@ export inequalities, basicsemialgebraicset -struct BasicSemialgebraicSet{T, PT<:APL{T}, AT <: AbstractAlgebraicSet} <: AbstractBasicSemialgebraicSet +struct BasicSemialgebraicSet{T,PT<:APL{T},AT<:AbstractAlgebraicSet} <: + AbstractBasicSemialgebraicSet V::AT p::Vector{PT} end -function BasicSemialgebraicSet{T, PT}() where {T, PT<:APL{T}} - BasicSemialgebraicSet(AlgebraicSet{T, PT}(), PT[]) +function BasicSemialgebraicSet{T,PT}() where {T,PT<:APL{T}} + return BasicSemialgebraicSet(AlgebraicSet{T,PT}(), PT[]) end -function BasicSemialgebraicSet(V::AlgebraicSet{T, PT, A, S}, p::Vector{PT}) where {T, PT<:APL{T}, A, S<:AbstractAlgebraicSolver} - BasicSemialgebraicSet{T, PT, typeof(V)}(V, p) +function BasicSemialgebraicSet( + V::AlgebraicSet{T,PT,A,S}, + p::Vector{PT}, +) where {T,PT<:APL{T},A,S<:AbstractAlgebraicSolver} + return BasicSemialgebraicSet{T,PT,typeof(V)}(V, p) end -function BasicSemialgebraicSet(V::AlgebraicSet{T, PT, A, SO, U}, p::Vector{PS}) where {T, PT<:APL{T}, S, PS<:APL{S}, A, SO<:AbstractAlgebraicSolver, U} +function BasicSemialgebraicSet( + V::AlgebraicSet{T,PT,A,SO,U}, + p::Vector{PS}, +) where {T,PT<:APL{T},S,PS<:APL{S},A,SO<:AbstractAlgebraicSolver,U} ST = promote_type(T, S) PST = promote_type(PT, PS) - BasicSemialgebraicSet(convert(AlgebraicSet{ST, PST, A, SO, U}, V), Vector{PST}(p)) + return BasicSemialgebraicSet( + convert(AlgebraicSet{ST,PST,A,SO,U}, V), + Vector{PST}(p), + ) end #BasicSemialgebraicSet{T, PT<:APL{T}}(V::AlgebraicSet{T, PT}, p::Vector{PT}) = BasicSemialgebraicSet{T, PT}(V, p) function basicsemialgebraicset(V, p) - BasicSemialgebraicSet(V, p) + return BasicSemialgebraicSet(V, p) end -MP.similar_type(::Type{BasicSemialgebraicSet{S, PS, AT}}, T::Type) where {S, PS, AT} = BasicSemialgebraicSet{T, MP.similar_type(PS, T), MP.similar_type(AT, T)} +function MP.similar_type( + ::Type{BasicSemialgebraicSet{S,PS,AT}}, + T::Type, +) where {S,PS,AT} + return BasicSemialgebraicSet{ + T, + MP.similar_type(PS, T), + MP.similar_type(AT, T), + } +end -function Base.convert(::Type{BasicSemialgebraicSet{T, PT, AT}}, set::BasicSemialgebraicSet) where {T, PT, AT} - return BasicSemialgebraicSet{T, PT, AT}(set.V, set.p) +function Base.convert( + ::Type{BasicSemialgebraicSet{T,PT,AT}}, + set::BasicSemialgebraicSet, +) where {T,PT,AT} + return BasicSemialgebraicSet{T,PT,AT}(set.V, set.p) end -MP.variables(S::BasicSemialgebraicSet{T, PT, FullSpace}) where {T, PT<:APL{T}} = MP.variables(S.p) -MP.variables(S::BasicSemialgebraicSet) = sort(union(MP.variables(S.V), MP.variables(S.p)), rev=true) +function MP.variables( + S::BasicSemialgebraicSet{T,PT,FullSpace}, +) where {T,PT<:APL{T}} + return MP.variables(S.p) +end +function MP.variables(S::BasicSemialgebraicSet) + return sort(union(MP.variables(S.V), MP.variables(S.p)); rev = true) +end nequalities(S::BasicSemialgebraicSet) = nequalities(S.V) equalities(S::BasicSemialgebraicSet) = equalities(S.V) addequality!(S::BasicSemialgebraicSet, p) = addequality!(S.V, p) @@ -35,21 +63,32 @@ ninequalities(S::BasicSemialgebraicSet) = length(S.p) inequalities(S::BasicSemialgebraicSet) = S.p addinequality!(S::BasicSemialgebraicSet, p) = push!(S.p, p) -Base.intersect(S::BasicSemialgebraicSet, T::BasicSemialgebraicSet) = BasicSemialgebraicSet(S.V ∩ T.V, [S.p; T.p]) -Base.intersect(S::BasicSemialgebraicSet, T::AbstractAlgebraicSet) = BasicSemialgebraicSet(S.V ∩ T, copy(S.p)) +function Base.intersect(S::BasicSemialgebraicSet, T::BasicSemialgebraicSet) + return BasicSemialgebraicSet(S.V ∩ T.V, [S.p; T.p]) +end +function Base.intersect(S::BasicSemialgebraicSet, T::AbstractAlgebraicSet) + return BasicSemialgebraicSet(S.V ∩ T, copy(S.p)) +end Base.intersect(S::BasicSemialgebraicSet, ::FullSpace) = S -Base.intersect(T::AbstractAlgebraicSet, S::BasicSemialgebraicSet) = intersect(S, T) +function Base.intersect(T::AbstractAlgebraicSet, S::BasicSemialgebraicSet) + return intersect(S, T) +end function Base.show(io::IO, V::BasicSemialgebraicSet) - print(io, "{ (", join(variables(V), ", "), ") | ", - join(string.(equalities(V)) .* " = 0", ", ")) + print( + io, + "{ (", + join(variables(V), ", "), + ") | ", + join(string.(equalities(V)) .* " = 0", ", "), + ) if nequalities(V) > 0 print(io, ", ") end - print(io, join(string.(inequalities(V)) .* " ≥ 0", ", "), " }") + return print(io, join(string.(inequalities(V)) .* " ≥ 0", ", "), " }") end function Base.show(io::IO, mime::MIME"text/plain", V::BasicSemialgebraicSet) print(io, "Basic semialgebraic Set defined by ") _show_els(io, "equalit", nequalities(V), equalities(V), "=") - _show_els(io, "inequalit", ninequalities(V), inequalities(V), "≥") + return _show_els(io, "inequalit", ninequalities(V), inequalities(V), "≥") end diff --git a/src/fix.jl b/src/fix.jl index e49cad2..61b8a46 100644 --- a/src/fix.jl +++ b/src/fix.jl @@ -1,24 +1,34 @@ export FixedVariablesSet -struct FixedVariablesIdeal{V<:AbstractVariable, T<:Number, MT<:AbstractMonomialLike} <: AbstractPolynomialIdeal - substitutions::Union{Nothing, Dict{V, T}} -end -function Base.convert(::Type{FixedVariablesIdeal{V, T, MT}}, ideal::FixedVariablesIdeal{V, T, MT}) where {V<:AbstractVariable, T<:Number, MT<:AbstractMonomialLike} +struct FixedVariablesIdeal{ + V<:AbstractVariable, + T<:Number, + MT<:AbstractMonomialLike, +} <: AbstractPolynomialIdeal + substitutions::Union{Nothing,Dict{V,T}} +end +function Base.convert( + ::Type{FixedVariablesIdeal{V,T,MT}}, + ideal::FixedVariablesIdeal{V,T,MT}, +) where {V<:AbstractVariable,T<:Number,MT<:AbstractMonomialLike} return ideal end -function Base.convert(::Type{FixedVariablesIdeal{V, T, MT}}, ideal::FixedVariablesIdeal{V, S, MT}) where {V, S, T, MT} +function Base.convert( + ::Type{FixedVariablesIdeal{V,T,MT}}, + ideal::FixedVariablesIdeal{V,S,MT}, +) where {V,S,T,MT} subs = ideal.substitutions if subs !== nothing - subs = convert(Dict{V, T}, subs) + subs = convert(Dict{V,T}, subs) end - return FixedVariablesIdeal{V, T, MT}(subs) + return FixedVariablesIdeal{V,T,MT}(subs) end -function MP.variables(ideal::FixedVariablesIdeal{V}) where V +function MP.variables(ideal::FixedVariablesIdeal{V}) where {V} if generate_nonzero_constant(ideal) return V[] else - return sort(collect(keys(ideal.substitutions)), rev=true) + return sort(collect(keys(ideal.substitutions)); rev = true) end end # In that case, the ideal can generate any polynomial. @@ -34,15 +44,26 @@ function Base.rem(p::AbstractPolynomialLike, I::FixedVariablesIdeal) end end -struct FixedVariablesSet{V, T, MT} <: AbstractAlgebraicSet - ideal::FixedVariablesIdeal{V, T, MT} -end -MP.similar_type(::Type{FixedVariablesSet{V, S, MT}}, T::Type) where {V, S, MT} = FixedVariablesSet{V, T, MT} -function Base.convert(::Type{FixedVariablesSet{V, T, MT}}, set::FixedVariablesSet{V, T, MT}) where {V, T, MT} +struct FixedVariablesSet{V,T,MT} <: AbstractAlgebraicSet + ideal::FixedVariablesIdeal{V,T,MT} +end +function MP.similar_type( + ::Type{FixedVariablesSet{V,S,MT}}, + T::Type, +) where {V,S,MT} + return FixedVariablesSet{V,T,MT} +end +function Base.convert( + ::Type{FixedVariablesSet{V,T,MT}}, + set::FixedVariablesSet{V,T,MT}, +) where {V,T,MT} return set end -function Base.convert(::Type{FixedVariablesSet{V, T, MT}}, set::FixedVariablesSet{V, S, MT}) where {V, S, T, MT} - return FixedVariablesSet(convert(FixedVariablesIdeal{V, T, MT}, set.ideal)) +function Base.convert( + ::Type{FixedVariablesSet{V,T,MT}}, + set::FixedVariablesSet{V,S,MT}, +) where {V,S,T,MT} + return FixedVariablesSet(convert(FixedVariablesIdeal{V,T,MT}, set.ideal)) end ideal(set::FixedVariablesSet, args...) = set.ideal @@ -54,14 +75,17 @@ function nequalities(set::FixedVariablesSet) return length(set.ideal.substitutions) end end -function equalities(set::FixedVariablesSet{V, T, MT}) where {V, T, MT} +function equalities(set::FixedVariablesSet{V,T,MT}) where {V,T,MT} if set.ideal.substitutions === nothing return [constant_term(one(T), MT)] else return [key - value for (key, value) in set.ideal.substitutions] end end -function Base.intersect(V::FixedVariablesSet{V1, T1, MT1}, W::FixedVariablesSet{V2, T2, MT2}) where {V1, V2, T1, T2, MT1, MT2} +function Base.intersect( + V::FixedVariablesSet{V1,T1,MT1}, + W::FixedVariablesSet{V2,T2,MT2}, +) where {V1,V2,T1,T2,MT1,MT2} # For `DynamicPolynomials`, they have the same type and for # `TypedPolynomials`, promoting would give `Monomial`. VT = V1 == V2 ? V1 : AbstractVariable @@ -76,14 +100,14 @@ function Base.intersect(V::FixedVariablesSet{V1, T1, MT1}, W::FixedVariablesSet{ end return a end - sub = Dict{VT, T}() + sub = Dict{VT,T}() merge!(combine, sub, V.ideal.substitutions) merge!(combine, sub, W.ideal.substitutions) if has_dup sub = nothing end end - ideal = FixedVariablesIdeal{VT, T, promote_type(MT1, MT2)}(sub) + ideal = FixedVariablesIdeal{VT,T,promote_type(MT1, MT2)}(sub) return FixedVariablesSet(ideal) end function Base.intersect(V::AlgebraicSet, W::FixedVariablesSet) @@ -103,12 +127,12 @@ end # Assumes `isempty(V)` function only_point(V::FixedVariablesSet) subs = collect(V.ideal.substitutions) - sort!(subs, rev = true, by = sub -> sub[1]) + sort!(subs; rev = true, by = sub -> sub[1]) return [sub[2] for sub in subs] end -Base.eltype(::FixedVariablesSet{V, T}) where {V, T} = Vector{T} -function Base.iterate(V::FixedVariablesSet, state=nothing) +Base.eltype(::FixedVariablesSet{V,T}) where {V,T} = Vector{T} +function Base.iterate(V::FixedVariablesSet, state = nothing) if state === nothing && !isempty(V) return only_point(V), true else @@ -117,7 +141,7 @@ function Base.iterate(V::FixedVariablesSet, state=nothing) end Base.length(V::FixedVariablesSet) = isempty(V) ? 0 : 1 -struct FixedVariable{V<:AbstractVariable, T} +struct FixedVariable{V<:AbstractVariable,T} variable::V value::T end @@ -130,6 +154,13 @@ end function Base.intersect(el::FixedVariable; kws...) subs = Dict(el.variable => el.value) - return FixedVariablesSet(FixedVariablesIdeal{ - typeof(el.variable), typeof(el.value), typeof(el.variable)}(subs)) + return FixedVariablesSet( + FixedVariablesIdeal{ + typeof(el.variable), + typeof(el.value), + typeof(el.variable), + }( + subs, + ), + ) end diff --git a/src/groebner.jl b/src/groebner.jl index a84f47d..c653ca3 100644 --- a/src/groebner.jl +++ b/src/groebner.jl @@ -52,9 +52,11 @@ lcmlm(F, i, j) = lcm(leading_monomial(F[i]), leading_monomial(F[j])) function criterion(F, B, i, j) m = lcmlm(F, i, j) for k in eachindex(F) - if k != i && k != j && - !(ext(i, k) in B) && !(ext(j, k) in B) && - divides(leading_monomial(F[k]), m) + if k != i && + k != j && + !(ext(i, k) in B) && + !(ext(j, k) in B) && + divides(leading_monomial(F[k]), m) return true end end @@ -71,13 +73,13 @@ end # Ideals, Varieties, and Algorithms # Cox, Little and O'Shea, Fourth edition p. 116 function presort!(F) - sort!(F, by=leading_monomial) + return sort!(F; by = leading_monomial) end # Selection of (i, j) ∈ B function dummyselection(F, B) - first(B) + return first(B) end # "BUCHBERGER (1985) suggests that there will often be some savings if we pick # (i, j) ∈ B such that lcm(LM(fi), LM(fj)) is as small as possible. @@ -99,7 +101,7 @@ function normalselection(F, B) m = cur end end - best + return best end # TODO sugar and double sugar selections @@ -112,14 +114,24 @@ struct Buchberger <: AbstractGröbnerBasisAlgorithm sel::Function end #Buchberger() = Buchberger(presort!, dummyselection) -Buchberger(ztol=Base.rtoldefault(Float64)) = Buchberger(ztol, presort!, normalselection) +function Buchberger(ztol = Base.rtoldefault(Float64)) + return Buchberger(ztol, presort!, normalselection) +end # Taken from Theorem 9 of # Ideals, Varieties, and Algorithms # Cox, Little and O'Shea, Fourth edition -function gröbnerbasis!(F::AbstractVector{<:APL}, algo=defaultgröbnerbasisalgorithm(F)) +function gröbnerbasis!( + F::AbstractVector{<:APL}, + algo = defaultgröbnerbasisalgorithm(F), +) algo.pre!(F) - B = Set{NTuple{2, Int}}(Iterators.filter(t -> t[1] < t[2], (i, j) for i in eachindex(F), j in eachindex(F))) + B = Set{NTuple{2,Int}}( + Iterators.filter( + t -> t[1] < t[2], + (i, j) for i in eachindex(F), j in eachindex(F) + ), + ) while !isempty(B) (i, j) = algo.sel(F, B) lmi = leading_monomial(F[i]) @@ -127,7 +139,7 @@ function gröbnerbasis!(F::AbstractVector{<:APL}, algo=defaultgröbnerbasisalgor if !isconstant(gcd(lmi, lmj)) && !criterion(F, B, i, j) #r = rem(spolynomial(F[i], F[j]), F; ztol=algo.ztol) r = rem(spolynomial(F[i], F[j]), F) - if !isapproxzero(r; ztol=algo.ztol) + if !isapproxzero(r; ztol = algo.ztol) I = eachindex(F) push!(F, r) n = last(eachindex(F)) @@ -138,12 +150,12 @@ function gröbnerbasis!(F::AbstractVector{<:APL}, algo=defaultgröbnerbasisalgor end delete!(B, (i, j)) end - reducebasis!(F, ztol=algo.ztol) + reducebasis!(F; ztol = algo.ztol) map!(monic, F, F) - F + return F end function gröbnerbasis(F::Vector{<:APL}, args...) T = Base.promote_op(rem, eltype(F), typeof(F)) - gröbnerbasis!(copyto!(similar(F, T), F), args...) + return gröbnerbasis!(copyto!(similar(F, T), F), args...) end const groebnerbasis = gröbnerbasis diff --git a/src/ideal.jl b/src/ideal.jl index 8cf9ab3..27d5b4f 100644 --- a/src/ideal.jl +++ b/src/ideal.jl @@ -3,7 +3,9 @@ export monomialbasis, ideal abstract type AbstractPolynomialIdeal end struct EmptyPolynomialIdeal <: AbstractPolynomialIdeal end -ideal(p::FullSpace, algo=defaultgröbnerbasisalgorithm(p)) = EmptyPolynomialIdeal() +function ideal(p::FullSpace, algo = defaultgröbnerbasisalgorithm(p)) + return EmptyPolynomialIdeal() +end Base.rem(p::AbstractPolynomialLike, I::EmptyPolynomialIdeal) = p promote_for_division(::Type{T}) where {T} = T @@ -13,39 +15,53 @@ function promote_for_division(::Type{Complex{T}}) where {T} return Complex{promote_for_division(T)} end -mutable struct PolynomialIdeal{T, PT<:APL{T}, A<:AbstractGröbnerBasisAlgorithm} <: AbstractPolynomialIdeal +mutable struct PolynomialIdeal{T,PT<:APL{T},A<:AbstractGröbnerBasisAlgorithm} <: + AbstractPolynomialIdeal p::Vector{PT} gröbnerbasis::Bool algo::A end -function PolynomialIdeal{T, PT}(p::Vector{PT}, algo::A) where {T, PT<:APL{T}, A<:AbstractGröbnerBasisAlgorithm} - PolynomialIdeal{T, PT, A}(p, false, algo) +function PolynomialIdeal{T,PT}( + p::Vector{PT}, + algo::A, +) where {T,PT<:APL{T},A<:AbstractGröbnerBasisAlgorithm} + return PolynomialIdeal{T,PT,A}(p, false, algo) end -function PolynomialIdeal(p::Vector{PT}, algo) where {T, PT<:APL{T}} +function PolynomialIdeal(p::Vector{PT}, algo) where {T,PT<:APL{T}} S = promote_for_division(T) - PolynomialIdeal{S, polynomial_type(PT, S)}(AbstractVector{polynomial_type(PT, S)}(p), algo) + return PolynomialIdeal{S,polynomial_type(PT, S)}( + AbstractVector{polynomial_type(PT, S)}(p), + algo, + ) end -function PolynomialIdeal{T, PT}() where {T, PT<:APL{T}} - PolynomialIdeal(PT[], defaultgröbnerbasisalgorithm(PT)) +function PolynomialIdeal{T,PT}() where {T,PT<:APL{T}} + return PolynomialIdeal(PT[], defaultgröbnerbasisalgorithm(PT)) end -function Base.convert(::Type{PolynomialIdeal{T, PT, A}}, I::PolynomialIdeal{T, PT, A}) where {T, PT<:APL{T}, A<:AbstractGröbnerBasisAlgorithm} +function Base.convert( + ::Type{PolynomialIdeal{T,PT,A}}, + I::PolynomialIdeal{T,PT,A}, +) where {T,PT<:APL{T},A<:AbstractGröbnerBasisAlgorithm} return I end -function Base.convert(::Type{PolynomialIdeal{T, PT, A}}, I::PolynomialIdeal) where {T, PT, A} - return PolynomialIdeal{T, PT, A}(I.p, I.gröbnerbasis, I.algo) +function Base.convert( + ::Type{PolynomialIdeal{T,PT,A}}, + I::PolynomialIdeal, +) where {T,PT,A} + return PolynomialIdeal{T,PT,A}(I.p, I.gröbnerbasis, I.algo) end +ideal(p, algo = defaultgröbnerbasisalgorithm(p)) = PolynomialIdeal(p, algo) -ideal(p, algo=defaultgröbnerbasisalgorithm(p)) = PolynomialIdeal(p, algo) - -Base.:+(I::PolynomialIdeal, J::PolynomialIdeal) = PolynomialIdeal([I.p; J.p], I.algo) +function Base.:+(I::PolynomialIdeal, J::PolynomialIdeal) + return PolynomialIdeal([I.p; J.p], I.algo) +end MP.variables(I::PolynomialIdeal) = variables(I.p) function Base.rem(p::AbstractPolynomialLike, I::PolynomialIdeal) computegröbnerbasis!(I) - rem(p, I.p) + return rem(p, I.p) end function computegröbnerbasis!(I::PolynomialIdeal) @@ -54,7 +70,7 @@ function computegröbnerbasis!(I::PolynomialIdeal) I.gröbnerbasis = true end end -function monomialbasis(I, vars=variables(I)) +function monomialbasis(I, vars = variables(I)) computegröbnerbasis!(I) if isempty(I.p) return false, monomial_type(eltype(I.p))[] @@ -77,5 +93,6 @@ function monomialbasis(I, vars=variables(I)) if any(lv .< 0) return false, monomial_type(eltype(I.p))[] end - return true, monomials(vars, 0:(sum(lv)-1), m -> !any(map(m2 -> divides(m2, m), mv))) + return true, + monomials(vars, 0:(sum(lv)-1), m -> !any(map(m2 -> divides(m2, m), mv))) end diff --git a/src/macro.jl b/src/macro.jl index 58f0eb0..5b14d10 100644 --- a/src/macro.jl +++ b/src/macro.jl @@ -13,7 +13,7 @@ end function equality(lhs, rhs) return PolynomialEquality(lhs - rhs) end -function Base.intersect(eq::PolynomialEquality; lib_or_solver=nothing) +function Base.intersect(eq::PolynomialEquality; lib_or_solver = nothing) if lib_or_solver === nothing return algebraicset([eq.p]) else @@ -21,32 +21,35 @@ function Base.intersect(eq::PolynomialEquality; lib_or_solver=nothing) end end -const Element = Union{PolynomialInequality, PolynomialEquality, FixedVariable} +const Element = Union{PolynomialInequality,PolynomialEquality,FixedVariable} -function Base.intersect(el::Element, - args...; kws...) +function Base.intersect(el::Element, args...; kws...) return intersect(intersect(el; kws...), args...; kws...) end -function Base.intersect(set::AbstractBasicSemialgebraicSet, - el::Element, args...; kws...) +function Base.intersect( + set::AbstractBasicSemialgebraicSet, + el::Element, + args...; + kws..., +) return intersect(set, intersect(el; kws...), args...; kws...) end # Taken from JuMP/macros.jl function _canonicalize_sense(sns::Symbol, _error) if sns == :(==) - return (:(==),false) + return (:(==), false) elseif sns == :(>=) || sns == :(≥) - return (:(>=),false) + return (:(>=), false) elseif sns == :(<=) || sns == :(≤) - return (:(<=),false) + return (:(<=), false) elseif sns == :(.==) - return (:(==),true) + return (:(==), true) elseif sns == :(.>=) || sns == :(.≥) - return (:(>=),true) + return (:(>=), true) elseif sns == :(.<=) || sns == :(.≤) - return (:(<=),true) + return (:(<=), true) else _error("Unrecognized sense $sns") end @@ -58,13 +61,37 @@ function appendconstraints!(elements, expr, _error) sense, vectorized = _canonicalize_sense(expr.args[1], _error) @assert !vectorized if sense == :(>=) - push!(elements, esc(:(SemialgebraicSets.PolynomialInequality($(expr.args[2]) - $(expr.args[3]))))) + push!( + elements, + esc( + :(SemialgebraicSets.PolynomialInequality( + $(expr.args[2]) - $(expr.args[3]), + )), + ), + ) elseif sense == :(<=) - push!(elements, esc(:(SemialgebraicSets.PolynomialInequality($(expr.args[3]) - $(expr.args[2]))))) + push!( + elements, + esc( + :(SemialgebraicSets.PolynomialInequality( + $(expr.args[3]) - $(expr.args[2]), + )), + ), + ) elseif sense == :(==) - push!(elements, esc(:(SemialgebraicSets.equality($(expr.args[2]), $(expr.args[3]))))) + push!( + elements, + esc( + :(SemialgebraicSets.equality( + $(expr.args[2]), + $(expr.args[3]), + )), + ), + ) else - _error("Unrecognized sense $(string(sense)) in domain specification") + _error( + "Unrecognized sense $(string(sense)) in domain specification", + ) end catch push!(elements, esc(expr)) @@ -77,12 +104,12 @@ function appendconstraints!(elements, expr, _error) return nothing end -macro set(expr, library=nothing) +macro set(expr, library = nothing) elements = [] appendconstraints!(elements, expr, msg -> error("In @set($expr: ", msg)) if library === nothing - return :( intersect($(elements...)) ) + return :(intersect($(elements...))) else - return :( intersect($(elements...); lib_or_solver=$(esc(library))) ) + return :(intersect($(elements...); lib_or_solver = $(esc(library)))) end end diff --git a/src/schur.jl b/src/schur.jl index 94d4cd3..e67a50f 100644 --- a/src/schur.jl +++ b/src/schur.jl @@ -8,7 +8,13 @@ function conditionnumber(sf::Schur, I) for i in I select[i] = 1 end - LinearAlgebra.LAPACK.trsen!('E', 'N', select, copy(sf.T), copy(sf.Z))[4] + return LinearAlgebra.LAPACK.trsen!( + 'E', + 'N', + select, + copy(sf.T), + copy(sf.Z), + )[4] end # Manocha, D. & Demmel, J. Algorithms for intersecting parametric and algebraic curves II: multiple intersections @@ -54,10 +60,10 @@ function _clusterordschur(M::AbstractMatrix{<:Real}, ɛ) i += 1 else @assert i < lastindex(v) && !isreal(v[i+1]) - pairatol = _atol([i, i+1]) + pairatol = _atol([i, i + 1]) if abs(v[i] - v[i+1]) / pairatol < ONE # Pair conjugate pairs into a real eigenvalue - push!(clusters, [i, i+1]) + push!(clusters, [i, i + 1]) push!(λ, real((v[i] + v[i+1]) / 2)) # The imaginary part should be zero anyway push!(atol, pairatol) end diff --git a/src/solve.jl b/src/solve.jl index f87bdd7..000dd10 100644 --- a/src/solve.jl +++ b/src/solve.jl @@ -13,7 +13,9 @@ abstract type AbstractAlgebraicSolver end Solve the algebraic equations for which `V` is the set of solutions using the algorithm `algo`. Returns a nullable which is `null` if `V` is not zero-dimensional and is the list of solutions otherwise. """ -solvealgebraicequations(V::AbstractAlgebraicSet) = solvealgebraicequations(V, defaultalgebraicsolver(V)) +function solvealgebraicequations(V::AbstractAlgebraicSet) + return solvealgebraicequations(V, defaultalgebraicsolver(V)) +end """ AbstractMultiplicationMatricesAlgorithm @@ -45,13 +47,18 @@ Returns the list of solutions. """ function solvemultiplicationmatrices end - -struct SolverUsingMultiplicationMatrices{A<:AbstractMultiplicationMatricesAlgorithm, S<:AbstractMultiplicationMatricesSolver} <: AbstractAlgebraicSolver +struct SolverUsingMultiplicationMatrices{ + A<:AbstractMultiplicationMatricesAlgorithm, + S<:AbstractMultiplicationMatricesSolver, +} <: AbstractAlgebraicSolver algo::A solver::S end -function solvealgebraicequations(V::AbstractAlgebraicSet, solver::SolverUsingMultiplicationMatrices) +function solvealgebraicequations( + V::AbstractAlgebraicSet, + solver::SolverUsingMultiplicationMatrices, +) Ms = multiplicationmatrices(V, solver.algo) if Ms === nothing nothing @@ -60,8 +67,8 @@ function solvealgebraicequations(V::AbstractAlgebraicSet, solver::SolverUsingMul end end -struct GröbnerBasisMultiplicationMatricesAlgorithm <: AbstractMultiplicationMatricesAlgorithm -end +struct GröbnerBasisMultiplicationMatricesAlgorithm <: + AbstractMultiplicationMatricesAlgorithm end function multiplicationmatrix(V::AbstractAlgebraicSet, v::AbstractVariable, B) M = Matrix{eltype(eltype(V))}(undef, length(B), length(B)) @@ -69,10 +76,13 @@ function multiplicationmatrix(V::AbstractAlgebraicSet, v::AbstractVariable, B) p = rem(v * B[i], equalities(V)) M[:, i] = coefficients(p, B) end - M + return M end -function multiplicationmatrices(V::AbstractAlgebraicSet, algo::GröbnerBasisMultiplicationMatricesAlgorithm) +function multiplicationmatrices( + V::AbstractAlgebraicSet, + algo::GröbnerBasisMultiplicationMatricesAlgorithm, +) vars = variables(V.I) iszd, B = monomialbasis(V.I, vars) if !iszd @@ -92,21 +102,33 @@ include("schur.jl") """ Corless, R. M.; Gianni, P. M. & Trager, B. M. A reordered Schur factorization method for zero-dimensional polynomial systems with multiple roots Proceedings of the 1997 international symposium on Symbolic and algebraic computation, 1997, 133-140 """ -struct ReorderedSchurMultiplicationMatricesSolver{T, RNGT <: Random.AbstractRNG} <: AbstractMultiplicationMatricesSolver +struct ReorderedSchurMultiplicationMatricesSolver{T,RNGT<:Random.AbstractRNG} <: + AbstractMultiplicationMatricesSolver ɛ::T rng::RNGT end -ReorderedSchurMultiplicationMatricesSolver(ɛ) = ReorderedSchurMultiplicationMatricesSolver(ɛ, Random.GLOBAL_RNG) -ReorderedSchurMultiplicationMatricesSolver{T}() where T = ReorderedSchurMultiplicationMatricesSolver(Base.rtoldefault(real(T))) +function ReorderedSchurMultiplicationMatricesSolver(ɛ) + return ReorderedSchurMultiplicationMatricesSolver(ɛ, Random.GLOBAL_RNG) +end +function ReorderedSchurMultiplicationMatricesSolver{T}() where {T} + return ReorderedSchurMultiplicationMatricesSolver(Base.rtoldefault(real(T))) +end -function solvemultiplicationmatrices(Ms::AbstractVector{<:AbstractMatrix{T}}, solver::ReorderedSchurMultiplicationMatricesSolver) where T +function solvemultiplicationmatrices( + Ms::AbstractVector{<:AbstractMatrix{T}}, + solver::ReorderedSchurMultiplicationMatricesSolver, +) where {T} λ = rand(solver.rng, length(Ms)) λ /= sum(λ) - _solvemultiplicationmatrices(Ms, λ, solver) + return _solvemultiplicationmatrices(Ms, λ, solver) end # Deterministic part -function _solvemultiplicationmatrices(Ms::AbstractVector{<:AbstractMatrix{T}}, λ, solver::ReorderedSchurMultiplicationMatricesSolver) where T<:Real +function _solvemultiplicationmatrices( + Ms::AbstractVector{<:AbstractMatrix{T}}, + λ, + solver::ReorderedSchurMultiplicationMatricesSolver, +) where {T<:Real} @assert length(Ms) == length(λ) n = length(λ) Z, clusters = clusterordschur(sum(λ .* Ms), solver.ɛ) @@ -124,24 +146,40 @@ function _solvemultiplicationmatrices(Ms::AbstractVector{<:AbstractMatrix{T}}, end end end - vals + return vals end -function algebraicsolver(algo::AbstractMultiplicationMatricesAlgorithm, - solver::AbstractMultiplicationMatricesSolver) - SolverUsingMultiplicationMatrices(algo, solver) +function algebraicsolver( + algo::AbstractMultiplicationMatricesAlgorithm, + solver::AbstractMultiplicationMatricesSolver, +) + return SolverUsingMultiplicationMatrices(algo, solver) end -defaultmultiplicationmatricesalgorithm(p) = GröbnerBasisMultiplicationMatricesAlgorithm() -defaultmultiplicationmatricessolver(::Type{T}) where T = ReorderedSchurMultiplicationMatricesSolver{T}() -defaultmultiplicationmatricessolver(::AbstractVector{PT}) where {T, PT<:APL{T}} = defaultmultiplicationmatricessolver(T) +function defaultmultiplicationmatricesalgorithm(p) + return GröbnerBasisMultiplicationMatricesAlgorithm() +end +function defaultmultiplicationmatricessolver(::Type{T}) where {T} + return ReorderedSchurMultiplicationMatricesSolver{T}() +end +function defaultmultiplicationmatricessolver( + ::AbstractVector{PT}, +) where {T,PT<:APL{T}} + return defaultmultiplicationmatricessolver(T) +end function defaultalgebraicsolver(p) - algebraicsolver(defaultmultiplicationmatricesalgorithm(p), defaultmultiplicationmatricessolver(p)) + return algebraicsolver( + defaultmultiplicationmatricesalgorithm(p), + defaultmultiplicationmatricessolver(p), + ) end -function defaultalgebraicsolver(p, algo::AbstractMultiplicationMatricesAlgorithm) - algebraicsolver(algo, defaultmultiplicationmatricessolver(p)) +function defaultalgebraicsolver( + p, + algo::AbstractMultiplicationMatricesAlgorithm, +) + return algebraicsolver(algo, defaultmultiplicationmatricessolver(p)) end function defaultalgebraicsolver(p, solver::AbstractMultiplicationMatricesSolver) - algebraicsolver(defaultmultiplicationmatricesalgorithm(p), solver) + return algebraicsolver(defaultmultiplicationmatricesalgorithm(p), solver) end diff --git a/src/variety.jl b/src/variety.jl index 8454dee..253c958 100644 --- a/src/variety.jl +++ b/src/variety.jl @@ -5,43 +5,117 @@ struct DefaultAlgebraicSetLibrary{S<:AbstractAlgebraicSolver} solver::S end -defaultalgebraicsetlibrary(::Vector{<:APL}, solver::AbstractAlgebraicSolver) = DefaultAlgebraicSetLibrary(solver) -defaultalgebraicsetlibrary(p::Vector{<:APL}, solveroralgo...) = defaultalgebraicsetlibrary(p, defaultalgebraicsolver(p, solveroralgo...)) +function defaultalgebraicsetlibrary( + ::Vector{<:APL}, + solver::AbstractAlgebraicSolver, +) + return DefaultAlgebraicSetLibrary(solver) +end +function defaultalgebraicsetlibrary(p::Vector{<:APL}, solveroralgo...) + return defaultalgebraicsetlibrary( + p, + defaultalgebraicsolver(p, solveroralgo...), + ) +end -mutable struct AlgebraicSet{T, PT<:APL{T}, A, S<:AbstractAlgebraicSolver, U} <: AbstractAlgebraicSet - I::PolynomialIdeal{T, PT, A} +mutable struct AlgebraicSet{T,PT<:APL{T},A,S<:AbstractAlgebraicSolver,U} <: + AbstractAlgebraicSet + I::PolynomialIdeal{T,PT,A} projective::Bool elements::Vector{Vector{U}} elementscomputed::Bool iszerodimensional::Bool solver::S end -AlgebraicSet{T, PT, A, S, U}(I::PolynomialIdeal{T, PT, A}, solver::S) where {T, PT, A, S, U} = AlgebraicSet{T, PT, A, S, U}(I, false, Vector{U}[], false, false, solver) -AlgebraicSet(I::PolynomialIdeal{T, PT, A}, solver::S) where {T, PT, A, S} = AlgebraicSet{T, PT, A, S, float(T)}(I, solver) -AlgebraicSet{T, PT}() where {T, PT} = AlgebraicSet(PolynomialIdeal{T, PT}(), defaultalgebraicsolver(T)) -AlgebraicSet(p::Vector, algo::AbstractGröbnerBasisAlgorithm, solver) = AlgebraicSet(ideal(p, algo), solver) +function AlgebraicSet{T,PT,A,S,U}( + I::PolynomialIdeal{T,PT,A}, + solver::S, +) where {T,PT,A,S,U} + return AlgebraicSet{T,PT,A,S,U}(I, false, Vector{U}[], false, false, solver) +end +function AlgebraicSet(I::PolynomialIdeal{T,PT,A}, solver::S) where {T,PT,A,S} + return AlgebraicSet{T,PT,A,S,float(T)}(I, solver) +end +function AlgebraicSet{T,PT}() where {T,PT} + return AlgebraicSet(PolynomialIdeal{T,PT}(), defaultalgebraicsolver(T)) +end +function AlgebraicSet(p::Vector, algo::AbstractGröbnerBasisAlgorithm, solver) + return AlgebraicSet(ideal(p, algo), solver) +end -MP.similar_type(::Type{AlgebraicSet{U, PU, A, S, UU}}, T::Type) where {U, PU, A, S, UU} = AlgebraicSet{T, MP.similar_type(PU, T), A, S, float(T)} -function Base.convert(::Type{AlgebraicSet{T, PT, A, S, U}}, set::AlgebraicSet{T, PT, A, S, U}) where {T, PT<:APL{T}, A, S<:AbstractAlgebraicSolver, U} +function MP.similar_type( + ::Type{AlgebraicSet{U,PU,A,S,UU}}, + T::Type, +) where {U,PU,A,S,UU} + return AlgebraicSet{T,MP.similar_type(PU, T),A,S,float(T)} +end +function Base.convert( + ::Type{AlgebraicSet{T,PT,A,S,U}}, + set::AlgebraicSet{T,PT,A,S,U}, +) where {T,PT<:APL{T},A,S<:AbstractAlgebraicSolver,U} return set end -function Base.convert(::Type{AlgebraicSet{T, PT, A, S, U}}, set::AlgebraicSet) where {T, PT, A, S, U} - return AlgebraicSet{T, PT, A, S, U}(set.I, set.projective, set.elements, set.elementscomputed, set.iszerodimensional, set.solver) +function Base.convert( + ::Type{AlgebraicSet{T,PT,A,S,U}}, + set::AlgebraicSet, +) where {T,PT,A,S,U} + return AlgebraicSet{T,PT,A,S,U}( + set.I, + set.projective, + set.elements, + set.elementscomputed, + set.iszerodimensional, + set.solver, + ) end -algebraicset(p::Vector, lib::DefaultAlgebraicSetLibrary) = AlgebraicSet(p, defaultgröbnerbasisalgorithm(p), lib.solver) -algebraicset(p::Vector, algo::AbstractGröbnerBasisAlgorithm=defaultgröbnerbasisalgorithm(p), lib::DefaultAlgebraicSetLibrary=defaultalgebraicsetlibrary(p)) = AlgebraicSet(p, algo, lib.solver) -algebraicset(p::Vector, solver) = algebraicset(p, defaultalgebraicsetlibrary(p, solver)) -algebraicset(p::Vector, algo::AbstractGröbnerBasisAlgorithm, solver) = algebraicset(p, algo, defaultalgebraicsetlibrary(p, solver)) +function algebraicset(p::Vector, lib::DefaultAlgebraicSetLibrary) + return AlgebraicSet(p, defaultgröbnerbasisalgorithm(p), lib.solver) +end +function algebraicset( + p::Vector, + algo::AbstractGröbnerBasisAlgorithm = defaultgröbnerbasisalgorithm(p), + lib::DefaultAlgebraicSetLibrary = defaultalgebraicsetlibrary(p), +) + return AlgebraicSet(p, algo, lib.solver) +end +function algebraicset(p::Vector, solver) + return algebraicset(p, defaultalgebraicsetlibrary(p, solver)) +end +function algebraicset(p::Vector, algo::AbstractGröbnerBasisAlgorithm, solver) + return algebraicset(p, algo, defaultalgebraicsetlibrary(p, solver)) +end -projectivealgebraicset(p::Vector, lib::DefaultAlgebraicSetLibrary) = projectivealgebraicset(p, defaultgröbnerbasisalgorithm(p), lib.solver) -function projectivealgebraicset(p::Vector, algo::AbstractGröbnerBasisAlgorithm=defaultgröbnerbasisalgorithm(p), lib::DefaultAlgebraicSetLibrary=defaultalgebraicsetlibrary(p)) +function projectivealgebraicset(p::Vector, lib::DefaultAlgebraicSetLibrary) + return projectivealgebraicset( + p, + defaultgröbnerbasisalgorithm(p), + lib.solver, + ) +end +function projectivealgebraicset( + p::Vector, + algo::AbstractGröbnerBasisAlgorithm = defaultgröbnerbasisalgorithm(p), + lib::DefaultAlgebraicSetLibrary = defaultalgebraicsetlibrary(p), +) V = AlgebraicSet(p, algo, lib.solver) V.projective = true - V + return V +end +function projectivealgebraicset(p::Vector, algo, solver) + return projectivealgebraicset( + p, + algo, + defaultalgebraicsetlibrary(p, solver), + ) +end +function projectivealgebraicset(p::Vector, solver) + return projectivealgebraicset( + p, + defaultgröbnerbasisalgorithm(p), + defaultalgebraicsetlibrary(p, solver), + ) end -projectivealgebraicset(p::Vector, algo, solver) = projectivealgebraicset(p, algo, defaultalgebraicsetlibrary(p, solver)) -projectivealgebraicset(p::Vector, solver) = projectivealgebraicset(p, defaultgröbnerbasisalgorithm(p), defaultalgebraicsetlibrary(p, solver)) ideal(V::AlgebraicSet) = V.I @@ -49,14 +123,22 @@ MP.variables(V::AlgebraicSet) = MP.variables(V.I) nequalities(V::AlgebraicSet) = length(V.I.p) equalities(V::AlgebraicSet) = V.I.p addequality!(V::AlgebraicSet, p) = push!(V.I.p, p) -Base.intersect(S::AlgebraicSet, T::AlgebraicSet) = AlgebraicSet(S.I + T.I, S.solver) +function Base.intersect(S::AlgebraicSet, T::AlgebraicSet) + return AlgebraicSet(S.I + T.I, S.solver) +end function algebraicset(set::AbstractAlgebraicSet, solver...) return algebraicset(equalities(set), solver...) end function Base.show(io::IO, V::AbstractAlgebraicSet) - print(io, "{ (", join(variables(V), ", "), ") | ", - join(string.(equalities(V)) .* " = 0", ", "), " }") + return print( + io, + "{ (", + join(variables(V), ", "), + ") | ", + join(string.(equalities(V)) .* " = 0", ", "), + " }", + ) end function _show_els(io::IO, name, n, els, sign) if n == 0 @@ -72,11 +154,11 @@ function _show_els(io::IO, name, n, els, sign) end function Base.show(io::IO, mime::MIME"text/plain", V::AbstractAlgebraicSet) print(io, "Algebraic Set defined by ") - _show_els(io, "equalit", nequalities(V), equalities(V), "=") + return _show_els(io, "equalit", nequalities(V), equalities(V), "=") end defaultalgebraicsolver(V::AlgebraicSet) = V.solver -function elements(V::AlgebraicSet{T, PT, A, S, U}) where {T, PT, A, S, U} +function elements(V::AlgebraicSet{T,PT,A,S,U}) where {T,PT,A,S,U} if V.projective I = V.I els = Vector{U}[] @@ -94,7 +176,7 @@ function elements(V::AlgebraicSet{T, PT, A, S, U}) where {T, PT, A, S, U} solvealgebraicequations(V, V.solver) end end -function computeelements!(V::AlgebraicSet{T}) where T +function computeelements!(V::AlgebraicSet{T}) where {T} if !V.elementscomputed els = elements(V) V.iszerodimensional = els !== nothing @@ -106,10 +188,10 @@ function computeelements!(V::AlgebraicSet{T}) where T end function iszerodimensional(V::AlgebraicSet) computeelements!(V) - V.iszerodimensional + return V.iszerodimensional end -Base.eltype(V::AlgebraicSet{T, PT, A, S, U}) where {T, PT, A, S, U} = Vector{U} +Base.eltype(V::AlgebraicSet{T,PT,A,S,U}) where {T,PT,A,S,U} = Vector{U} for f in [:length, :iterate, :lastindex, :getindex] @eval begin @@ -118,7 +200,7 @@ for f in [:length, :iterate, :lastindex, :getindex] if !iszerodimensional(V) error("A non zero-dimensional algebraic set is not iterable") end - $f(V.elements, args...) + return $f(V.elements, args...) end end end diff --git a/test/bug.jl b/test/bug.jl index d921ff9..59d6dce 100644 --- a/test/bug.jl +++ b/test/bug.jl @@ -1,12 +1,15 @@ -function test(ztol=Base.rtoldefault(Float64)) +function test(ztol = Base.rtoldefault(Float64)) X = [0.5459627556242905, 1.7288950489507429, 0.7167681447476535] Y = [-54.06002080721971, 173.77393714162503, 71.48154370522498] n = 3 @polyvar W[1:n] α β - I = @set(sum([W[i] * X[i] * (β * X[i] + α - Y[i]) for i = 1:n]) == 0, library = Buchberger(ztol)) + I = @set( + sum([W[i] * X[i] * (β * X[i] + α - Y[i]) for i in 1:n]) == 0, + library = Buchberger(ztol) + ) for i in 1:n addequality!(I, W[i] - W[i]^2) end SemialgebraicSets.computegröbnerbasis!(I.I) - I + return I end diff --git a/test/groebner.jl b/test/groebner.jl index 33e9df4..b45b463 100644 --- a/test/groebner.jl +++ b/test/groebner.jl @@ -13,9 +13,10 @@ @testset "S-polynomial" begin Mod.@polyvar x y z - @test spolynomial(x^3*y^2 - x^2*y^3 + x, 3x^4*y + y^2) == -x^3*y^3 + x^2 - y^3/3 - @test spolynomial(y - x^2, z - x^3) == z - x*y - @test spolynomial(x^3 - 2x*y, x^2*y - 2y^2 + x) == -x^2 + @test spolynomial(x^3 * y^2 - x^2 * y^3 + x, 3x^4 * y + y^2) == + -x^3 * y^3 + x^2 - y^3 / 3 + @test spolynomial(y - x^2, z - x^3) == z - x * y + @test spolynomial(x^3 - 2x * y, x^2 * y - 2y^2 + x) == -x^2 end @testset "Groebner basis" begin @@ -26,21 +27,26 @@ end @test all(_isz.(G - H)) end function testf(F, H) - testg(gröbnerbasis(F), H) + return testg(gröbnerbasis(F), H) end testf([4x^2 + 5x, 3x^3], [x]) - testf([y - x^2, z - x^3], [x*z - y^2, x*y - z, x^2 - y, y^3 - z^2]) - testf([x^3 - 2x*y, x^2*y - 2y^2 + x], [y^2 - 0.5x, x*y, x^2]) + testf([y - x^2, z - x^3], [x * z - y^2, x * y - z, x^2 - y, y^3 - z^2]) + testf([x^3 - 2x * y, x^2 * y - 2y^2 + x], [y^2 - 0.5x, x * y, x^2]) #p = 9x^11 + 27x^14 + x # root x=>-0.83005, y=>-3*(-0.83005)^4) - F = [x^3*y^2 - x^2*y^3 + x, 3x^4*y + y^2] - H = [x*y^3 - y^4 - 3x^3, - x^3*y^2 - x^2*y^3 + x, - x^4*y + y^2/3, - x^5 + x*y/3, - y^6 + 3y^5 + 9x^4 + 9x^3*y + y^3/3 - x^2 - x*y - y^2 - 3x] + F = [x^3 * y^2 - x^2 * y^3 + x, 3x^4 * y + y^2] + H = [ + x * y^3 - y^4 - 3x^3, + x^3 * y^2 - x^2 * y^3 + x, + x^4 * y + y^2 / 3, + x^5 + x * y / 3, + y^6 + 3y^5 + 9x^4 + 9x^3 * y + y^3 / 3 - x^2 - x * y - y^2 - 3x, + ] function testb(pre, sel) - testg(groebnerbasis(F, Buchberger(Base.rtoldefault(Float64), pre, sel)), H) + return testg( + groebnerbasis(F, Buchberger(Base.rtoldefault(Float64), pre, sel)), + H, + ) end testb(identity, dummyselection) testb(identity, normalselection) @@ -51,7 +57,7 @@ end @testset "Reduce" begin Mod.@polyvar x y V1 = @set x == 1 && y == x^2 - @test rem(x^2 + 3x*y + 2y, ideal(V1)) == 6 + @test rem(x^2 + 3x * y + 2y, ideal(V1)) == 6 p = x^2 + x V2 = FullSpace() @test rem(p, ideal(V2)) === p diff --git a/test/macro.jl b/test/macro.jl index 877a327..08c330b 100644 --- a/test/macro.jl +++ b/test/macro.jl @@ -5,37 +5,47 @@ struct DummySolver <: SemialgebraicSets.AbstractAlgebraicSolver end @test isa(FullSpace(), FullSpace) V = @set x * y == 1 @test V isa AlgebraicSet{Rational{BigInt}} - @test_throws ArgumentError addinequality!(V, x*y) + @test_throws ArgumentError addinequality!(V, x * y) @testset "Basic" begin - S = @set x - y == 0 && x^2*y >= 1 + S = @set x - y == 0 && x^2 * y >= 1 addequality!(S, x^2 - y) addinequality!(S, x + y - 1) # Algebraic set forces `Rational{BigInt}` @test S isa BasicSemialgebraicSet{Rational{BigInt}} @test S == basicsemialgebraicset(S.V, S.p) - @test sprint(show, S) == "{ (x, y) | -y + x = 0, -y + x^2 = 0, -1//1 + x^2*y ≥ 0, -1//1 + y + x ≥ 0 }" - @test sprint(show, MIME"text/plain"(), S) == "Basic semialgebraic Set defined by 2 equalities\n -y + x = 0\n -y + x^2 = 0\n2 inequalities\n -1//1 + x^2*y ≥ 0\n -1//1 + y + x ≥ 0\n" + @test sprint(show, S) == + "{ (x, y) | -y + x = 0, -y + x^2 = 0, -1//1 + x^2*y ≥ 0, -1//1 + y + x ≥ 0 }" + @test sprint(show, MIME"text/plain"(), S) == + "Basic semialgebraic Set defined by 2 equalities\n -y + x = 0\n -y + x^2 = 0\n2 inequalities\n -1//1 + x^2*y ≥ 0\n -1//1 + y + x ≥ 0\n" @test S.V isa AlgebraicSet{Rational{BigInt}} @test sprint(show, S.V) == "{ (x, y) | -y + x = 0, -y + x^2 = 0 }" - @test sprint(show, MIME"text/plain"(), S.V) == "Algebraic Set defined by 2 equalities\n -y + x = 0\n -y + x^2 = 0\n" + @test sprint(show, MIME"text/plain"(), S.V) == + "Algebraic Set defined by 2 equalities\n -y + x = 0\n -y + x^2 = 0\n" @test S === similar(S, Rational{BigInt}) @test S.V === similar(S.V, Rational{BigInt}) @test S.V.I === convert(typeof(S.V.I), S.V.I) - @test BasicSemialgebraicSet{Int, polynomial_type(x, Int)}() isa BasicSemialgebraicSet{Rational{BigInt}, polynomial_type(x, Rational{BigInt})} - @test Int32(2)*x^2*y isa MultivariatePolynomials.AbstractTerm{Int32} + @test BasicSemialgebraicSet{Int,polynomial_type(x, Int)}() isa + BasicSemialgebraicSet{ + Rational{BigInt}, + polynomial_type(x, Rational{BigInt}), + } + @test Int32(2) * x^2 * y isa MultivariatePolynomials.AbstractTerm{Int32} Sf = similar(S, Float32) @test Sf isa BasicSemialgebraicSet{Float32} @test Sf.V isa AlgebraicSet{Float32} @testset "Mixed types" begin - S = (@set Int32(2)*x^2*y == 0 && 1.0x^2*y >= 0 && (6//3)*x^2*y == -y && 1.5x+y >= 0) + S = (@set Int32(2) * x^2 * y == 0 && + 1.0x^2 * y >= 0 && + (6 // 3) * x^2 * y == -y && + 1.5x + y >= 0) S2 = S ∩ V S3 = V ∩ S @test inequalities(S2) == inequalities(S3) == S.p @test equalities(S2) == S3.V.I.p end - T = (@set x*y^2 == -1 && x^2 + y^2 <= 1) + T = (@set x * y^2 == -1 && x^2 + y^2 <= 1) V2 = @set T.V && V && x + y == 2.0 @test V2 isa AlgebraicSet @test V2.I.p == [equalities(T); equalities(V); x + y - 2.0] @@ -46,7 +56,8 @@ struct DummySolver <: SemialgebraicSets.AbstractAlgebraicSolver end @testset "Different variables" begin T = (@set x == x^2 && y <= y^2) @test sprint(show, T) == "{ (x, y) | x - x^2 = 0, -y + y^2 ≥ 0 }" - @test sprint(show, MIME"text/plain"(), T) == "Basic semialgebraic Set defined by 1 equalitty\n x - x^2 = 0\n1 inequalitty\n -y + y^2 ≥ 0\n" + @test sprint(show, MIME"text/plain"(), T) == + "Basic semialgebraic Set defined by 1 equalitty\n x - x^2 = 0\n1 inequalitty\n -y + y^2 ≥ 0\n" end end @testset "Basic with no equality" begin @@ -71,46 +82,50 @@ Algebraic Set defined by no equality @test S ∩ FullSpace() === S @test S.V ∩ FullSpace() === S.V - @test FullSpace() ∩ S === S - @test FullSpace() ∩ S.V === S.V + @test FullSpace() ∩ S === S + @test FullSpace() ∩ S.V === S.V @test S.V === similar(S.V, Float64) end @testset "Basic with fixed variables" begin S = @set x == 1 && x ≤ x^2 @test sprint(show, S) == "{ (x) | -1 + x = 0, -x + x^2 ≥ 0 }" - @test sprint(show, MIME"text/plain"(), S) == "Basic semialgebraic Set defined by 1 equalitty\n -1 + x = 0\n1 inequalitty\n -x + x^2 ≥ 0\n" + @test sprint(show, MIME"text/plain"(), S) == + "Basic semialgebraic Set defined by 1 equalitty\n -1 + x = 0\n1 inequalitty\n -x + x^2 ≥ 0\n" @test sprint(show, S.V) == "{ (x) | -1 + x = 0 }" - @test sprint(show, MIME"text/plain"(), S.V) == "Algebraic Set defined by 1 equalitty\n -1 + x = 0\n" + @test sprint(show, MIME"text/plain"(), S.V) == + "Algebraic Set defined by 1 equalitty\n -1 + x = 0\n" S = @set x == 1 && x ≤ y && 2 == y @test S isa BasicSemialgebraicSet{Int} - @test S.V isa FixedVariablesSet{<:AbstractVariable, Int} + @test S.V isa FixedVariablesSet{<:AbstractVariable,Int} @test rem(x + y, ideal(S.V)) == 3 @test S === similar(S, Int) @test S.V === similar(S.V, Int) @test S.V.ideal === convert(typeof(S.V.ideal), S.V.ideal) Sf = similar(S, Float64) @test Sf isa BasicSemialgebraicSet{Float64} - @test Sf.V isa FixedVariablesSet{<:AbstractVariable, Float64} + @test Sf.V isa FixedVariablesSet{<:AbstractVariable,Float64} @test rem(x + y, ideal(Sf.V)) == 3 S = @set x == 1 && x ≤ y && 2 == y && 1 == x @test S isa BasicSemialgebraicSet{Int} - @test S.V isa FixedVariablesSet{<:AbstractVariable, Int} + @test S.V isa FixedVariablesSet{<:AbstractVariable,Int} @test rem(x + y, ideal(S.V)) == 3 Sf = similar(S, Float64) @test Sf isa BasicSemialgebraicSet{Float64} - @test Sf.V isa FixedVariablesSet{<:AbstractVariable, Float64} + @test Sf.V isa FixedVariablesSet{<:AbstractVariable,Float64} @test rem(x + y, ideal(Sf.V)) == 3 @testset "Empty" begin - for S in (@set(x == 1 && x ≤ y && 2 == y && 2 == x), - @set(x == 1 && x ≤ y && 2 == x), - @set(x == 1 && x ≤ y && 2 == x && y == 3), - @set(x ≤ y && y == 3 && x == 1 && y == 1), - @set(x ≤ y && y == 3 && x == 1 && y == 3 && x == 2)) + for S in ( + @set(x == 1 && x ≤ y && 2 == y && 2 == x), + @set(x == 1 && x ≤ y && 2 == x), + @set(x == 1 && x ≤ y && 2 == x && y == 3), + @set(x ≤ y && y == 3 && x == 1 && y == 1), + @set(x ≤ y && y == 3 && x == 1 && y == 3 && x == 2) + ) @test S isa BasicSemialgebraicSet{Int} - @test S.V isa FixedVariablesSet{<:AbstractVariable, Int} + @test S.V isa FixedVariablesSet{<:AbstractVariable,Int} @test isempty(S.V) @test iszero(length(S.V)) @test isempty(collect(S.V)) @@ -119,23 +134,27 @@ Algebraic Set defined by no equality end end @testset "Basic mixed fixed variables and equalities" begin - for S in [@set(x == 1 && x ≤ y && 2 + y == x), - @set(x == 1 && 2 + y == x && x ≤ y), - @set(x ≤ y && 2 + y == x && x == 1), - @set(x ≤ y && x == 1 && 2 + y == x), - @set(2 + y == x && x ≤ y && x == 1), - @set(2 + y == x && x ≤ y && x == 1)] + for S in [ + @set(x == 1 && x ≤ y && 2 + y == x), + @set(x == 1 && 2 + y == x && x ≤ y), + @set(x ≤ y && 2 + y == x && x == 1), + @set(x ≤ y && x == 1 && 2 + y == x), + @set(2 + y == x && x ≤ y && x == 1), + @set(2 + y == x && x ≤ y && x == 1) + ] @test S isa BasicSemialgebraicSet{Rational{BigInt}} @test S.V isa AlgebraicSet{Rational{BigInt}} end solver = DummySolver() - for S in [@set(x == 1 && x ≤ y && 2 + y == x, solver), - @set(x == 1 && 2 + y == x && x ≤ y, solver), - @set(x ≤ y && 2 + y == x && x == 1, solver), - @set(x ≤ y && x == 1 && 2 + y == x, solver), - @set(2 + y == x && x ≤ y && x == 1, solver), - @set(2 + y == x && x ≤ y && x == 1, solver)] + for S in [ + @set(x == 1 && x ≤ y && 2 + y == x, solver), + @set(x == 1 && 2 + y == x && x ≤ y, solver), + @set(x ≤ y && 2 + y == x && x == 1, solver), + @set(x ≤ y && x == 1 && 2 + y == x, solver), + @set(2 + y == x && x ≤ y && x == 1, solver), + @set(2 + y == x && x ≤ y && x == 1, solver) + ] @test S isa BasicSemialgebraicSet{Rational{BigInt}} @test S.V isa AlgebraicSet{Rational{BigInt}} @test S.V.solver === solver diff --git a/test/solve.jl b/test/solve.jl index 4d80fd0..ba1821b 100644 --- a/test/solve.jl +++ b/test/solve.jl @@ -13,23 +13,27 @@ using LinearAlgebra # for I A = [zeros(10, 1) Matrix(I, 10, 10); zeros(1, 10) 0.5] A[10, 11] = 0 A[10, 1] = η - @test sort.(SemialgebraicSets.clusterordschur(A, sqrt(eps(Float64)))[2]) == [[1, 2, 3, 4, 5, 6, 7, 8, 9, 10], [11]] + @test sort.(SemialgebraicSets.clusterordschur(A, sqrt(eps(Float64)))[2]) == + [[1, 2, 3, 4, 5, 6, 7, 8, 9, 10], [11]] end @testset "Example 4.1 MD95" begin - A = [ 0 0 0 1 0 0; - 0 0 0 0 1 0; - 0 0 0 0 0 1; - -1 0 1 -1 0 0; - -1 0 1 -2 -1 0; - -1 0 1 -2 -2 -1] - @test sort.(SemialgebraicSets.clusterordschur(A, sqrt(eps(Float64)))[2]) == [[2], [1, 5, 6]] + A = [ + 0 0 0 1 0 0 + 0 0 0 0 1 0 + 0 0 0 0 0 1 + -1 0 1 -1 0 0 + -1 0 1 -2 -1 0 + -1 0 1 -2 -2 -1 + ] + @test sort.(SemialgebraicSets.clusterordschur(A, sqrt(eps(Float64)))[2]) == + [[2], [1, 5, 6]] end -function testelements(X, Y; atol=Base.rtoldefault(Float64), kwargs...) +function testelements(X, Y; atol = Base.rtoldefault(Float64), kwargs...) @test length(X) == length(Y) for y in Y - @test any(x -> isapprox(x, y; atol=atol, kwargs...), X) + @test any(x -> isapprox(x, y; atol = atol, kwargs...), X) end end function testelementstypes(X, Y; kwargs...) @@ -49,7 +53,10 @@ end # We use a fixed RNG in the tests to decrease nondeterminism. There is still nondeterminism in LAPACK though using Random -solver = ReorderedSchurMultiplicationMatricesSolver(sqrt(eps(Float64)), MersenneTwister(0)) +solver = ReorderedSchurMultiplicationMatricesSolver( + sqrt(eps(Float64)), + MersenneTwister(0), +) @testset "Zero-dimensional ideal" begin Mod.@polyvar x y z @@ -63,7 +70,7 @@ solver = ReorderedSchurMultiplicationMatricesSolver(sqrt(eps(Float64)), Mersenne testelementstypes(V, [[0]]) V = @set y == x^2 && z == x^3 solver @test !iszerodimensional(V) - V = @set x^3 == 2x*y && x^2*y == 2y^2 + x solver + V = @set x^3 == 2x * y && x^2 * y == 2y^2 + x solver @test iszerodimensional(V) testelementstypes(V, [[0, 0]]) V = @set x == 1 solver @@ -75,7 +82,7 @@ solver = ReorderedSchurMultiplicationMatricesSolver(sqrt(eps(Float64)), Mersenne V = @set x == 4 && y^2 == x solver @test iszerodimensional(V) testelementstypes(V, [[4, 2], [4, -2]]) - V = @set x^2 + x == 6 && y == x+1 solver + V = @set x^2 + x == 6 && y == x + 1 solver @test iszerodimensional(V) testelements(V, [[2, 3], [-3, -2]]) V = @set x^2 + x == 6 && y^2 == x solver @@ -106,17 +113,33 @@ end @testset "Example 5.1 of CGT97" begin ɛ = 1e-4 - Iɛ = [1 - ɛ 0 - 0 1 + ɛ] - J = [0 1 - 1 0] + Iɛ = [ + 1-ɛ 0 + 0 1+ɛ + ] + J = [ + 0 1 + 1 0 + ] Z = zeros(2, 2) - A = [Iɛ Z - Z J] - B = [J Z - Z Iɛ] + A = [ + Iɛ Z + Z J + ] + B = [ + J Z + Z Iɛ + ] α = 0.219 - testelements(SemialgebraicSets._solvemultiplicationmatrices([A, B], [α, 1-α], ReorderedSchurMultiplicationMatricesSolver{Float64}()), [[1.0, -1.0], [1.0, 1.0], [-1.0, 1.0]]; rtol=1e-7) + testelements( + SemialgebraicSets._solvemultiplicationmatrices( + [A, B], + [α, 1 - α], + ReorderedSchurMultiplicationMatricesSolver{Float64}(), + ), + [[1.0, -1.0], [1.0, 1.0], [-1.0, 1.0]]; + rtol = 1e-7, + ) end @testset "Example 4.3 of MD95" begin @@ -125,17 +148,50 @@ end # This test is tricky because in the schur decomposition, the 4 last eigenvalues are e.g. 3.4e-7, -1.7e-7+3e-7im, -1.7e-7-3e-7im, -6e-16 # the second and third do not seem that close but when the three first are averaged it is very close to zero. @test iszerodimensional(V) - testelementstypes(V, [[0.66209555, 0.935259169], [0.66209555, -0.935259169], [0.0516329456, -0.025825086], [0.0516329456, 0.025825086], [0, 0]]) + testelementstypes( + V, + [ + [0.66209555, 0.935259169], + [0.66209555, -0.935259169], + [0.0516329456, -0.025825086], + [0.0516329456, 0.025825086], + [0, 0], + ], + ) end @testset "Example 5.2 of CGT97" begin Mod.@polyvar x y z - V = @set x^2 + y^2 == 1 && x^3 + (2 + z)*x*y + y^3 == 1 && z^2 == 2 solver + V = + @set x^2 + y^2 == 1 && x^3 + (2 + z) * x * y + y^3 == 1 && z^2 == 2 solver @test iszerodimensional(V) iszd, B = monomialbasis(V.I) @test iszd - @test B == [1, z, y, x, y*z, y^2, x*z, x*y, y^2*z, y^3, x*y*z, y^3*z] - testelements(V, [[0, 1, √2], [0, 1, -√2], [1, 0, -√2], [1, 0, √2], [-√2/2, -√2/2, √2], [√2/2, √2/2, -√2]]) + @test B == [ + 1, + z, + y, + x, + y * z, + y^2, + x * z, + x * y, + y^2 * z, + y^3, + x * y * z, + y^3 * z, + ] + testelements( + V, + [ + [0, 1, √2], + [0, 1, -√2], + [1, 0, -√2], + [1, 0, √2], + [-√2 / 2, -√2 / 2, √2], + [√2 / 2, √2 / 2, -√2], + ], + ) end #@testset "Example 4.4 of MD95 and 5.3 of CGT97" begin diff --git a/test/variety.jl b/test/variety.jl index de6c1ae..35752f7 100644 --- a/test/variety.jl +++ b/test/variety.jl @@ -6,8 +6,11 @@ const MP = MultivariatePolynomials using SemialgebraicSets struct DummyAlgo <: SemialgebraicSets.AbstractGröbnerBasisAlgorithm end -function SemialgebraicSets.gröbnerbasis!(::AbstractVector{<:MP.APL}, ::DummyAlgo) - error("Dummy algo") +function SemialgebraicSets.gröbnerbasis!( + ::AbstractVector{<:MP.APL}, + ::DummyAlgo, +) + return error("Dummy algo") end struct DummySolver <: SemialgebraicSets.AbstractAlgebraicSolver end @@ -15,7 +18,7 @@ function SemialgebraicSets.solvealgebraicequations( ::SemialgebraicSets.AlgebraicSet, ::DummySolver, ) - error("Dummy solver") + return error("Dummy solver") end function algo_and_solver_test(func)