Turn Groebner.jl into an extension

Project.toml

@@ -39,10 +39,12 @@ SymbolicUtils = "d1185830-fcd6-423d-90d6-eec64667417b"
 
 [weakdeps]
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
+Groebner = "0b43b601-686d-58a3-8a1c-6623616c7cd4"
 PreallocationTools = "d236fae5-4411-538c-8e31-a6e3d9e00b46"
 SymPy = "24249f21-da20-56a4-8eb1-6a02cf4ae2e6"
 
 [extensions]
+SymbolicsGroebnerExt = "Groebner"
 SymbolicsPreallocationToolsExt = ["ForwardDiff", "PreallocationTools"]
 SymbolicsSymPyExt = "SymPy"
 
@@ -56,7 +58,7 @@ Distributions = "0.25"
 DocStringExtensions = "0.9"
 DomainSets = "0.6"
 DynamicPolynomials = "0.5"
-Groebner = "0.5"
+Groebner = "0.5, 0.6"
 IfElse = "0.1"
 LaTeXStrings = "1.3"
 LambertW = "0.4.5"
@@ -81,6 +83,7 @@ julia = "1.10"
 [extras]
 BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
+Groebner = "0b43b601-686d-58a3-8a1c-6623616c7cd4"
 Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 PkgBenchmark = "32113eaa-f34f-5b0d-bd6c-c81e245fc73d"
 PreallocationTools = "d236fae5-4411-538c-8e31-a6e3d9e00b46"
@@ -91,4 +94,4 @@ SymPy = "24249f21-da20-56a4-8eb1-6a02cf4ae2e6"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Test", "SafeTestsets", "Pkg", "PkgBenchmark", "PreallocationTools", "ForwardDiff", "BenchmarkTools", "ReferenceTests", "SymPy", "Random"]
+test = ["Test", "SafeTestsets", "Pkg", "PkgBenchmark", "PreallocationTools", "ForwardDiff", "Groebner", "BenchmarkTools", "ReferenceTests", "SymPy", "Random"]

ext/SymbolicsGroebnerExt.jl

@@ -0,0 +1,73 @@
+module SymbolicsGroebnerExt
+
+using Groebner
+
+if isdefined(Base, :get_extension)
+    using Symbolics
+    using Symbolics: Num, symtype
+else
+    using ..Symbolics
+    using ..Symbolics: Num, symtype
+end
+
+"""
+    groebner_basis(polynomials; kwargs...)
+
+Computes a Groebner basis of the ideal generated by the given `polynomials`
+using Groebner.jl as the backend.
+
+The basis is guaranteed to be unique.
+The algorithm is randomized, and the output is correct with high probability.
+
+If a coefficient in the resulting basis becomes too large to be represented
+exactly, `DomainError` is thrown.
+
+## Optional Arguments
+
+The Groebner.jl backend provides a number of useful keyword arguments, which are
+also available for this function. See `?Groebner.groebner`.
+
+## Example
+
+```jldoctest
+julia> using Symbolics, Groebner
+
+julia> @variables x y;
+
+julia> groebner_basis([x*y^2 + x, x^2*y + y])
+```
+"""
+function Symbolics.groebner_basis(polynomials::Vector{Num}; kwargs...)
+    polynoms, pvar2sym, sym2term = Symbolics.symbol_to_poly(polynomials)
+    basis = Groebner.groebner(polynoms; kwargs...)
+    PolyType = symtype(first(polynomials))
+    Symbolics.poly_to_symbol(basis, pvar2sym, sym2term, PolyType)
+end
+
+"""
+    is_groebner_basis(polynomials; kwargs...)
+
+Checks whether the given `polynomials` forms a Groebner basis using Groebner.jl
+as the backend.
+
+## Optional Arguments
+
+The Groebner.jl backend provides a number of useful keyword arguments, which are
+also available for this function. See `?Groebner.isgroebner`.
+
+## Example
+
+```jldoctest
+julia> using Symbolics, Groebner
+
+julia> @variables x y;
+
+julia> is_groebner_basis([x^2 - y^2, x*y^2 + x, y^3 + y])
+```
+"""
+function Symbolics.is_groebner_basis(polynomials::Vector{Num}; kwargs...)
+    polynoms, _, _ = Symbolics.symbol_to_poly(polynomials)
+    Groebner.isgroebner(polynoms; kwargs...)
+end
+
+end # module

src/Symbolics.jl

@@ -78,6 +78,8 @@ include("equations.jl")
 export Inequality, ≲, ≳
 include("inequality.jl")
 
+using Bijections
+import DynamicPolynomials
 include("utils.jl")
 
 using ConstructionBase
@@ -116,9 +118,8 @@ include("logexpfunctions-lib.jl")
 
 include("linear_algebra.jl")
 
-using Groebner
 include("groebner_basis.jl")
-export groebner_basis
+export groebner_basis, is_groebner_basis
 
 import Libdl
 include("build_function.jl")
@@ -182,6 +183,9 @@ end
 
 @static if !isdefined(Base,:get_extension)
     function __init__()
+        @require Groebner="0b43b601-686d-58a3-8a1c-6623616c7cd4" begin
+            include("../ext/SymbolicsGroebnerExt.jl")
+        end
         @require SymPy="24249f21-da20-56a4-8eb1-6a02cf4ae2e6" begin
             include("../ext/SymbolicsSymPyExt.jl")
         end

src/groebner_basis.jl

@@ -1,104 +1,25 @@
-import DynamicPolynomials
-using Bijections
 
-const DP = DynamicPolynomials
-# extracting underlying polynomial and coefficient type from Polyforms
-underlyingpoly(x::Number) = x
-underlyingpoly(pf::PolyForm) = pf.p
-coefftype(x::Number) = typeof(x)
-coefftype(pf::PolyForm) = DP.coefficienttype(underlyingpoly(pf))
-
-#=
-Converts an array of symbolic polynomials
-into an array of DynamicPolynomials.Polynomials
-=#
-function symbol_to_poly(sympolys::AbstractArray)
-    @assert !isempty(sympolys) "Empty input."
-
-    # standardize input
-    stdsympolys = map(unwrap, sympolys)
-    sort!(stdsympolys, lt=(<ₑ))
-
-    pvar2sym = Bijections.Bijection{Any,Any}()
-    sym2term = Dict{BasicSymbolic,Any}()
-    polyforms = map(f -> PolyForm(f, pvar2sym, sym2term), stdsympolys)
-
-    # Discover common coefficient type
-    commontype = mapreduce(coefftype, promote_type, polyforms, init=Int)
-    @assert commontype <: Union{Integer,Rational} "Only integer and rational coefficients are supported as input."
-
-    # Convert all to DP.Polynomial, so that coefficients are of same type,
-    # and constants are treated as polynomials
-    # We also need this because Groebner does not support abstract types as input
-    polynoms = Vector{DP.Polynomial{DP.Commutative{DP.CreationOrder}, DP.Graded{DP.LexOrder}, commontype}}(undef, length(sympolys))
-    for (i, pf) in enumerate(polyforms)
-        polynoms[i] = underlyingpoly(pf)
-    end
-
-    polynoms, pvar2sym, sym2term
-end
-
-#=
-Converts an array of AbstractPolynomialLike`s into an array of
-symbolic expressions mapping variables w.r.t pvar2sym
-=#
-function poly_to_symbol(polys, pvar2sym, sym2term, ::Type{T}) where {T}
-    map(f -> PolyForm{T}(f, pvar2sym, sym2term), polys)
-end
+@noinline __throw_absent_groebner_engine() = throw(
+    """Groebner bases engine is required. Execute `using Groebner` to enable this functionality.""")
 
 """
-
-groebner_basis(polynomials)
+    groebner_basis(polynomials)
 
 Computes a Groebner basis of the ideal generated by the given `polynomials`.
-The basis is reduced, thus guaranteed to be unique.
-
-# Example
-
-```jldoctest
-julia> using Symbolics
-
-julia> @variables x y;
-
-julia> groebner_basis([x*y^2 + x, x^2*y + y])
-```
-
-The coefficients in the resulting basis are in the same domain as for input polynomials.
-Hence, if the coefficient becomes too large to be represented exactly, `DomainError` is thrown.
-
-The algorithm is randomized, so the basis will be correct with high probability.
 
+This function requires a Groebner bases backend (such as Groebner.jl) to be loaded.
 """
-function groebner_basis(polynomials)
-    polynoms, pvar2sym, sym2term = symbol_to_poly(polynomials)
-
-    basis = groebner(polynoms)
-
-    # polynomials is nonemtpy
-    T = symtype(first(polynomials))
-    poly_to_symbol(basis, pvar2sym, sym2term, T)
+function groebner_basis(args; kwargs...)
+    __throw_absent_groebner_engine()
 end
 
 """
-
-Do not document this for now.
-
-is_groebner_basis(polynomials)
+    groebner_basis(polynomials)
 
 Checks whether the given `polynomials` forms a Groebner basis.
 
-# Example
-
-```jldoctest
-julia> using Symbolics
-
-julia> @variables x y;
-
-julia> is_groebner_basis([x^2 - y^2, x*y^2 + x, y^3 + y])
-```
-
+This function requires a Groebner bases backend (such as Groebner.jl) to be loaded.
 """
-function is_groebner_basis(polynomials)
-    polynoms, _, _ = symbol_to_poly(polynomials)
-    isgroebner(polynoms)
+function is_groebner_basis(args; kwargs...)
+    __throw_absent_groebner_engine()
 end

src/utils.jl

@@ -331,3 +331,50 @@ function coeff(p, sym=nothing)
         throw(DomainError(p, "Datatype $(typeof(p)) not accepted."))
     end
 end
+
+### Nums <--> Polys
+
+const DP = DynamicPolynomials
+# extracting underlying polynomial and coefficient type from Polyforms
+underlyingpoly(x::Number) = x
+underlyingpoly(pf::PolyForm) = pf.p
+coefftype(x::Number) = typeof(x)
+coefftype(pf::PolyForm) = DP.coefficient_type(underlyingpoly(pf))
+
+#=
+Converts an array of symbolic polynomials
+into an array of DynamicPolynomials.Polynomials
+=#
+function symbol_to_poly(sympolys::AbstractArray)
+    @assert !isempty(sympolys) "Empty input."
+
+    # standardize input
+    stdsympolys = map(unwrap, sympolys)
+    sort!(stdsympolys, lt=(<ₑ))
+
+    pvar2sym = Bijections.Bijection{Any,Any}()
+    sym2term = Dict{BasicSymbolic,Any}()
+    polyforms = map(f -> PolyForm(f, pvar2sym, sym2term), stdsympolys)
+
+    # Discover common coefficient type
+    commontype = mapreduce(coefftype, promote_type, polyforms, init=Int)
+    @assert commontype <: Union{Integer,Rational} "Only integer and rational coefficients are supported as input."
+
+    # Convert all to DP.Polynomial, so that coefficients are of same type,
+    # and constants are treated as polynomials
+    # We also need this because Groebner does not support abstract types as input
+    polynoms = Vector{DP.Polynomial{DP.Commutative{DP.CreationOrder},DP.Graded{DP.LexOrder},commontype}}(undef, length(sympolys))
+    for (i, pf) in enumerate(polyforms)
+        polynoms[i] = underlyingpoly(pf)
+    end
+
+    polynoms, pvar2sym, sym2term
+end
+
+#=
+Converts an array of AbstractPolynomialLike`s into an array of
+symbolic expressions mapping variables w.r.t pvar2sym
+=#
+function poly_to_symbol(polys, pvar2sym, sym2term, ::Type{T}) where {T}
+    map(f -> PolyForm{T}(f, pvar2sym, sym2term), polys)
+end

test/extensions/groebner.jl

@@ -0,0 +1,76 @@
+using Symbolics
+using Symbolics: symbol_to_poly, poly_to_symbol
+using Groebner
+using Test
+
+@variables x y z
+
+syms = [
+    [1], [1, BigInt(2)], [x], [x^2 + 2], [x + (5 // 8)y],
+    [1 - x * y * z], [(x + y + z)^5], [0, BigInt(0)y^2, Rational(1)z^3],
+    [x, sin((44 // 31)y) * z]
+]
+for sym in syms
+    polynoms, pvar2sym, sym2term = Symbolics.symbol_to_poly(sym)
+    sym2 = Symbolics.poly_to_symbol(polynoms, pvar2sym, sym2term, Real)
+    @test isequal(expand.(sym2), expand.(sym))
+end
+
+@test isequal(expand.(groebner_basis([x, y])), [y, x])
+@test isequal(expand.(groebner_basis([y, x])), [y, x])
+@test isequal(expand.(groebner_basis([x])), [x])
+@test isequal(expand.(groebner_basis([x, x^2])), [x])
+@test isequal(expand.(groebner_basis([BigInt(1)x + 2 // 3])), [x + 2 // 3])
+
+@test Symbolics.is_groebner_basis([x, y, z])
+@test Symbolics.is_groebner_basis([x^2 - x, y^2 - y])
+@test !Symbolics.is_groebner_basis([x^2 + y, x * y^3 - 1, y^4 - 1])
+
+@variables x1 x2 x3 x4
+@test isequal(expand.(groebner_basis([x1, x, y])), [y, x1, x])
+
+# input unchanged
+f1 = [x2, x1, x4, x3]
+f2 = [x2, x1, x4, x3]
+groebner_basis(f1)
+@test isequal(expand.(f1), expand.(f2))
+
+@variables x1 x2 x3 x4 x5
+system = [
+    x1 + x2 + x3 + x4 + x5,
+    x1 * x2 + x1 * x3 + x1 * x4 + x1 * x5 + x2 * x3 + x2 * x4 + x2 * x5 + x3 * x4 + x3 * x5 + x4 * x5,
+    x1 * x2 * x3 + x1 * x2 * x4 + x1 * x2 * x5 + x1 * x3 * x4 + x1 * x3 * x5 + x1 * x4 * x5 + x2 * x3 * x4 + x2 * x3 * x5 + x2 * x4 * x5 + x3 * x4 * x5,
+    x1 * x2 * x3 * x4 + x1 * x2 * x3 * x5 + x1 * x2 * x4 * x5 + x1 * x3 * x4 * x5 + x2 * x3 * x4 * x5,
+    x1 * x2 * x3 * x4 * x5 - 1
+]
+truebasis = [
+    x1 + x2 + x3 + x4 + x5,
+    x2^2 + x2 * x3 + x2 * x4 + x2 * x5 + x3^2 + x3 * x4 + x3 * x5 + x4^2 + x4 * x5 + x5^2,
+    x3^3 + x3 * (x4^2) + x3 * (x5^2) + x4^3 + x4 * (x3^2) + x5 * (x4^2) + x4 * (x5^2) + x5^3 + x5 * (x3^2) + x3 * x4 * x5,
+    x4^4 + x4 * (x5^3) + x5^4 + x5 * (x4^3) + (x4^2) * (x5^2),
+    x5^5 - 1
+]
+basis = expand.(groebner_basis(system))
+@test isequal(basis, truebasis)
+
+basis = expand.(groebner_basis(system, linalg=:deterministic))
+@test isequal(basis, truebasis)
+
+N = 45671930739135174346839766056203605080877915151
+system = [
+    x1 + x2 + x3 + x4,
+    x1 * x2 + x1 * x3 + x1 * x4 + x2 * x3 + x2 * x4 + x3 * x4,
+    x1 * x2 * x3 + x1 * x2 * x4 + x1 * x3 * x4 + x2 * x3 * x4,
+    x1 * x2 * x3 * x4 + N
+]
+truebasis = [
+    x1 + x2 + x3 + x4,
+    x2^2 + x2 * x3 + x3^2 + x2 * x4 + x3 * x4 + x4^2,
+    x3^3 + x3^2 * x4 + x3 * x4^2 + x4^3,
+    x4^4 - N
+]
+basis = groebner_basis(system)
+@test isequal(expand.(basis), truebasis)
+
+# Groebner does not yet work with constant ideals
+@test_broken groebner_basis([1])