Merge pull request #54 from JuliaAlgebra/bl/copy

Copy internal vectors instead of giving the same instance

Author: blegat

Project.toml

@@ -15,7 +15,7 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [compat]
 DataStructures = "0.17"
-MultivariatePolynomials = "~0.3.3"
+MultivariatePolynomials = "0.3.4"
 MutableArithmetics = "0.1.1"
 Reexport = "0.2"
 julia = "1"

src/cmult.jl

@@ -57,8 +57,7 @@ function multdivmono!(output_variables::Vector{PolyVar{true}},
 end
 function multdivmono(v::Vector{PolyVar{true}}, x::Monomial{true}, op)
     if v == x.vars
-        # /!\ no copy done here for efficiency, do not mess up with vars
-        w = v
+        w = copy(v)
         updatez = z -> op.(z, x.z)
     else
         w, maps = mergevars([v, x.vars])

src/mult.jl

@@ -4,8 +4,7 @@ function multiplyexistingvar(v::Vector{PolyVar{C}}, x::PolyVar{C}, i::Int) where
         newz[i] += 1
         newz
     end
-    # /!\ v not copied for efficiency, do not mess up with vars
-    v, updatez
+    copy(v), updatez
 end
 function insertvar(v::Vector{PolyVar{C}}, x::PolyVar{C}, i::Int) where {C}
     n = length(v)
@@ -29,8 +28,8 @@ end
 include("cmult.jl")
 include("ncmult.jl")
 
-MP.multconstant(α, x::Monomial)   = Term(α, x)
-MP.mapcoefficientsnz(f::Function, p::Polynomial) = Polynomial(f.(p.a), p.x)
+MP.multconstant(α, x::Monomial)   = Term(α, MA.mutable_copy(x))
+MP.mapcoefficientsnz(f::Function, p::Polynomial) = Polynomial(map(f, p.a), MA.mutable_copy(p.x))
 function MP.mapcoefficientsnz_to!(output::Polynomial, f::Function, t::MP.AbstractTermLike)
     MP.mapcoefficientsnz_to!(output, f, polynomial(t))
 end
@@ -48,17 +47,14 @@ end
 
 # I do not want to cast x to TermContainer because that would force the promotion of eltype(q) with Int
 function Base.:(*)(x::DMonomialLike, p::Polynomial)
-    # /!\ No copy of a is done
-    Polynomial(p.a, x*p.x)
+    Polynomial(MA.mutable_copy(p.a), x*p.x)
 end
 function Base.:(*)(x::DMonomialLike{false}, p::Polynomial)
-    # /!\ No copy of a is done
     # Order may change, e.g. y * (x + y) = y^2 + yx
-    Polynomial(monovec(p.a, [x*m for m in p.x])...)
+    Polynomial(monovec(MA.mutable_copy(p.a), [x*m for m in p.x])...)
 end
 function Base.:(*)(p::Polynomial, x::DMonomialLike)
-    # /!\ No copy of a is done
-    Polynomial(p.a, p.x*x)
+    Polynomial(MA.mutable_copy(p.a), p.x*x)
 end
 
 function _term_poly_mult(t::Term{C, S}, p::Polynomial{C, T}, op::Function) where {C, S, T}
@@ -87,7 +83,7 @@ end
 function _mul(::Type{T}, p::Polynomial{true}, q::Polynomial{true}) where T
     samevars = _vars(p) == _vars(q)
     if samevars
-        allvars = _vars(p)
+        allvars = copy(_vars(p))
     else
         allvars, maps = mergevars([_vars(p), _vars(q)])
     end

src/operators.jl

@@ -15,10 +15,6 @@ end
 Base.:(+)(x::DMonomialLike, y::DMonomialLike) = Term(x) + Term(y)
 Base.:(-)(x::DMonomialLike, y::DMonomialLike) = Term(x) - Term(y)
 
-# `+(::UniformScaling)` is not defined...
-_unary(::typeof(+), x) = MA.copy_if_mutable(x)
-_unary(::typeof(-), x) = -x
-
 _getindex(p::Polynomial, i) = p[i]
 _getindex(t::Term, i) = t
 function _plusorminus_to!(a::Vector{U}, Z::Vector{Vector{Int}}, op::Function, p::TermPoly{C}, q::TermPoly{C}, maps, nvars) where {C, U}
@@ -72,6 +68,11 @@ function MA.mutable_operate_to!(output::Polynomial{C}, op::Union{typeof(+), type
     return output
 end
 
+# TODO replace by MA.operate with MA v0.2
+# `+(::UniformScaling)` is not defined...
+_unary(::typeof(+), x) = MA.copy_if_mutable(x)
+_unary(::typeof(-), x) = -x
+
 function MA.mutable_operate!(op::Union{typeof(+), typeof(-)}, p::Polynomial,
                              q::Union{PolyVar, Monomial, Term})
     return MA.mutable_operate!(op, p, polynomial(q))
@@ -98,7 +99,7 @@ function MA.mutable_operate!(op::Union{typeof(+), typeof(-)}, p::Polynomial{true
         return MA.mutable_operate!(op, p, rhs)
     end
     get1(i) = (p.a[i], p.x.Z[i])
-    get2(i) = (MA.scaling_convert(eltype(p.a), MA.copy_if_mutable(q.a[i])), copy(q.x.Z[i]))
+    get2(i) = (MA.scaling_convert(eltype(p.a), _unary(op, q.a[i])), copy(q.x.Z[i]))
     function set(i, t::_NoVarTerm)
         p.a[i] = t[1]
         p.x.Z[i] = t[2]
@@ -133,6 +134,11 @@ Base.:(+)(x::Union{Monomial,PolyVar}, y::TermPoly{C}) where C = Term{C}(x) + y
 Base.:(-)(x::TermPoly{T}, y::DMonomialLike) where T = x - Term{T}(y)
 Base.:(-)(x::DMonomialLike, y::TermPoly{T}) where T = Term{T}(x) - y
 
-Base.:(-)(p::Polynomial) = Polynomial(-p.a, p.x)
+# `MA.operate(-, p)` redirects to `-p` as it assumes that `-p` can be modified
+# through the MA API without modifying `p`. We should either copy `p.x` here
+# or implement a `MA.operate(-, p)` that copies it. We choose the first option.
+Base.:(-)(p::Polynomial) = Polynomial(-p.a, copy(p.x))
+# TODO use that with MA v0.2
+#Base.:(-)(p::Polynomial) = Polynomial(map(α -> MA.operate(-, α), p.a), copy(p.x))
 
 include("mult.jl")