mapexponents for polynomial (#124)

* mapexponents for polynomial

* Simplify

* Fixes

Author: blegat

src/cmult.jl

@@ -53,22 +53,36 @@ function _operate_exponents_to!(output::Vector{Int}, op::F, z1::Vector{Int}, z2:
     end
     return
 end
+function _operate_exponents_to!(output::Vector{Vector{Int}}, op::F, z1::Vector{Vector{Int}}, z2::Vector{Int}, maps) where {F<:Function}
+    for i in eachindex(output)
+        _operate_exponents_to!(output[i], op, z1[i], z2, maps)
+    end
+    return
+end
 # Not used yet
 #function _operate_exponents!(op::F, Z::Vector{Vector{Int}}, z2::Vector{Int}, args::Vararg{Any,N}) where {F<:Function,N}
 #    return Vector{Int}[_operate_exponents!(op, z, z2, args...) for z in Z]
 #end
+function _resize!(output::Vector, n)
+    resize!(output, n)
+end
+function _resize!(output::Vector{<:Vector}, n)
+    for out in output
+        _resize!(out, n)
+    end
+end
 function _multdivmono!(output, output_variables::Vector{PolyVar{true}},
                       v::Vector{PolyVar{true}}, x::Monomial{true}, op, z)
     if v == x.vars
         if output_variables != v
             resize!(output_variables, length(v))
             copyto!(output_variables, v)
-            resize!(output, length(output_variables))
+            _resize!(output, length(output_variables))
         end
         _operate_exponents_to!(output, op, z, x.z)
     else
         maps = mergevars_to!(output_variables, [v, x.vars])
-        resize!(output, length(output_variables))
+        _resize!(output, length(output_variables))
         _operate_exponents_to!(output, op, z, x.z, maps)
     end
     return
@@ -135,13 +149,23 @@ function MP.mapexponents(f::Function, x::Monomial{true}, y::Monomial{true})
     w, z = multdivmono(x.vars, y, f, x.z)
     return Monomial{true}(w, z)
 end
-function Base.:(*)(x::Monomial{true}, y::MonomialVector{true})
-    w, Z = multdivmono(y.vars, x, +, y.Z)
+function MP.mapexponents(f::Function, x::MonomialVector{true}, y::Monomial{true})
+    w, Z = multdivmono(x.vars, y, f, x.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))
+function MP.mapexponents!(f::Function, x::MonomialVector{true}, y::Monomial{true})
+    _multdivmono!(x.Z, x.vars, copy(x.vars), y, f, copy.(x.Z))
     return x
 end
-Base.:(*)(y::MonomialVector{true}, x::Monomial{true}) = x * y
+function MP.mapexponents(f::Function, x::MonomialVector{true}, y::PolyVar{true})
+    return MP.mapexponents(f, x, MP.monomial(y))
+end
+function MP.mapexponents!(f::Function, x::MonomialVector{true}, y::PolyVar{true})
+    return MP.mapexponents!(f, x, MP.monomial(y))
+end
+function MA.operate!(::typeof(*), x::MonomialVector{true}, y::DMonomialLike{true})
+    return MP.mapexponents!(+, x, y)
+end
+Base.:(*)(y::MonomialVector{true}, x::DMonomialLike{true}) = MP.mapexponents(+, y, x)
+Base.:(*)(x::DMonomialLike{true}, y::MonomialVector{true}) = y * x
 Base.:(*)(x::Monomial{true}, y::PolyVar{true}) = y * x

src/mult.jl

@@ -50,16 +50,6 @@ function Base.:(*)(x::DMonomialLike{false}, p::Polynomial)
     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)
@@ -140,5 +130,11 @@ function MA.operate_to!(p::Polynomial{true, T}, ::typeof(*), q1::MP.AbstractPoly
     end
 end
 function MA.operate!(::typeof(*), p::Polynomial{C}, q::Polynomial{C}) where C
-    return MA.operate_to!(p, *, p, q)
+    if iszero(q)
+        return MA.operate!(zero, p)
+    elseif nterms(q) == 1
+        return MA.operate!(*, p, MP.leadingterm(q))
+    else
+        return MA.operate_to!(p, *, p, q)
+    end
 end

src/poly.jl

@@ -306,3 +306,12 @@ function MP.mapcoefficients_to!(output::Polynomial, f::Function, p::Polynomial;
     end
     return output
 end
+
+function MP.mapexponents(f::Function, p::Polynomial, m::DMonomialLike)
+    return Polynomial(MA.mutable_copy(p.a), MP.mapexponents(f, p.x, m))
+end
+
+function MP.mapexponents!(f::Function, p::Polynomial, m::DMonomialLike)
+    MP.mapexponents!(f, p.x, m)
+    return p
+end