Use promote_operation for subs in isolate_variable (#218)

* Use promote_operation for subs in isolate_variable

* Fix

* Fix degree of constant monomial with no variable

* Missing end

Author: blegat

src/gcd.jl

@@ -175,8 +175,8 @@ function deflation(p::AbstractPolynomialLike)
     end
     @assert all(d -> d >= 0, shift)
     @assert all(d -> d >= 0, defl)
-    s = prod(variables(p).^shift)::monomialtype(p)
-    d = prod(variables(p).^defl)::monomialtype(p)
+    s = prod(variables(p).^shift; init = constantmonomial(p))::monomialtype(p)
+    d = prod(variables(p).^defl; init = constantmonomial(p))::monomialtype(p)
     return s, d
 end
 
@@ -357,7 +357,8 @@ The output can be mutated without affecting `poly` if `mutability` is
 """
 function isolate_variable(poly::APL, var::AbstractVariable, mutability::MA.MutableTrait)
     old_terms = sort!(_vector(terms(_copy(poly, mutability))), by = Base.Fix2(degree, var), rev=true)
-    T = termtype(var, typeof(substitute(Subs(), zero(poly), (var,) => (1,))))
+    U = MA.promote_operation(substitute, Subs, typeof(poly), Pair{typeof(var),Int})
+    T = termtype(var, U)
     new_terms = T[]
     i = firstindex(old_terms)
     while i <= lastindex(old_terms)

src/monomial.jl

@@ -69,7 +69,7 @@ Calling `degree(x^2*y)` should return 3 which is ``2 + 1``.
 Calling `degree(x^2*y, x)` should return 2 and calling `degree(x^2*y, y)` should return 1.
 
 """
-degree(t::AbstractTermLike) = sum(exponents(t))
+degree(t::AbstractTermLike) = sum(exponents(t); init = 0)
 
 degree(t::AbstractTermLike, var::AbstractVariable) = degree(monomial(t), var)
 degree(v::AbstractVariable, var::AbstractVariable) = (v == var ? 1 : 0)

src/promote.jl

@@ -116,3 +116,9 @@ function MA.promote_operation(::typeof(*), PT::Type{<:APL{S}}, QT::Type{<:APL{T}
     U = MA.promote_operation(+, ST, ST)
     return polynomialtype(promote_type(monomialtype(PT), monomialtype(QT)), U)
 end
+function MA.promote_operation(::typeof(*), ::Type{T}, ::Type{P}) where {T, U, P<:APL{U}}
+    return changecoefficienttype(P, MA.promote_operation(*, T, U))
+end
+function MA.promote_operation(::typeof(*), ::Type{P}, ::Type{T}) where {T, U, P<:APL{U}}
+    return changecoefficienttype(P, MA.promote_operation(*, U, T))
+end

src/substitution.jl

@@ -50,6 +50,12 @@ substitute(st::AST, m::AbstractMonomial, s::Substitutions) = powersubstitute(st,
 ## Terms
 substitute(st::AST, t::AbstractTerm, s::Substitutions) = coefficient(t) * substitute(st, monomial(t), s)
 
+function MA.promote_operation(::typeof(substitute), ::Type{Subs}, ::Type{T}, args::Vararg{Type,N}) where {T<:AbstractTerm,N}
+    M = MA.promote_operation(substitute, Subs, monomialtype(T), args...)
+    U = coefficienttype(T)
+    return MA.promote_operation(*, U, M)
+end
+
 ## Polynomials
 _polynomial(α) = α
 _polynomial(p::APL) = polynomial(p)
@@ -59,7 +65,7 @@ function substitute(st::AST, p::AbstractPolynomial, s::Substitutions)
     else
         ts = terms(p)
         r1 = substitute(st, ts[1], s)
-        R = Base.promote_op(+, typeof(r1), typeof(r1))
+        R = MA.promote_operation(+, typeof(r1), typeof(r1))
         result::R = convert(R, r1)
         for i in 2:length(ts)
             result += substitute(st, ts[i], s)
@@ -68,6 +74,11 @@ function substitute(st::AST, p::AbstractPolynomial, s::Substitutions)
     end
 end
 
+function MA.promote_operation(::typeof(substitute), ::Type{Subs}, ::Type{P}, args::Vararg{Type,N}) where {P<:AbstractPolynomial,N}
+    T = MA.promote_operation(substitute, Subs, termtype(P), args...)
+    return MA.promote_operation(+, T, T)
+end
+
 ## Fallbacks
 substitute(st::AST, p::APL, s::Substitutions) = substitute(st, polynomial(p), s)
 substitute(st::AST, q::RationalPoly, s::Substitutions) = substitute(st, q.num, s) / substitute(st, q.den, s)

test/allocations.jl

@@ -30,6 +30,15 @@ function test_polynomial_merge()
     end
 end
 
+
+function test_promotion()
+    @polyvar x y
+    p = x * y + x + y + 1
+    alloc_test(0) do
+        promote_operation(MP.substitute, MP.Subs, typeof(p), Pair{typeof(x),Int})
+    end
+end
+
 function test_isapproxzero()
     @polyvar x y
     p = x * y + x + y + 1
@@ -57,6 +66,7 @@ end
 
 function test_gcd()
     _test_gcd(Int)
+
 end
 
 end

test/promote.jl

@@ -104,4 +104,6 @@ Base.:+(::C, ::C) = C()
 @testset "promote_operation with tricky types" begin
     Mod.@polyvar x
     @test MA.promote_operation(*, polynomialtype(x, A), polynomialtype(x, A)) == polynomialtype(x, C)
+    @test MA.promote_operation(*, polynomialtype(x, A), A) == polynomialtype(x, B)
+    @test MA.promote_operation(*, A, polynomialtype(x, A)) == polynomialtype(x, B)
 end