Keep variables for constant zero Polynomials

src/mult.jl

@@ -29,6 +29,19 @@ include("cmult.jl")
 include("ncmult.jl")
 
 MP.multconstant(α, 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
+
 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))
@@ -60,7 +73,7 @@ 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)
-        zero(Polynomial{C, U})
+        zero(Polynomial{C,U})
     else
         n = nterms(p)
         allvars, maps = mergevars([t.x.vars, p.x.vars])
@@ -138,3 +151,18 @@ end
 function MA.mutable_operate!(::typeof(*), p::Polynomial{C}, q::Polynomial{C}) where C
     return MA.mutable_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.mutable_operate!(zero, output)
+        Future.copy!(output.x.vars, variables(p))
+        return output
+    else
+        MP.mapcoefficientsnz_to!(output, f, p)
+    end
+end

test/runtests.jl

@@ -12,6 +12,25 @@ using LinearAlgebra
     @test p(x0 => y0, x1 => y1) == y1 * y0 * y1
 end
 
+@testset "Issue #79 and Issue #80" begin
+    @polyvar x[1:2]
+    p1 = x[1] * 0.0 + x[2] * 0
+    p2 = ( x[1] + x[2] ) * 0.0
+    @test variables(p1) == x 
+    @test variables(p1) == variables(p2)
+
+    @polyvar χ[1:4]
+    poly_array_int = [1 0 1 0; 0 1 0 1]*χ;
+    poly_array_float = Float64.([1 0 1 0; 0 1 0 1])*χ;
+
+    vars_array_int = variables.(poly_array_int);
+    vars_array_float = variables.(poly_array_float)
+    @test isa(vars_array_int, Vector{<:Vector{<:PolyVar}})
+    @test isa(vars_array_float, Vector{<:Vector{<:PolyVar}})
+    @test vars_array_int == vars_array_float
+    @test vars_array_int[1] == vars_array_float[2] == χ
+end
+
 include("mono.jl")
 include("poly.jl")
 include("comp.jl")