Fix variable with subs of zero power (#135)

* Fix variable with subs of zero power

* Fix

* Fix

Author: blegat

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}

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

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,

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

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

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