diff --git a/Project.toml b/Project.toml
index 158d2d4a..3679242e 100644
--- a/Project.toml
+++ b/Project.toml
@@ -11,9 +11,9 @@ MutableArithmetics = "d8a4904e-b15c-11e9-3269-09a3773c0cb0"
 
 [compat]
 DataStructures = "0.17.7, 0.18"
-DynamicPolynomials = "0.4"
+DynamicPolynomials = "0.4.1"
 MutableArithmetics = "0.3"
-TypedPolynomials = "0.3"
+TypedPolynomials = "0.3.1"
 julia = "1"
 
 [extras]
diff --git a/src/division.jl b/src/division.jl
index 2e0160c2..e23aa55e 100644
--- a/src/division.jl
+++ b/src/division.jl
@@ -1,14 +1,16 @@
 export divides
 
-Base.round(t::AbstractTermLike; args...) = term(round(coefficient(t); args...), monomial(t))
 function Base.round(p::APL; args...)
-    # round(0.1) is zero so we cannot use SortedUniqState
-    polynomial!(round.(terms(p); args...), SortedState())
+    # round(0.1) is zero so we cannot use `mapcoefficientsnz`
+    return mapcoefficients(p) do term
+        round(term; args...)
+    end
 end
 
-Base.div(t::AbstractTermLike, α::Number, args...) = term(div(coefficient(t), α, args...), monomial(t))
 function Base.div(p::APL, α::Number, args...)
-    polynomial!(div.(terms(p), α, args...), SortedState())
+    return mapcoefficients(p) do term
+        div(term, α, args...)
+    end
 end
 
 """
diff --git a/src/operators.jl b/src/operators.jl
index 6b3f5401..01fc071a 100644
--- a/src/operators.jl
+++ b/src/operators.jl
@@ -216,7 +216,7 @@ Base.isapprox(α, p::APL; kwargs...) = isapprox(promote(p, α)...; kwargs...)
 # option.
 Base.:-(m::AbstractMonomialLike) = _term(-1, MA.copy_if_mutable(m))
 Base.:-(t::AbstractTermLike) = _term(MA.operate(-, coefficient(t)), monomial(t))
-Base.:-(p::APL) = polynomial!((-).(terms(p)))
+Base.:-(p::APL) = mapcoefficients(-, p)
 Base.:+(p::Union{APL, RationalPoly}) = p
 Base.:*(p::Union{APL, RationalPoly}) = p
 
@@ -243,39 +243,25 @@ end
 multconstant(α, v::AbstractVariable) = multconstant(α, monomial(v)) # TODO linear term
 multconstant(m::AbstractMonomialLike, α) = multconstant(α, m)
 
-_multconstant(α, f, t::AbstractTermLike) = mapcoefficientsnz(f, t)
-function _multconstant(α::T, f, p::AbstractPolynomial{S}) where {S, T}
-    if iszero(α)
-        zero(polynomialtype(p, MA.promote_operation(*, T, S)))
-    else
-        mapcoefficientsnz(f, p)
-    end
-end
-_multconstant(α, f, p::AbstractPolynomialLike) = _multconstant(α, f, polynomial(p))
+# TODO delete once DynamicPolynomials stops using it
+function _multconstant end
+function _multconstant_to! end
 
-multconstant(α, p::AbstractPolynomialLike) = _multconstant(α, β -> α*β, p)
-multconstant(p::AbstractPolynomialLike, α) = _multconstant(α, β -> β*α, p)
+multconstant(α, p::AbstractPolynomialLike) = mapcoefficients(Base.Fix1(*, α), p)
+multconstant(p::AbstractPolynomialLike, α) = mapcoefficients(Base.Fix2(*, α), p)
 
-# TODO delete once DynamicPolynomials stops using it
-function _multconstant_to!(output, α, f, p)
-    if iszero(α)
-        MA.operate!(zero, output)
-    else
-        mapcoefficientsnz_to!(output, f, p)
-    end
+function MA.operate_to!(output, ::typeof(multconstant), α, p::APL)
+    return mapcoefficients_to!(output, Base.Fix1(*, α), p)
 end
 function MA.operate_to!(output, ::typeof(multconstant), p::APL, α)
-    return _multconstant_to!(output, α, β -> β*α, p)
+    return mapcoefficients_to!(output, Base.Fix2(*, α), p)
 end
-function MA.operate_to!(output, ::typeof(multconstant), α, p::APL)
-    return _multconstant_to!(output, α, β -> α*β, p)
+function MA.operate!(::typeof(multconstant), α, p::APL)
+    return mapcoefficients!(Base.Fix1(*, α), p)
 end
 function MA.operate!(::typeof(multconstant), p::APL, α)
     return mapcoefficients!(Base.Fix2(*, α), p)
 end
-function MA.operate!(::typeof(multconstant), α, p::APL)
-    return mapcoefficients!(Base.Fix1(*, α), p)
-end
 
 MA.operate_to!(output::AbstractMonomial, ::typeof(*), m1::AbstractMonomialLike, m2::AbstractMonomialLike) = mapexponents_to!(output, +, m1, m2)
 MA.operate!(::typeof(*), m1::AbstractMonomial, m2::AbstractMonomialLike) = mapexponents!(+, m1, m2)
diff --git a/src/polynomial.jl b/src/polynomial.jl
index 7402513e..6f3e3e21 100644
--- a/src/polynomial.jl
+++ b/src/polynomial.jl
@@ -409,7 +409,7 @@ See also [`mapcoefficients!`](@ref) and [`mapcoefficients_to!`](@ref).
 Calling `mapcoefficients(α -> mod(3α, 6), 2x*y + 3x + 1)` should return `3x + 3`.
 """
 function mapcoefficients end
-function mapcoefficients(f::Function, p::AbstractPolynomialLike; nonzero = false) # TODO remove
+function mapcoefficients(f::Function, p::AbstractPolynomialLike; nonzero = false) # Not used by either TypedPolynomials or DynamicPolynomials but used by CustomPoly in tests. FIXME Remove in a breaking release
     # Invariant: p has only nonzero coefficient
     # therefore f(α) will be nonzero for every coefficient α of p
     # hence we can use Uniq
diff --git a/test/polynomial.jl b/test/polynomial.jl
index 4e4ed8a3..5cfd7eb0 100644
--- a/test/polynomial.jl
+++ b/test/polynomial.jl
@@ -192,4 +192,21 @@ const MP = MultivariatePolynomials
         q = MA.add!!(q, 2y)
         @test p == 2x + 1
     end
+
+    @testset "mapcoefficients $nz" for nz in [false, true]
+        p = 2x + 1
+        @test mapcoefficients(x -> x / 2, p, nonzero = nz) == 1.0x + 0.5
+        @test mapcoefficients(x -> x / 2, 3x, nonzero = nz) == 1.5x
+        @test p === mapcoefficients!(x -> x + 1, p, nonzero = nz)
+        @test p == 3x + 2
+        q = zero(p)
+        @test q === mapcoefficients_to!(q, x -> 2x, p, nonzero = nz)
+        @test q == 6x + 4
+        @test q === mapcoefficients_to!(q, x -> 2x, 3x, nonzero = nz)
+        @test q == 6x
+        @test q === mapcoefficients_to!(q, x -> 2x, x, nonzero = nz)
+        @test q == 2x
+        @test q === mapcoefficients_to!(q, x -> 2x, x^2, nonzero = nz)
+        @test q == 2x^2
+    end
 end