-
Notifications
You must be signed in to change notification settings - Fork 156
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
- Loading branch information
Showing
8 changed files
with
223 additions
and
169 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -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 |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -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 |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -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]) |
Oops, something went wrong.