Refactoring (#37)

* Use CommonSolve and rename + format

* Fix format

* monomial_basis -> standard_monomials

* standard monomials docstring

Author: blegat

.github/workflows/format_check.yml

@@ -1,4 +1,4 @@
-ame: format-check
+name: format-check
 on:
   push:
     branches:

Project.toml

@@ -4,12 +4,14 @@ repo = "https://github.com/JuliaAlgebra/SemialgebraicSets.jl.git"
 version = "0.3.0"
 
 [deps]
+CommonSolve = "38540f10-b2f7-11e9-35d8-d573e4eb0ff2"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MultivariatePolynomials = "102ac46a-7ee4-5c85-9060-abc95bfdeaa3"
 MutableArithmetics = "d8a4904e-b15c-11e9-3269-09a3773c0cb0"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 
 [compat]
+CommonSolve = "0.2"
 MultivariatePolynomials = "0.4.2"
 MutableArithmetics = "0.3.1, 1"
 julia = "1.6"

README.md

@@ -33,9 +33,9 @@ using TypedPolynomials
 Once the algebraic set has been created, you can check whether it is zero-dimensional and if it is the case, you can get the finite number of elements of the set simply by iterating over it, or by using `collect` to transform it into an array containing the solutions.
 ```julia
 V = @set y == x^2 && z == x^3
-iszerodimensional(V) # should return false
+is_zero_dimensional(V) # should return false
 V = @set x^2 + x == 6 && y == x+1
-iszerodimensional(V) # should return true
+is_zero_dimensional(V) # should return true
 collect(V) # should return [[2, 3], [-3, -2]]
 ```
 The code sample above solves the system of algbraic equations by first

src/SemialgebraicSets.jl

@@ -10,18 +10,20 @@ const MP = MultivariatePolynomials
 
 const APL = AbstractPolynomialLike
 
+import CommonSolve: solve
+
 export AbstractSemialgebraicSet,
     AbstractBasicSemialgebraicSet, AbstractAlgebraicSet
 export FullSpace,
-    AlgebraicSet, BasicSemialgebraicSet, addequality!, addinequality!
+    AlgebraicSet, BasicSemialgebraicSet, add_equality!, add_inequality!
 
 # Semialgebraic set described by polynomials with coefficients in T
 abstract type AbstractSemialgebraicSet end
 
 abstract type AbstractBasicSemialgebraicSet <: AbstractSemialgebraicSet end
 abstract type AbstractAlgebraicSet <: AbstractBasicSemialgebraicSet end
 
-function addinequality!(S::AbstractAlgebraicSet, p)
+function add_inequality!(S::AbstractAlgebraicSet, p)
     throw(ArgumentError("Cannot add inequality to an algebraic set"))
 end
 
@@ -63,4 +65,6 @@ include("basic.jl")
 include("fix.jl")
 include("macro.jl")
 
+include("deprecate.jl")
+
 end # module

src/basic.jl

@@ -1,4 +1,4 @@
-export inequalities, basicsemialgebraicset
+export inequalities, basic_semialgebraic_set
 
 struct BasicSemialgebraicSet{T,PT<:APL{T},AT<:AbstractAlgebraicSet} <:
        AbstractBasicSemialgebraicSet
@@ -26,7 +26,7 @@ function BasicSemialgebraicSet(
     )
 end
 #BasicSemialgebraicSet{T, PT<:APL{T}}(V::AlgebraicSet{T, PT}, p::Vector{PT}) = BasicSemialgebraicSet{T, PT}(V, p)
-function basicsemialgebraicset(V, p)
+function basic_semialgebraic_set(V, p)
     return BasicSemialgebraicSet(V, p)
 end
 
@@ -58,10 +58,10 @@ function MP.variables(S::BasicSemialgebraicSet)
 end
 nequalities(S::BasicSemialgebraicSet) = nequalities(S.V)
 equalities(S::BasicSemialgebraicSet) = equalities(S.V)
-addequality!(S::BasicSemialgebraicSet, p) = addequality!(S.V, p)
+add_equality!(S::BasicSemialgebraicSet, p) = add_equality!(S.V, p)
 ninequalities(S::BasicSemialgebraicSet) = length(S.p)
 inequalities(S::BasicSemialgebraicSet) = S.p
-addinequality!(S::BasicSemialgebraicSet, p) = push!(S.p, p)
+add_inequality!(S::BasicSemialgebraicSet, p) = push!(S.p, p)
 
 function Base.intersect(S::BasicSemialgebraicSet, T::BasicSemialgebraicSet)
     return BasicSemialgebraicSet(S.V ∩ T.V, [S.p; T.p])

src/deprecate.jl

@@ -0,0 +1,22 @@
+@deprecate iszerodimensional is_zero_dimensional
+@deprecate addequality! add_equality!
+@deprecate addinequality! add_inequality!
+@deprecate basicsemialgebraicset basic_semialgebraic_set
+@deprecate reducebasis! reduce_basis!
+@deprecate dummyselection dummy_selection
+@deprecate normalselection normal_selection
+@deprecate groebnerbasis groebner_basis
+@deprecate gröbnerbasis gröbner_basis
+@deprecate defaultgröbnerbasisalgorithm default_gröbner_basis_algorithm
+@deprecate computegröbnerbasis! compute_gröbner_basis!
+@deprecate algebraicset algebraic_set
+@deprecate algebraicsolver algebraic_solver
+@deprecate computeelements! compute_elements!
+@deprecate conditionnumber condition_number
+@deprecate multiplicationmatrices multiplication_matrices
+@deprecate monomialbasis standard_monomials
+@deprecate solvealgebraicequations solve
+@deprecate defaultmultiplicationmatricesalgorithm default_multiplication_matrices_algorithm
+@deprecate defaultalgebraicsolver default_algebraic_solver
+@deprecate projectivealgebraicset projective_algebraic_set
+@deprecate defaultalgebraicsetlibrary default_algebraic_set_library

src/fix.jl

@@ -111,15 +111,15 @@ function Base.intersect(
     return FixedVariablesSet(ideal)
 end
 function Base.intersect(V::AlgebraicSet, W::FixedVariablesSet)
-    return V ∩ algebraicset(W)
+    return V ∩ algebraic_set(W)
 end
 function Base.intersect(V::FixedVariablesSet, W::AlgebraicSet)
     # Need `W` to be first so that the intersection uses its solver.
-    # and not the default one taken in `algebraicset(W)`.
+    # and not the default one taken in `algebraic_set(W)`.
     return W ∩ V
 end
 
-iszerodimensional(::FixedVariablesSet) = true
+is_zero_dimensional(::FixedVariablesSet) = true
 function Base.isempty(V::FixedVariablesSet)
     return generate_nonzero_constant(V.ideal)
 end

src/groebner.jl

@@ -1,5 +1,5 @@
-export spolynomial, gröbnerbasis, groebnerbasis
-export Buchberger, presort!, dummyselection, normalselection
+export spolynomial, gröbner_basis, groebner_basis
+export Buchberger, presort!, dummy_selection, normal_selection
 
 """
     spolynomial(p::AbstractPolynomialLike, q::AbstractPolynomialLike)
@@ -23,7 +23,7 @@ function spolynomial(p::APL, q::APL)
     return MA.operate!!(-, ad, bd)
 end
 
-function reducebasis!(F::AbstractVector{<:APL}; kwargs...)
+function reduce_basis!(F::AbstractVector{<:APL}; kwargs...)
     changed = true
     while changed
         changed = false
@@ -78,7 +78,7 @@ end
 
 # Selection of (i, j) ∈ B
 
-function dummyselection(F, B)
+function dummy_selection(F, B)
     return first(B)
 end
 # "BUCHBERGER (1985) suggests that there will often be some savings if we pick
@@ -91,7 +91,7 @@ end
 #
 # Ideals, Varieties, and Algorithms
 # Cox, Little and O'Shea, Fourth edition p. 116
-function normalselection(F, B)
+function normal_selection(F, B)
     best = first(B)
     m = lcmlm(F, best...)
     for (i, j) in Iterators.drop(B, 1)
@@ -105,25 +105,25 @@ function normalselection(F, B)
 end
 # TODO sugar and double sugar selections
 
-defaultgröbnerbasisalgorithm(p) = Buchberger()
+default_gröbner_basis_algorithm(p) = Buchberger()
 abstract type AbstractGröbnerBasisAlgorithm end
 
 struct Buchberger <: AbstractGröbnerBasisAlgorithm
     ztol::Float64
     pre!::Function
     sel::Function
 end
-#Buchberger() = Buchberger(presort!, dummyselection)
+#Buchberger() = Buchberger(presort!, dummy_selection)
 function Buchberger(ztol = Base.rtoldefault(Float64))
-    return Buchberger(ztol, presort!, normalselection)
+    return Buchberger(ztol, presort!, normal_selection)
 end
 
 # Taken from Theorem 9 of
 # Ideals, Varieties, and Algorithms
 # Cox, Little and O'Shea, Fourth edition
-function gröbnerbasis!(
+function gröbner_basis!(
     F::AbstractVector{<:APL},
-    algo = defaultgröbnerbasisalgorithm(F),
+    algo = default_gröbner_basis_algorithm(F),
 )
     algo.pre!(F)
     B = Set{NTuple{2,Int}}(
@@ -150,12 +150,12 @@ function gröbnerbasis!(
         end
         delete!(B, (i, j))
     end
-    reducebasis!(F; ztol = algo.ztol)
+    reduce_basis!(F; ztol = algo.ztol)
     map!(monic, F, F)
     return F
 end
-function gröbnerbasis(F::Vector{<:APL}, args...)
+function gröbner_basis(F::Vector{<:APL}, args...)
     T = Base.promote_op(rem, eltype(F), typeof(F))
-    return gröbnerbasis!(copyto!(similar(F, T), F), args...)
+    return gröbner_basis!(copyto!(similar(F, T), F), args...)
 end
-const groebnerbasis = gröbnerbasis
+const groebner_basis = gröbner_basis

src/ideal.jl

@@ -1,12 +1,12 @@
-export monomialbasis, ideal
+export standard_monomials, ideal
 
 abstract type AbstractPolynomialIdeal end
 
 struct EmptyPolynomialIdeal <: AbstractPolynomialIdeal end
-function ideal(p::FullSpace, algo = defaultgröbnerbasisalgorithm(p))
+function ideal(p::FullSpace, _ = default_gröbner_basis_algorithm(p))
     return EmptyPolynomialIdeal()
 end
-Base.rem(p::AbstractPolynomialLike, I::EmptyPolynomialIdeal) = p
+Base.rem(p::AbstractPolynomialLike, ::EmptyPolynomialIdeal) = p
 
 promote_for_division(::Type{T}) where {T} = T
 promote_for_division(::Type{T}) where {T<:Integer} = Rational{big(T)}
@@ -18,7 +18,7 @@ end
 mutable struct PolynomialIdeal{T,PT<:APL{T},A<:AbstractGröbnerBasisAlgorithm} <:
                AbstractPolynomialIdeal
     p::Vector{PT}
-    gröbnerbasis::Bool
+    gröbner_basis::Bool
     algo::A
 end
 function PolynomialIdeal{T,PT}(
@@ -35,7 +35,7 @@ function PolynomialIdeal(p::Vector{PT}, algo) where {T,PT<:APL{T}}
     )
 end
 function PolynomialIdeal{T,PT}() where {T,PT<:APL{T}}
-    return PolynomialIdeal(PT[], defaultgröbnerbasisalgorithm(PT))
+    return PolynomialIdeal(PT[], default_gröbner_basis_algorithm(PT))
 end
 
 function Base.convert(
@@ -48,10 +48,10 @@ 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)
+    return PolynomialIdeal{T,PT,A}(I.p, I.gröbner_basis, I.algo)
 end
 
-ideal(p, algo = defaultgröbnerbasisalgorithm(p)) = PolynomialIdeal(p, algo)
+ideal(p, algo = default_gröbner_basis_algorithm(p)) = PolynomialIdeal(p, algo)
 
 function Base.:+(I::PolynomialIdeal, J::PolynomialIdeal)
     return PolynomialIdeal([I.p; J.p], I.algo)
@@ -60,20 +60,39 @@ end
 MP.variables(I::PolynomialIdeal) = variables(I.p)
 
 function Base.rem(p::AbstractPolynomialLike, I::PolynomialIdeal)
-    computegröbnerbasis!(I)
+    compute_gröbner_basis!(I)
     return rem(p, I.p)
 end
 
-function computegröbnerbasis!(I::PolynomialIdeal)
-    if !I.gröbnerbasis
-        gröbnerbasis!(I.p, I.algo)
-        I.gröbnerbasis = true
+function compute_gröbner_basis!(I::PolynomialIdeal)
+    if !I.gröbner_basis
+        gröbner_basis!(I.p, I.algo)
+        I.gröbner_basis = true
     end
 end
-function monomialbasis(I, vars = variables(I))
-    computegröbnerbasis!(I)
+
+"""
+    standard_monomials(I::AbstractPolynomialIdeal, vars = variables(I))
+
+Return the vector of *standard monomials* in the variables `var` of the ideal
+`I`. A monomial is a *standard monomial* if it is not the leading monomial of
+any polynomial of the ideal. If there are infinitely such monomials, `nothing`
+is returned.
+
+See [CLO5, p. 38], where it is also called *basis monomial*, for more
+information.
+
+[CLO05] Cox, A. David & Little, John & O'Shea, Donal
+*Using Algebraic Geometry*.
+Graduate Texts in Mathematics, **2005**.
+https://doi.org/10.1007/b138611
+"""
+function standard_monomials end
+
+function standard_monomials(I, vars = variables(I))
+    compute_gröbner_basis!(I)
     if isempty(I.p)
-        return false, monomial_type(eltype(I.p))[]
+        return nothing # Nonzero dimensional
     end
     mv = monomial_vector(leading_monomial.(I.p))
     # monomial_vector makes sure all monomials have the same variables
@@ -90,9 +109,8 @@ function monomialbasis(I, vars = variables(I))
             end
         end
     end
-    if any(lv .< 0)
-        return false, monomial_type(eltype(I.p))[]
+    if any(Base.Fix2(<, 0), lv)
+        return nothing # Nonzero dimensional
     end
-    return true,
-    monomials(vars, 0:(sum(lv)-1), m -> !any(map(m2 -> divides(m2, m), mv)))
+    return monomials(vars, 0:(sum(lv)-1), m -> !any(Base.Fix2(divides, m), mv))
 end

src/macro.jl

@@ -15,9 +15,9 @@ function equality(lhs, rhs)
 end
 function Base.intersect(eq::PolynomialEquality; lib_or_solver = nothing)
     if lib_or_solver === nothing
-        return algebraicset([eq.p])
+        return algebraic_set([eq.p])
     else
-        return algebraicset([eq.p], lib_or_solver)
+        return algebraic_set([eq.p], lib_or_solver)
     end
 end