diff --git a/src/mono.jl b/src/mono.jl
index 699ba1b..e25d566 100644
--- a/src/mono.jl
+++ b/src/mono.jl
@@ -80,3 +80,10 @@ MP.monomial(m::Monomial) = m
 #    end
 #    i > length(m1.z)
 #end
+function _add_variables!(mono::Monomial, allvars, map)
+    Future.copy!(mono.vars, allvars)
+    tmp = copy(mono.z)
+    resize!(mono.z, length(allvars))
+    fill!(mono.z, 0)
+    mono.z[map] = tmp
+end
diff --git a/src/operators.jl b/src/operators.jl
index dc2a077..271a9a5 100644
--- a/src/operators.jl
+++ b/src/operators.jl
@@ -68,9 +68,54 @@ function MA.mutable_operate_to!(output::Polynomial{C}, op::Union{typeof(+), type
     return output
 end
 
+_copy_monos(m::Monomial) = copy(m)
+_copy_monos(t::Term) = Term(t.α, copy(t.x))
+_copy_monos(p::Polynomial) = Polynomial(p.a, copy(p.x))
+function _merge_vars!(p, q)
+    if _vars(p) == _vars(q)
+        return q
+    end
+    varsvec = [_vars(p), _vars(q)]
+    allvars, maps = mergevars(varsvec)
+    if length(allvars) != length(_vars(p))
+        _add_variables!(p.x, allvars, maps[1])
+    end
+    if length(allvars) != length(_vars(q))
+        # We could avoid promoting `q` to the same variables
+        # like in `plusorminus` to avoid extra allocation but it then
+        # gives slower comparison. There is a tradeoff and the approach used here
+        # should be better of `q` has less terms and then the same term is compared
+        # many times.
+        q = _copy_monos(q)
+        _add_variables!(q, allvars, maps[2])
+    end
+    return q
+end
 function MA.mutable_operate!(op::Union{typeof(+), typeof(-)}, p::Polynomial,
-                             q::Union{PolyVar, Monomial, Term})
-    return MA.mutable_operate!(op, p, polynomial(q))
+                             v::PolyVar)
+    return MA.mutable_operate!(op, p, monomial(v))
+end
+function MA.mutable_operate!(op::Union{typeof(+), typeof(-)}, p::Polynomial,
+                             t::Union{Monomial, Term})
+    t = _merge_vars!(p, t)
+    z = MP.exponents(t)
+    i = findfirst(eachindex(p.a)) do i
+        return samevars_grlex(p.x.Z[i], z) <= 0
+    end
+    if i === nothing
+        push!(p.a, MA.operate(op, MP.coefficient(t)))
+        push!(p.x.Z, copy(z))
+    else
+        comp = samevars_grlex(p.x.Z[i], z)
+        if iszero(comp)
+            p.a[i] = MA.operate!(op, p.a[i], coefficient(t))
+        else
+            @assert comp < 0
+            insert!(p.a, i, MA.operate(op, MP.coefficient(t)))
+            insert!(p.x.Z, i, copy(z))
+        end
+    end
+    return p
 end
 const _NoVarTerm{T} = Tuple{T, Vector{Int}}
 function MA.mutable_operate!(op::Union{typeof(+), typeof(-)}, p::Polynomial{false},
@@ -80,25 +125,7 @@ end
 # TODO need to check that this also works for non-commutative
 function MA.mutable_operate!(op::Union{typeof(+), typeof(-)}, p::Polynomial{true},
                              q::Polynomial{true})
-    if _vars(p) != _vars(q)
-        varsvec = [_vars(p), _vars(q)]
-        allvars, maps = mergevars(varsvec)
-        if length(allvars) != length(_vars(p))
-            _add_variables!(p.x, allvars, maps[1])
-        end
-        if length(allvars) == length(_vars(q))
-            rhs = q
-        else
-            # We could avoid promoting `q` to the same variables
-            # like in `plusorminus` to avoid extra allocation but it then
-            # gives slower comparison. There is a tradeoff and the approach used here
-            # 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])
-        end
-        return MA.mutable_operate!(op, p, rhs)
-    end
+    q = _merge_vars!(p, q)
     get1(i) = (p.a[i], p.x.Z[i])
     get2(i) = (MA.scaling_convert(eltype(p.a), MA.operate(op, q.a[i])), copy(q.x.Z[i]))
     function set(i, t::_NoVarTerm)
diff --git a/src/poly.jl b/src/poly.jl
index a0eb367..2d406c1 100644
--- a/src/poly.jl
+++ b/src/poly.jl
@@ -262,3 +262,4 @@ function MA.mutable_operate!(::typeof(one), p::Polynomial{C, T}) where {C, T}
     end
     return p
 end
+_add_variables!(p::Polynomial, allvars, map) = _add_variables!(p.x, allvars, map)
diff --git a/src/term.jl b/src/term.jl
index 2a923b2..b7ff5ef 100644
--- a/src/term.jl
+++ b/src/term.jl
@@ -68,3 +68,4 @@ function MA.mutable_operate!(::typeof(one), t::Term)
     MA.mutable_operate!(zero, t.x.z)
     return t
 end
+_add_variables!(t::Term, allvars, map) = _add_variables!(t.x, allvars, map)