Changes need for mutable pseudo_rem (#122)

* Changes need for mutable pseudo_rem

* Fix

* Back to master branch

Author: blegat

src/cmult.jl

@@ -28,6 +28,12 @@ function _operate_exponents_to!(output::Vector{Int}, op::F, z1::Vector{Int}, z2:
     @. output = op(z1, z2)
     return
 end
+function _operate_exponents_to!(output::Vector{Vector{Int}}, op::F, z1::Vector{Vector{Int}}, z2::Vector{Int}) where {F<:Function}
+    for i in eachindex(output)
+        _operate_exponents_to!(output[i], op, z1[i], z2)
+    end
+    return
+end
 function _operate_exponents_to!(output::Vector{Int}, op::F, z1::Vector{Int}, z2::Vector{Int}, maps) where {F<:Function}
     I = maps[1]; i = 1; lI = length(I)
     J = maps[2]; j = 1; lJ = length(J)
@@ -133,5 +139,9 @@ function Base.:(*)(x::Monomial{true}, y::MonomialVector{true})
     w, Z = multdivmono(y.vars, x, +, y.Z)
     return MonomialVector{true}(w, Z)
 end
+function MA.operate!(::typeof(*), x::MonomialVector{true}, y::Monomial{true})
+    _multdivmono!(x.Z, x.vars, copy(x.vars), y, +, copy(x.Z))
+    return x
+end
 Base.:(*)(y::MonomialVector{true}, x::Monomial{true}) = x * y
 Base.:(*)(x::Monomial{true}, y::PolyVar{true}) = y * x

src/mult.jl

@@ -35,20 +35,12 @@ end
 include("cmult.jl")
 include("ncmult.jl")
 
-MP.multconstant(α, x::Monomial)   = MP.term(α, MA.mutable_copy(x))
+MP.left_constant_mult(α, x::Monomial)   = MP.term(α, MA.mutable_copy(x))
 
 function zero_with_variables(::Type{Polynomial{C,T}}, vars::Vector{PolyVar{C}}) where{C, T}
     Polynomial(T[], emptymonovec(vars))
 end
 
-function MP._multconstant(α::T, f, p::Polynomial{C,S} ) where {T, C, S}
-    if iszero(α)
-        zero_with_variables(polynomialtype(p, MA.promote_operation(*, T, S)), variables(p))
-    else
-        MP.mapcoefficientsnz(f, p)
-    end
-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)
     Polynomial(MA.mutable_copy(p.a), x*p.x)
@@ -57,10 +49,17 @@ function Base.:(*)(x::DMonomialLike{false}, p::Polynomial)
     # Order may change, e.g. y * (x + y) = y^2 + yx
     Polynomial(monovec(MA.mutable_copy(p.a), [x*m for m in p.x])...)
 end
+
 function Base.:(*)(p::Polynomial, x::DMonomialLike)
     Polynomial(MA.mutable_copy(p.a), p.x*x)
 end
 
+function MA.operate!(::typeof(*), p::Polynomial, t::DMonomialLike)
+    MA.operate!(*, p.x, monomial(t))
+    return p
+end
+
+
 function _term_poly_mult(t::Term{C, S}, p::Polynomial{C, T}, op::Function) where {C, S, T}
     U = MA.promote_operation(op, S, T)
     if iszero(t)
@@ -143,18 +142,3 @@ end
 function MA.operate!(::typeof(*), p::Polynomial{C}, q::Polynomial{C}) where C
     return MA.operate_to!(p, *, p, q)
 end
-
-# Overwrite this method for monomial-like terms because
-# otherwise it would check `iszero(α)` and in that case
-# dismiss of the variable of `p` by performing
-# `operate_to!(zero, output :: Polynomial )` which only
-# respects the variables that are stored already
-function MP._multconstant_to!(output::Polynomial, α, f, p :: DMonomialLike)
-    if iszero(α)
-        MA.operate!(zero, output)
-        Future.copy!(output.x.vars, variables(p))
-        return output
-    else
-        MP.mapcoefficientsnz_to!(output, f, p)
-    end
-end

src/operators.jl

@@ -53,6 +53,14 @@ function plusorminus(p::TermPoly{C, S}, q::TermPoly{C, T}, op::Function) where {
     Polynomial(a, MonomialVector{C}(allvars, Z))
 end
 
+function MA.operate!(::typeof(*), p::Polynomial, t::Term)
+    # In case `coefficient(t)` is a polynomial (e.g. in `gcd` algorithm)
+    # We cannot do `MA.operate!(*, ...)`, we need `MP.right_constant_mult`
+    MA.operate!(MP.right_constant_mult, p, coefficient(t))
+    MA.operate!(*, p, monomial(t))
+    return p
+end
+
 function MA.operate_to!(output::Polynomial{C}, op::Union{typeof(+), typeof(-)},
                                 p::TermPoly{C}, q::TermPoly{C}) where C
     if output === p || output === q

src/poly.jl

@@ -270,7 +270,7 @@ function MA.operate!(::typeof(one), p::Polynomial{C, T}) where {C, T}
         push!(p.x.Z, zeros(Int, length(p.x.vars)))
     else
         resize!(p.a, 1)
-        MA.operate!(one, p.a[1])
+        p.a[1] = MA.operate!!(one, p.a[1])
         resize!(p.x.Z, 1)
         MA.operate!(zero, p.x.Z[1])
     end