diff --git a/src/DynamicPolynomials.jl b/src/DynamicPolynomials.jl index 68b326d..895172e 100644 --- a/src/DynamicPolynomials.jl +++ b/src/DynamicPolynomials.jl @@ -13,8 +13,12 @@ include("var.jl") #const NonCommutativeVariable{O,M} = Variable{NonCommutative{O},M} include("mono.jl") const DMonomialLike{V,M} = Union{Monomial{V,M},Variable{V,M}} -MA.mutability(::Type{<:Monomial}) = MA.IsMutable() +MA.mutability(::Type{<:Monomial{<:Commutative}}) = MA.IsMutable() +MA.mutability(::Type{<:Monomial{<:NonCommutative}}) = MA.IsNotMutable() const _Term{V,M,T} = MP.Term{T,Monomial{V,M}} +function __add_variables!(t::_Term, allvars, map) + return __add_variables!(MP.monomial(t), allvars, map) +end include("monomial_vector.jl") include("poly.jl") MA.mutability(::Type{<:Polynomial}) = MA.IsMutable() @@ -28,7 +32,7 @@ function MP.variable_union_type( end MP.constant_monomial(::Type{<:PolyType{V,M}}) where {V,M} = Monomial{V,M}() function MP.constant_monomial(p::PolyType) - return Monomial(MP.variables(p), zeros(Int, nvariables(p))) + return Monomial(copy(MP.variables(p)), zeros(Int, nvariables(p))) end MP.monomial_type(::Type{<:PolyType{V,M}}) where {V,M} = Monomial{V,M} MP.monomial_type(::PolyType{V,M}) where {V,M} = Monomial{V,M} diff --git a/src/mono.jl b/src/mono.jl index 2cf7183..b3d97c6 100644 --- a/src/mono.jl +++ b/src/mono.jl @@ -100,3 +100,16 @@ MP.monomial(m::Monomial) = m # end # i > length(m1.z) #end + +function __add_variables!( + mono::Monomial{V,M}, + allvars::Vector{Variable{V,M}}, + map, +) where {V,M} + Future.copy!(mono.vars, allvars) + tmp = copy(mono.z) + resize!(mono.z, length(allvars)) + fill!(mono.z, 0) + mono.z[map] = tmp + return +end diff --git a/src/monomial_vector.jl b/src/monomial_vector.jl index a01e993..c636d57 100644 --- a/src/monomial_vector.jl +++ b/src/monomial_vector.jl @@ -355,7 +355,7 @@ function MP.merge_monomial_vectors(ms::Vector{MonomialVector{V,M}}) where {V,M} return MonomialVector{V,M}(buildZvarsvec(Variable{V,M}, X)...) end -function _add_variables!( +function __add_variables!( monos::MonomialVector{V,M}, allvars::Vector{Variable{V,M}}, map, diff --git a/src/operators.jl b/src/operators.jl index 35f1470..1a7d0cb 100644 --- a/src/operators.jl +++ b/src/operators.jl @@ -125,11 +125,7 @@ function MA.operate!( q::Polynomial{V}, ) where {V<:Commutative} if MP.variables(p) != MP.variables(q) - varsvec = [MP.variables(p), MP.variables(q)] - allvars, maps = mergevars(varsvec) - if length(allvars) != length(MP.variables(p)) - _add_variables!(p.x, allvars, maps[1]) - end + allvars, maps = ___add_variables!(p, q) if length(allvars) == length(MP.variables(q)) rhs = q else @@ -139,7 +135,7 @@ function MA.operate!( # should be better of `q` has less terms and then the same term is compared # many times. rhs = Polynomial(q.a, copy(q.x)) - _add_variables!(rhs.x, allvars, maps[2]) + __add_variables!(rhs, allvars, maps[2]) end return MA.operate!(op, p, rhs) end diff --git a/src/poly.jl b/src/poly.jl index 7db0896..0c8f727 100644 --- a/src/poly.jl +++ b/src/poly.jl @@ -402,3 +402,7 @@ function MP.map_exponents!(f::Function, p::Polynomial, m::DMonomialLike) MP.map_exponents!(f, p.x, m) return p end + +function __add_variables!(p::Polynomial, allvars, map) + return __add_variables!(p.x, allvars, map) +end diff --git a/src/subs.jl b/src/subs.jl index e52f0f6..216ca00 100644 --- a/src/subs.jl +++ b/src/subs.jl @@ -94,15 +94,38 @@ function subsmap(st, vars, s::Tuple{MP.VectorMultiSubstitution}) end end +_add_variables!(α, β) = α +_add_variables!(p::PolyType, α) = p +_add_variables!(α, p::PolyType) = MP.operate!!(*, α, one(p)) +function _add_variables!(x::Variable, p::PolyType) + return MP.operate!!(*, x, one(p)) +end +function ___add_variables!(p, q) + varsvec = [MP.variables(p), MP.variables(q)] + allvars, maps = mergevars(varsvec) + if length(allvars) != length(MP.variables(p)) + __add_variables!(p, allvars, maps[1]) + end + return allvars, maps +end +function _add_variables!(p::PolyType, q::PolyType) + ___add_variables!(p, q) + return p +end + function monoeval(z::Vector{Int}, vals::AbstractVector) @assert length(z) == length(vals) if isempty(z) return one(eltype(vals))^1 end + # `Base.power_by_squaring` does a `copy` if `z[1]` is `1` + # which is redirected to `MA.mutable_copy` val = vals[1]^z[1] for i in 2:length(vals) - if z[i] > 0 - val *= vals[i]^z[i] + if iszero(z[i]) + val = _add_variables!(val, vals[i]) + else + val = MA.operate!!(*, val, vals[i]^z[i]) end end return val @@ -122,7 +145,7 @@ function _subs( if iszero(p) zero(Base.promote_op(*, S, T)) else - sum(i -> p.a[i] * monoeval(p.x.Z[i], vals), 1:length(p)) + sum(i -> p.a[i] * monoeval(p.x.Z[i], vals), eachindex(p.a)) end end function _subs( @@ -135,7 +158,7 @@ function _subs( Polynomial{V,M,Tout}, mergevars_of(Variable{V,M}, vals)[1], ) - for i in 1:length(p.a) + for i in eachindex(p.a) MA.operate!(+, q, p.a[i] * monoeval(p.x.Z[i], vals)) end return q diff --git a/src/term.jl b/src/term.jl deleted file mode 100644 index 22de195..0000000 --- a/src/term.jl +++ /dev/null @@ -1,77 +0,0 @@ -export Term - -struct Term{C,T} <: AbstractTerm{T} - α::T - x::Monomial{C} -end -MP.term(α, mono::Monomial) = Term(α, mono) - -iscomm(::Type{Term{C,T}}) where {C,T} = C - -Base.convert(::Type{Term{C,T}}, t::Term{C,T}) where {C,T} = t -function Base.convert(::Type{Term{C,T}}, t::Term{C}) where {C,T} - return Term{C,T}(T(t.α), t.x) -end -Term(t::Term) = t - -function Base.convert(::Type{Term{C,T}}, x::Monomial{C}) where {C,T} - return Term{C,T}(one(T), x) -end -Term{C}(x::Monomial{C}) where {C} = convert(Term{C,Int}, x) -Term(x::Monomial{C}) where {C} = Term{C}(x) - -function Base.convert(::Type{Term{C,T}}, x::Variable{C}) where {C,T} - return convert(Term{C,T}, convert(Monomial{C}, x)) -end -Term{C}(x::Variable{C}) where {C} = Term{C}(convert(Monomial{C}, x)) -Term(x::Variable{C}) where {C} = Term{C}(x) - -function MP.convert_constant(::Type{Term{C,T}}, α) where {C,T} - return Term{C}(convert(T, α)) -end -Term{C}(α::T) where {C,T} = Term{C,T}(α, Monomial{C}()) - -Base.broadcastable(t::Term) = Ref(t) -#(::Type{TermContainer{C}}){C}(x::Variable{C}) = Term(x) -#(::Type{TermContainer{C}}){C}(x::Monomial{C}) = Term(x) -#(::Type{TermContainer{C}}){C}(t::TermContainer{C}) = t -#TermContainer(x::Variable) = Term(x) -#TermContainer(x::Monomial) = Term(x) -#TermContainer(t::TermContainer) = t - -#Base.convert{C, T}(::Type{TermContainer{C, T}}, x::Union{Monomial{C},Variable{C}}) = Term{C, T}(x) -#Base.convert{C, T}(::Type{TermContainer{C, T}}, α::T) = Term{C, T}(α, Monomial{C}()) -#Base.convert{C, S, T}(::Type{TermContainer{C, T}}, α::S) = TermContainer{C, T}(T(α)) -#(::Type{TermContainer{C}}){C, T}(α::T) = TermContainer{C, T}(α) - -#Base.convert{C, T}(::Type{TermContainer{C, T}}, t::Term{C}) = Term{C, T}(t) - -function MA.mutable_copy(t::T) where {T<:Term} - return T(MA.copy_if_mutable(t.α), MA.mutable_copy(t.x)) -end -Base.copy(t::Term) = MA.mutable_copy(t) - -MP.coefficient(t::Term) = t.α -MP.monomial(t::Term) = t.x -_vars(t) = _vars(t.x) - -function MA.operate_to!( - t::Term, - ::typeof(*), - t1::MP.AbstractTermLike, - t2::MP.AbstractTermLike, -) - MA.operate_to!(t.α, *, coefficient(t1), coefficient(t2)) - MA.operate_to!(t.x, *, monomial(t1), monomial(t2)) - return t -end -function MA.operate!(::typeof(*), t1::Term, t2::MP.AbstractTermLike) - MA.operate!(*, t1.α, coefficient(t2)) - MA.operate!(*, t1.x, monomial(t2)) - return t1 -end -function MA.operate!(::typeof(one), t::Term) - MA.operate!(one, t.α) - MA.operate!(zero, t.x.z) - return t -end