Groebner.jl integration

src/interface.jl

@@ -19,6 +19,7 @@ Base.:*(c::CoeffType, m::AbstractMonomial) = Term(c, m)
 Base.:*(m::AbstractMonomial, c::CoeffType) = c * m
 Base.:^(y::T, n::Integer) where {T <: AbstractMonomial} = iszero(n) ? T() : Base.power_by_squaring(y, n)
 monomialtype(x::AbstractMonomial) = typeof(x)
+nvariables(x::AbstractMonomial) = error("nvariables not implemented for $(typeof(x))")
 
 ###
 ### AbstractTerm: `coeff` and `monomial`
@@ -41,6 +42,7 @@ Base.one(t::T) where {T<:AbstractTerm} = parameterless_type(T)(one(coeff(t)), on
 Base.isinteger(x::AbstractTerm) = isinteger(coeff(x))
 
 monomialtype(x::AbstractTerm) = monomialtype(monomial(x))
+exps(x::AbstractTerm) = exps(monomial(x))
 
 Base.:*(x::T, y::T) where {T<:AbstractTerm} = T(coeff(x) * coeff(y), monomial(x) * monomial(y))
 # TODO
@@ -124,12 +126,14 @@ Term(x::M) where {M<:AbstractMonomial} = Term(1, x)
 Term{A,B}(x::M) where {A,B,M<:AbstractMonomial} = Term(1, x)
 Term{A,B}(x) where {A,B} = Term(x, B())
 #const EMPTY_TERM = Term[]
+coefftype(::Type{<:Term{T}}) where T = T
 
 ###
 ### AbstractPolynomial: terms, copy
 ###
 abstract type AbstractPolynomial <: Number end
 
+coefftype(p::AbstractPolynomial) = coefftype(eltype(terms(p)))
 Base.promote_rule(p::Type{<:AbstractPolynomial}, ::Type{<:AbstractMonomial}) = p
 Base.promote_rule(p::Type{<:AbstractPolynomial}, ::Type{<:AbstractTerm}) = p
 Base.promote_rule(p::Type{<:AbstractPolynomial}, ::Type{<:CoeffType}) = p
@@ -138,6 +142,7 @@ Base.iszero(x::AbstractPolynomial) = isempty(terms(x))
 Base.isone(x::AbstractPolynomial) = (ts = terms(x); length(ts) == 1 && isone(only(ts)))
 Base.:(==)(x::T, y::T) where {T<:AbstractPolynomial} = x === y || (terms(x) == terms(y))
 monomialtype(p::AbstractPolynomial) = monomialtype(lt(p))
+monomials(p::AbstractPolynomial) = (monomial(t) for t in terms(p))
 
 function Base.show(io::IO, p::AbstractPolynomial)
     ts = terms(p)

src/monomial.jl

@@ -8,6 +8,18 @@ end
 Monomial() = Monomial(EMPTY_IDS)
 degree(x::Monomial) = length(x.ids)
 Base.copy(x::Monomial) = Monomial(copy(x.ids))
+#TODO: optimize
+function nvariables(x::Monomial)
+    nvar = 0
+    lastvar::Int = -1
+    for id in x.ids
+        if id != lastvar
+            lastvar = id
+            nvar += 1
+        end
+    end
+    return nvar
+end
 
 firstid(m::Monomial) = m.ids[1]
 degree(m::Monomial, id) = count(isequal(id), m.ids)

src/mpoly.jl

@@ -9,6 +9,14 @@ MPoly{T}(x::AbstractTerm) where T = MPoly([x])
 MPoly{T}(x::M) where {T,M<:AbstractMonomial} = MPoly(Term(x))
 MPoly{T}(x::CoeffType) where {T} = MPoly(eltype(T)(x))
 terms(x::MPoly) = x.terms
+function nvariables(p::MPoly)
+    if monomialtype(p) <: PackedMonomial
+        nvariables(monomial(lt(p)))
+    else
+        error()
+    end
+end
+coeffs(x::MPoly) = (coeff(t) for t in terms(x))
 
 Base.copy(x::MPoly) = MPoly(map(copy, terms(x)))
 Base.one(::Type{<:MPoly{T}}) where T = MPoly(eltype(T)(1))

src/packedmonomial.jl

@@ -49,6 +49,27 @@ struct PackedMonomial{L,E,K} <: AbstractMonomial
     PackedMonomial{L,E}(bits::NTuple{K,UInt64}) where {L,E,K} = new{L,new_E(Val(E)),K}(bits)
     PackedMonomial{L,E,K}(bits::NTuple{K,UInt64}) where {L,E,K} = new{L,new_E(Val(E)),K}(bits)
 end
+nvariables(x::PackedMonomial{L}) where L = L
+function exps(m::PackedMonomial{L,E,K}) where {L,E,K}
+    tup = m.bits
+    k = 0
+    vpu = var_per_UInt64(Val(E))
+    es = ntuple(Ret(0), L)
+    isone(m) && (return es)
+    for i in eachindex(tup)
+        int = tup[i] << ((i == 1) * (E+1))
+        for j in 1:vpu - (i == 1)
+            exponent = int >> ((E+1)*(vpu-1))
+            if exponent > 0
+                es = Base.setindex(es, Int(exponent), k+1)
+            end
+            k += 1
+            int <<= (E+1)
+            L == k && break
+        end
+    end
+    return es
+end
 
 PackedMonomial{L,E,K}(i::Integer) where {L,E,K} = PackedMonomial{L,E}(i)
 function PackedMonomial{L,OE}(i::Integer) where {L,OE} # x_i