diff --git a/src/interface.jl b/src/interface.jl index b834bbd..67e7be1 100644 --- a/src/interface.jl +++ b/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) diff --git a/src/monomial.jl b/src/monomial.jl index 84335c6..c68f6be 100644 --- a/src/monomial.jl +++ b/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) diff --git a/src/mpoly.jl b/src/mpoly.jl index 546eee9..17bb57c 100644 --- a/src/mpoly.jl +++ b/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)) diff --git a/src/packedmonomial.jl b/src/packedmonomial.jl index 3c4090b..bd0719c 100644 --- a/src/packedmonomial.jl +++ b/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