diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index 47be6240..e5935259 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -41,7 +41,7 @@ jobs:
         run: |
           using Pkg
           Pkg.add([
-              PackageSpec(name="DynamicPolynomials", rev="master"),
+              PackageSpec(name="DynamicPolynomials", rev="release-0.7"),
               PackageSpec(name="TypedPolynomials", rev="master"),
               PackageSpec(name="StarAlgebras", rev="bl/term"),
           ])
diff --git a/AGENTS.md b/AGENTS.md
new file mode 100644
index 00000000..0e6c6b3b
--- /dev/null
+++ b/AGENTS.md
@@ -0,0 +1,57 @@
+
+This repository as well as DynamicPolynomials have accumulated a lot of technical depth due to bad design decisions and it is time to get rid of it.
+The bad design decisions were
+1. Try to support `*(::MP.AbstractPolynomialLike, ::Any)` this creates a lot of invalidation and promotion rules are a mess. We should only support `*(::MP.AbstractPolynomialLike{T}, ::T)` and have `*(::MP.AbstractPolynomialLike, ::Number)` that convert the number to `T`. Same for other operators like `+`, `-` etc..
+2. Implement `+(p::MP.AbstractPolynomialLike, q::MP.AbstractPolynomialLike)` where `p` and `q` have different variables. Not simply by first promoting them to the same set of variables and then use a simple implementation of the sum but instead to have some complicated implementation dealing with monomials over different variables. We should just directly promote with `StarAlgebras.promote_bases` now and only implement `+` over polynomials of the same algebra
+3. Having different implementation of polynomials, `MP.Polynomial` using a list of terms and `DynamicPolynomials.Polynomial` having a separate list of coefficients and list of monomials. We should just use StarAlgebras.AlgebraElement now.
+
+So the plan is:
+1. Move `MP.Term` to `StarAlgebras.Term`, so `MP.AbstractTermLike` won't be an abstract type anymore but a union of `AbstractMonomialLike` and `SA.Term`
+2. Remove `MP.Polynomial`, and replace it by the implementation of polynomials done in MultivariateBases using StarAlgebras. Doing so, we can just merge MultivariateBases into MultivariatePolynomials, `MP.AbstractPolynomialLike` will then be a union of that StarAlgebras.AlgebraElement and `AbstractTermLike`.
+
+Terminology:
+
+MOI: MathOptInterface
+MA: MutableArithmetics
+DP: DynamicPolynomials
+TP: TypedPolynomials
+MP: MultivariatePolynomials
+SA: StarAlgebras
+MB: MultivariateBases
+
+All packages should be in ~/.julia/dev
+
+MA has a large test suite that is used by JuMP, MOI, MP and therefore also DP and TP. It is in MP/test. It's nice but it takes a lot of time to run. First makes sure everything is done before running it for the final touches.
+This refactor is supposed to get rid of a lot of code and technical depth, don't start writing a lot of code before first asking me.
+
+The promotion rules are a mess. Because we had to support promotion with Any, this messed up badly with Julia's internals.
+In Julia, normally you just need to implement promote_rule(::A, ::B), not promote_rule(::B, ::A). But because I implemented promotion with ::Any, I had to do both and it created such a mess!!
+It was also a rabbithole where I then needed also APL{<:Any} and things like that, quite complicated...
+I want things to be easier now. Before, I belived if x is a variable, promote_type(typeof(x), typeof("x")) would be something like Term{String}. We don't want that anymore, let it just be Any[].
+You see the vibe ? Let's start simple and then be really careful for test we used to have whether we really want this test to still pass or we just want to simplify.
+
+About backward compat, we are going to make a breaking release so I prefer being more breaking and having a simpler code than the opposite.
+
+A lot of code in MP are like sum(t::Vector{<:AbstractTerm}) = dot(coefficient.(t), monomial.(t)) and then dot(c::AbstractVector, m::AbstractVector{<:AbstractMonomialLike}) = sum(c .* m).
+These create stackoverflow in case DP or TP forget to implement one of the two! This was because, as explained above, TP used the default polynomials that were a vector of terms and DP
+was a vector of coefficients separated to a vector of monomials. Now it all gets simpler because it will be a separate vector of coefficients and monomials since we'll just be using SA.AlgebraElement!!
+So we can simplify. Also, because SA.AlgebraElement <: MA.AbstractMutable, we already have a fallback for Base.sum so we might just be able to remove these!
+These were written before MutableArithmetics existed. I am maintaining MA, SA, JuMP, MOI, MP, DP, TP and I don't want to maintain duplicates. All these packages define mutable objects.
+MA allows me to have already a lot of code in common of all of them. Then, JuMP, SA and MOI will have different implementations of basically "sum of terms" but that's fine.
+What I don't want is MP, TP and DP to also have their own version, they should just use SA.AlgebraElement!
+
+It has already been started by another AI agent on the branches:
+DP: bl/sa_term
+MP: bl/sa_poly
+SA: bl/term
+
+Just continue on these branches, feel free to simplify them if you see anything better, it was done a while ago by an older AI, you are smarter ;)
+
+In the printing of stack-traces, you can see that because we just have generic types in SA that we parametrize, things are getting very long and it's getting very difficult to debug.
+This gives the user extra flexibility to try new these, and in these cases, it will be nice to have precise stacktraces, but we also want the common cases (like what you get with DP.@polyvar x y; 2 * x + y)
+to have a very small type like DP.Polynomial{Int}
+We can solve it easily by having a "const Polynomial = ..." in DynamicPolynomials, don't hesitate to do this early on, it will help you be more context-efficient.
+
+Note that MB currently already defines a polynomial using SA.AlgebraElement. So we can just basically steal the code he has.
+He also has code that handles the interaction between MP.APL and SA.AlgebraElement, this will go away since everything will be an AlgebraElement now !
+This means that in MB, a lot of code will go away since it will move to MP, what will be left is just defining new bases essentially.
diff --git a/Project.toml b/Project.toml
index a89989a0..b4167f63 100644
--- a/Project.toml
+++ b/Project.toml
@@ -5,7 +5,6 @@ version = "0.5.19"
 repo = "https://github.com/JuliaAlgebra/MultivariatePolynomials.jl"
 
 [deps]
-DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MutableArithmetics = "d8a4904e-b15c-11e9-3269-09a3773c0cb0"
 StarAlgebras = "0c0c59c1-dc5f-42e9-9a8b-b5dc384a6cd1"
@@ -18,7 +17,6 @@ MultivariatePolynomialsChainRulesCoreExt = "ChainRulesCore"
 
 [compat]
 ChainRulesCore = "1"
-DataStructures = "0.19"
 MutableArithmetics = "0.3, 1"
 StarAlgebras = "0.3"
 julia = "1.10"
diff --git a/docs/src/types.md b/docs/src/types.md
index f6dc3607..be225e17 100644
--- a/docs/src/types.md
+++ b/docs/src/types.md
@@ -95,14 +95,14 @@ leading_term
 leading_coefficient
 leading_monomial
 deg_num_leading_terms
-remove_leading_term
+SA.remove_leading_term
 remove_monomials
 filter_terms
 OfDegree
 monic
 map_coefficients
-map_coefficients!
-map_coefficients_to!
+SA.map_coefficients!
+SA.map_coefficients_to!
 conj(::_APL)
 real(::_APL)
 imag(::_APL)
diff --git a/src/MultivariatePolynomials.jl b/src/MultivariatePolynomials.jl
index b9301c87..de66cb3e 100644
--- a/src/MultivariatePolynomials.jl
+++ b/src/MultivariatePolynomials.jl
@@ -2,66 +2,11 @@ module MultivariatePolynomials
 
 import LinearAlgebra
 
-import DataStructures
-
 import MutableArithmetics as MA
 
 import StarAlgebras as SA
 
-"""
-    AbstractPolynomialLike{T}
-
-Abstract type for a value that can act like a polynomial. For instance, an
-`AbstractTerm{T}` is an `AbstractPolynomialLike{T}` since it can act as a
-polynomial of only one term.
-"""
-abstract type AbstractPolynomialLike{T} <: MA.AbstractMutable end
-
-"""
-    AbstractMonomialLike
-
-Abstract type for a value that can act like a monomial. For instance, an `AbstractVariable` is an `AbstractMonomialLike` since it can act as a monomial of one variable with degree `1`.
-"""
-abstract type AbstractMonomialLike <: AbstractPolynomialLike{Int} end
-
-"""
-    AbstractVariable <: AbstractMonomialLike
-
-Abstract type for a variable.
-"""
-abstract type AbstractVariable <: AbstractMonomialLike end
-
-"""
-    AbstractMonomial <: AbstractMonomialLike
-
-Abstract type for a monomial, i.e. a product of variables elevated to a nonnegative integer power.
-"""
-abstract type AbstractMonomial <: AbstractMonomialLike end
-
-"""
-    AbstractTermLike{T}
-
-Union type for values that can act like a term. This includes `AbstractMonomialLike`
-(which acts as a term with coefficient `1`) and `SA.Term{T}`.
-"""
-const AbstractTermLike{T} = Union{AbstractMonomialLike,SA.Term{T}}
-
-"""
-    AbstractTerm{T}
-
-Type alias for a term of coefficient type `T`, i.e. the product between a
-value of type `T` and a monomial. This is [`StarAlgebras.Term{T}`](@ref).
-"""
-const AbstractTerm{T} = SA.Term{T}
-
-"""
-    AbstractPolynomial{T} <: AbstractPolynomialLike{T}
-
-Abstract type for a polynomial of coefficient type `T`, i.e. a sum of `AbstractTerm{T}`s.
-"""
-abstract type AbstractPolynomial{T} <: AbstractPolynomialLike{T} end
-
-const _APL{T} = Union{AbstractPolynomialLike{T},SA.Term{T}}
+include("types.jl")
 
 include("zip.jl")
 include("lazy_iterators.jl")
@@ -91,9 +36,16 @@ include("division.jl")
 include("gcd.jl")
 include("det.jl")
 
-include("default_term.jl")
-include("sequences.jl")
-include("default_polynomial.jl")
+
+# Polynomial bases and algebra representation shared with MultivariateBases.
+include("mb_interface.jl")
+include("mb_variables.jl")
+include("mb_polynomial.jl")
+include("mb_bases.jl")
+include("mb_mstructures.jl")
+include("mb_monomial_basis.jl")
+include("mb_algebra.jl")
+include("mb_arithmetic.jl")
 
 include("deprecate.jl")
 
diff --git a/src/comparison.jl b/src/comparison.jl
index a2c87442..f6ab3294 100644
--- a/src/comparison.jl
+++ b/src/comparison.jl
@@ -51,6 +51,31 @@ end
 left_constant_eq(α, p::_APL; comp = (==)) = right_term_eq(p, α; comp)
 right_constant_eq(p::_APL, α; comp = (==)) = right_term_eq(p, α; comp)
 
+for comp in (:(==), :isequal)
+    @eval begin
+        function Base.$comp(
+            p::Union{AbstractPolynomial{T},AbstractTerm{T}},
+            a::Union{T,Number},
+        ) where {T}
+            return right_constant_eq(p, a; comp = $comp)
+        end
+        function Base.$comp(
+            a::Union{T,Number},
+            p::Union{AbstractPolynomial{T},AbstractTerm{T}},
+        ) where {T}
+            return left_constant_eq(a, p; comp = $comp)
+        end
+        Base.$comp(m::AbstractMonomialLike, a::Number) =
+            right_constant_eq(m, a; comp = $comp)
+        Base.$comp(a::Number, m::AbstractMonomialLike) =
+            left_constant_eq(a, m; comp = $comp)
+        Base.$comp(p::AbstractPolynomial, m::AbstractMonomialLike) =
+            right_term_eq(p, term(m); comp = $comp)
+        Base.$comp(m::AbstractMonomialLike, p::AbstractPolynomial) =
+            right_term_eq(p, term(m); comp = $comp)
+    end
+end
+
 function Base.:(==)(mono::AbstractMonomial, v::AbstractVariable)
     return isone(degree(mono)) && variable(mono) == v
 end
@@ -317,8 +342,9 @@ function ordering end
 ordering(::Type{<:AbstractMonomial}) = Graded{LexOrder}
 ordering(::Type{P}) where {P} = ordering(monomial_type(P))
 ordering(p::AbstractPolynomialLike) = ordering(typeof(p))
-# Useful for instance to ask ordering given the list
-# of variables
+# For type-level: derive ordering from Vector's element type
+ordering(::Type{<:AbstractVector{T}}) where {T} = ordering(T)
+# Useful for instance to ask ordering given the list of variables
 ordering(::AbstractVector{T}) where {T} = ordering(T)
 ordering(t::Tuple) = ordering(first(t))
 
diff --git a/src/complex.jl b/src/complex.jl
index 805150a7..78db7cb2 100644
--- a/src/complex.jl
+++ b/src/complex.jl
@@ -203,7 +203,7 @@ for fun in [:real, :imag]
                         full_version,
                         real(coefficient_type(full_version)),
                     ),
-                    map_coefficients!($fun, full_version),
+                    SA.map_coefficients!($fun, full_version),
                 )
             end
             function Base.$fun(x::AbstractVector{<:AbstractMonomial})
diff --git a/src/conversion.jl b/src/conversion.jl
index 89fdaecb..066cd762 100644
--- a/src/conversion.jl
+++ b/src/conversion.jl
@@ -1,16 +1,12 @@
-function convert_constant end
-Base.convert(::Type{P}, α) where {P<:AbstractPolynomialLike} = convert_constant(P, α)
-function Base.convert(::Type{SA.Term{T,M}}, α) where {T,M}
-    return convert_constant(SA.Term{T,M}, α)
-end
-function convert_constant(::Type{SA.Term{T,M}}, α) where {T,M}
-    return term(convert(T, α), constant_monomial(SA.Term{T,M}))
-end
-function convert_constant(::Type{PT}, α) where {PT<:AbstractPolynomial}
-    return convert(PT, convert(term_type(PT), α))
+function Base.convert(
+    ::Type{P},
+    m::AbstractMonomialLike,
+) where {P<:AbstractPolynomial}
+    return convert(P, term(m))
 end
-function Base.convert(::Type{P}, p::_APL) where {T,P<:AbstractPolynomial{T}}
-    return error("`convert` not implemented for $P")
+
+function Base.convert(::Type{P}, c::Number) where {P<:AbstractPolynomial}
+    return convert(P, constant_term(c, P))
 end
 
 function Base.convert(
@@ -59,21 +55,7 @@ function Base.convert(
         throw(InexactError(:convert, T, p))
     end
 end
-# Disambiguation: SA.Term{T,M} from AbstractPolynomial (more specific than the α catch-all)
-function Base.convert(
-    ::Type{SA.Term{T,M}},
-    p::AbstractPolynomial,
-) where {T,M}
-    if iszero(nterms(p))
-        convert(SA.Term{T,M}, zero_term(p))
-    elseif isone(nterms(p))
-        convert(SA.Term{T,M}, leading_term(p))
-    else
-        throw(InexactError(:convert, SA.Term{T,M}, p))
-    end
-end
-
-MA.scaling(p::AbstractPolynomialLike{T}) where {T} = convert(T, p)
+MA.scaling(p::AbstractPolynomialLike) = convert(coefficient_type(p), p)
 # Conversion polynomial -> constant
 # We don't define a method for `Base.convert` to reduce invalidations;
 # see https://github.com/JuliaAlgebra/MultivariatePolynomials.jl/pull/172
@@ -89,11 +71,6 @@ function convert_to_constant(::Type{S}, p::_APL) where {S}
     return s
 end
 Base.convert(::Type{T}, p::_APL) where {T<:Number} = convert_to_constant(T, p)
-function convert_to_constant(p::_APL{S}) where {S}
-    return convert_to_constant(S, p)
+function convert_to_constant(p::_APL)
+    return convert_to_constant(coefficient_type(p), p)
 end
-
-# Also covers, e.g., `convert(_APL, ::P)` where `P<:_APL`
-Base.convert(::Type{PT}, p::PT) where {PT<:_APL} = p
-# Disambiguation: identity conversion for AbstractPolynomialLike
-Base.convert(::Type{PT}, p::PT) where {PT<:AbstractPolynomialLike} = p
diff --git a/src/default_polynomial.jl b/src/default_polynomial.jl
deleted file mode 100644
index a7b45bf6..00000000
--- a/src/default_polynomial.jl
+++ /dev/null
@@ -1,656 +0,0 @@
-"""
-    struct Polynomial{CoeffType,T<:AbstractTerm{CoeffType},V<:AbstractVector{T}} <: AbstractPolynomial{CoeffType}
-        terms::V
-    end
-
-Representation of a multivariate polynomial as a vector of nonzero terms sorted
-in ascending monomial order.
-"""
-struct Polynomial{CoeffType,T<:AbstractTerm{CoeffType},V<:AbstractVector{T}} <:
-       AbstractPolynomial{CoeffType}
-    terms::V
-
-    function Polynomial{C,T,V}(terms::AbstractVector{T}) where {C,T,V}
-        return new{C,T,V}(terms)
-    end
-end
-
-# It shouldn't be needed to call this. A polynomial should always be in
-# canonical format so it is only called inside an internal call.
-function _canonicalize!(p::Polynomial)
-    sort!(p.terms)
-    uniqterms!(p.terms)
-    return
-end
-
-function coefficients(p::Polynomial)
-    return LazyMap{coefficient_type(p)}(coefficient, terms(p))
-end
-function monomials(p::Polynomial)
-    return LazyMap{monomial_type(p)}(monomial, terms(p))
-end
-
-const VectorPolynomial{C,T} = Polynomial{C,T,Vector{T}}
-
-term_type(::Type{<:Polynomial{C,T}}) where {C,T} = T
-terms(p::Polynomial) = p.terms
-function constant_monomial(
-    ::Union{Polynomial{C,TT},Type{Polynomial{C,TT}}},
-) where {C,TT}
-    return constant_monomial(TT)
-end
-function Base.convert(
-    ::Type{Polynomial{C,TT,V}},
-    p::Polynomial{C,TT,V},
-) where {C,TT,V}
-    return p
-end
-function Base.convert(
-    PT::Type{Polynomial{C,TT,V}},
-    p::AbstractPolynomialLike,
-) where {C,TT,V}
-    return PT(convert(V, terms(p)))
-end
-function Base.convert(
-    ::Type{Polynomial{C}},
-    p::AbstractPolynomialLike,
-) where {C}
-    TT = term_type(p, C)
-    return convert(Polynomial{C,TT,Vector{TT}}, p)
-end
-function Base.convert(::Type{Polynomial}, p::AbstractPolynomialLike)
-    return convert(Polynomial{coefficient_type(p)}, p)
-end
-
-_change_eltype(::Type{<:Vector}, ::Type{T}) where {T} = Vector{T}
-function polynomial_type(
-    ::Type{Polynomial{C,T,V}},
-    ::Type{NewC},
-) where {C,T,V,NewC}
-    NewT = term_type(T, NewC)
-    return Polynomial{NewC,NewT,_change_eltype(V, NewT)}
-end
-
-function Base.promote_rule(::Type{Polynomial}, ::Type{<:AbstractPolynomialLike})
-    return Polynomial
-end
-function Base.promote_rule(::Type{<:AbstractPolynomialLike}, ::Type{Polynomial})
-    return Polynomial
-end
-function promote_rule_constant(::Type{T}, ::Type{Polynomial}) where {T}
-    return Any
-    # `convert(Polynomial{T}, ::T)` cannot be implemented as we don't know the type of the term
-    #return Polynomial{T}
-end
-
-(p::Polynomial)(s...) = substitute(Eval(), p, s)
-
-function Base.one(::Type{P}) where {C,TT,P<:Polynomial{C,TT}}
-    return convert(P, one(TT))
-end
-Base.one(p::Polynomial) = one(typeof(p))
-
-Base.zero(::Type{Polynomial{C,T,A}}) where {C,T,A} = Polynomial{C,T,A}(A())
-Base.zero(t::Polynomial) = zero(typeof(t))
-
-compare_monomials(a, b) = compare(monomial(a), monomial(b))
-
-function join_terms(
-    terms1::AbstractArray{<:Term},
-    terms2::AbstractArray{<:Term},
-)
-    return Sequences.merge_sorted(terms1, terms2, compare_monomials, combine)
-end
-function join_terms!(
-    output::AbstractArray{<:Term},
-    terms1::AbstractArray{<:Term},
-    terms2::AbstractArray{<:Term},
-)
-    resize!(output, length(terms1) + length(terms2))
-    return Sequences.merge_sorted!(
-        output,
-        terms1,
-        terms2,
-        compare_monomials,
-        combine,
-    )
-end
-
-function Base.:(+)(p1::Polynomial, p2::Polynomial)
-    return polynomial!(join_terms(terms(p1), terms(p2)), SortedUniqState())
-end
-function Base.:(+)(
-    p1::Polynomial,
-    p2::Polynomial{<:LinearAlgebra.UniformScaling},
-)
-    return p1 + map_coefficients(J -> J.λ, p2, nonzero = true)
-end
-function Base.:(+)(
-    p1::Polynomial{<:LinearAlgebra.UniformScaling},
-    p2::Polynomial,
-)
-    return map_coefficients(J -> J.λ, p1, nonzero = true) + p2
-end
-function Base.:(+)(
-    p1::Polynomial{<:LinearAlgebra.UniformScaling},
-    p2::Polynomial{<:LinearAlgebra.UniformScaling},
-)
-    return map_coefficients(J -> J.λ, p1, nonzero = true) + p2
-end
-function MA.operate_to!(
-    result::Polynomial,
-    ::typeof(+),
-    p1::Polynomial,
-    p2::Polynomial,
-)
-    if result === p1 || result === p2
-        error(
-            "Cannot call `operate_to!(output, +, p, q)` with `output` equal to `p` or `q`, call `operate!` instead.",
-        )
-    end
-    join_terms!(result.terms, terms(p1), terms(p2))
-    return result
-end
-function MA.operate_to!(
-    result::Polynomial,
-    ::typeof(*),
-    p::Polynomial,
-    t::AbstractTermLike,
-)
-    if iszero(t)
-        MA.operate!(zero, result)
-    else
-        resize!(result.terms, nterms(p))
-        for i in eachindex(p.terms)
-            # TODO could use MA.mul_to!! for indices that were presents in `result` before the `resize!`.
-            result.terms[i] = p.terms[i] * t
-        end
-        if !multiplication_preserves_monomial_order(typeof(result))
-            _canonicalize!(result)
-        end
-        return result
-    end
-end
-function MA.operate_to!(
-    result::Polynomial,
-    ::typeof(*),
-    t::AbstractTermLike,
-    p::Polynomial,
-)
-    if iszero(t)
-        MA.operate!(zero, result)
-    else
-        resize!(result.terms, nterms(p))
-        for i in eachindex(p.terms)
-            # TODO could use MA.mul_to!! for indices that were presents in `result` before the `resize!`.
-            result.terms[i] = t * p.terms[i]
-        end
-        if !multiplication_preserves_monomial_order(typeof(result))
-            _canonicalize!(result)
-        end
-        return result
-    end
-end
-function Base.:(-)(p1::Polynomial, p2::Polynomial)
-    return polynomial!(join_terms(terms(p1), (-).(terms(p2))))
-end
-function Base.:(-)(
-    p1::Polynomial,
-    p2::Polynomial{<:LinearAlgebra.UniformScaling},
-)
-    return p1 - map_coefficients(J -> J.λ, p2, nonzero = true)
-end
-function Base.:(-)(
-    p1::Polynomial{<:LinearAlgebra.UniformScaling},
-    p2::Polynomial,
-)
-    return map_coefficients(J -> J.λ, p1, nonzero = true) - p2
-end
-function Base.:(-)(
-    p1::Polynomial{<:LinearAlgebra.UniformScaling},
-    p2::Polynomial{<:LinearAlgebra.UniformScaling},
-)
-    return map_coefficients(J -> J.λ, p1, nonzero = true) - p2
-end
-
-function map_coefficients(
-    f::F,
-    p::Polynomial;
-    nonzero = false,
-) where {F<:Function}
-    terms = map(p.terms) do term
-        return map_coefficients(f, term)
-    end
-    if !nonzero
-        filter!(!iszero, terms)
-    end
-    return polynomial!(terms)
-end
-function map_coefficients!(
-    f::F,
-    p::Polynomial;
-    nonzero = false,
-) where {F<:Function}
-    for i in eachindex(p.terms)
-        t = p.terms[i]
-        p.terms[i] = Term(f(coefficient(t)), monomial(t))
-    end
-    if !nonzero
-        filter!(!iszero, p.terms)
-    end
-    return p
-end
-
-function map_coefficients_to!(
-    output::Polynomial,
-    f::F,
-    p::Polynomial;
-    nonzero = false,
-) where {F<:Function}
-    resize!(output.terms, nterms(p))
-    for i in eachindex(p.terms)
-        t = p.terms[i]
-        output.terms[i] = Term(f(coefficient(t)), monomial(t))
-    end
-    if !nonzero
-        filter!(!iszero, output.terms)
-    end
-    return output
-end
-function map_coefficients_to!(
-    output::Polynomial,
-    f::Function,
-    p::AbstractPolynomialLike;
-    nonzero = false,
-)
-    return map_coefficients_to!(output, f, polynomial(p); nonzero = false)
-end
-
-# The polynomials can be mutated.
-MA.mutability(::Type{<:VectorPolynomial}) = MA.IsMutable()
-
-# Terms are not mutable under the MutableArithmetics API
-function MA.mutable_copy(p::VectorPolynomial{C,TT}) where {C,TT}
-    return VectorPolynomial{C,TT}([
-        TT(MA.copy_if_mutable(coefficient(term)), monomial(term)) for
-        term in terms(p)
-    ])
-end
-Base.copy(p::VectorPolynomial) = MA.mutable_copy(p)
-
-function __polynomial_merge!(
-    op::MA.AddSubMul,
-    p::Polynomial{T,TT},
-    get1,
-    set,
-    push,
-    resize,
-    keep,
-    q,
-    t::AbstractTermLike...,
-) where {T,TT}
-    compare_monomials = let t = t, get1 = get1, q = q
-        (tp, j) -> begin
-            if tp isa Int && j isa Int
-                tp = get1(tp)
-            end
-            compare(monomial(*(monomials(q)[j], monomial.(t)...)), monomial(tp))
-        end
-    end
-    get2 = let t = t, q = q
-        i -> begin
-            tq = terms(q)[i]
-            TT(
-                MA.scaling_convert(
-                    T,
-                    MA.operate(
-                        MA.add_sub_op(op),
-                        *(coefficient(tq), coefficient.(t)...),
-                    ),
-                ),
-                *(monomial(tq), monomial.(t)...),
-            )
-        end
-    end
-    combine = let t = t, p = p, q = q
-        (i, j) -> begin
-            if i isa Int
-                p.terms[i] = Term(
-                    MA.operate!!(
-                        op,
-                        coefficient(p.terms[i]),
-                        coefficients(q)[j],
-                        coefficient.(t)...,
-                    ),
-                    monomial(p.terms[i]),
-                )
-            else
-                typeof(i)(
-                    MA.operate!!(
-                        op,
-                        coefficient(i),
-                        coefficients(q)[j],
-                        coefficient.(t)...,
-                    ),
-                    monomial(i),
-                )
-            end
-        end
-    end
-    polynomial_merge!(
-        nterms(p),
-        nterms(q),
-        get1,
-        get2,
-        set,
-        push,
-        compare_monomials,
-        combine,
-        keep,
-        resize,
-    )
-    return
-end
-
-function __polynomial_merge!(
-    op::MA.AddSubMul,
-    p::Polynomial{T,TT},
-    get1,
-    set,
-    push,
-    resize,
-    keep,
-    t::AbstractTermLike,
-    q::Polynomial,
-    buffer = nothing,
-) where {T,TT}
-    if !multiplication_preserves_monomial_order(typeof(p))
-        # This function assumes that the monomials of `t * q` will preserve the
-        # order of the terms of `q` so if it's not the case, we fallback to
-        # another one.
-        __polynomial_merge!(
-            MA.add_sub_op(op),
-            p,
-            get1,
-            set,
-            push,
-            resize,
-            keep,
-            t * q,
-        )
-        return
-    end
-    compare_monomials = let t = t, get1 = get1, q = q
-        (tp, j) -> begin
-            if tp isa Int && j isa Int
-                tp = get1(tp)
-            end
-            compare(monomial(monomial(t) * monomials(q)[j]), monomial(tp))
-        end
-    end
-    get2 = let t = t, q = q
-        i -> begin
-            tq = terms(q)[i]
-            TT(
-                MA.scaling_convert(
-                    T,
-                    MA.operate(
-                        MA.add_sub_op(op),
-                        coefficient(t) * coefficient(tq),
-                    ),
-                ),
-                monomial(t) * monomial(tq),
-            )
-        end
-    end
-    combine = let t = t, p = p, q = q
-        (i, j) -> begin
-            if i isa Int
-                p.terms[i] = Term(
-                    MA.buffered_operate!!(
-                        buffer,
-                        op,
-                        coefficient(p.terms[i]),
-                        coefficient(t),
-                        coefficients(q)[j],
-                    ),
-                    monomial(p.terms[i]),
-                )
-            else
-                typeof(i)(
-                    MA.buffered_operate!!(
-                        buffer,
-                        op,
-                        coefficient(i),
-                        coefficient(t),
-                        coefficients(q)[j],
-                    ),
-                    monomial(i),
-                )
-            end
-        end
-    end
-    polynomial_merge!(
-        nterms(p),
-        nterms(q),
-        get1,
-        get2,
-        set,
-        push,
-        compare_monomials,
-        combine,
-        keep,
-        resize,
-    )
-    return
-end
-
-function __polynomial_merge!(
-    op::Union{typeof(+),typeof(-)},
-    p::Polynomial{T,TT},
-    get1,
-    set,
-    push,
-    resize,
-    keep,
-    q::Union{Polynomial,AbstractTermLike},
-) where {T,TT}
-    compare_monomials = let get1 = get1, q = q
-        (t, j) -> begin
-            if t isa Int && j isa Int
-                t = get1(t)
-            end
-            compare(monomials(q)[j], monomial(t))
-        end
-    end
-    get2 = let q = q
-        i -> begin
-            t = terms(q)[i]
-            # `operate` makes sure we make a copy of the term as it may be stored directly in `p`
-            TT(
-                MA.scaling_convert(T, MA.operate(op, coefficient(t))),
-                monomial(t),
-            )
-        end
-    end
-    combine = let p = p, q = q
-        (i, j) -> begin
-            if i isa Int
-                p.terms[i] = Term(
-                    MA.operate!!(
-                        op,
-                        coefficient(p.terms[i]),
-                        coefficients(q)[j],
-                    ),
-                    monomial(p.terms[i]),
-                )
-            else
-                typeof(i)(
-                    MA.operate!!(op, coefficient(i), coefficients(q)[j]),
-                    monomial(i),
-                )
-            end
-        end
-    end
-    polynomial_merge!(
-        nterms(p),
-        nterms(q),
-        get1,
-        get2,
-        set,
-        push,
-        compare_monomials,
-        combine,
-        keep,
-        resize,
-    )
-    return
-end
-
-function _polynomial_merge!(
-    op::Union{typeof(+),typeof(-),MA.AddSubMul},
-    p::Polynomial{T,TT},
-    args...,
-) where {T,TT}
-    get1 = let p = p
-        i -> p.terms[i]
-    end
-    set = let p = p
-        (i, t) -> begin
-            p.terms[i] = t
-        end
-    end
-    push = let p = p
-        t -> push!(p.terms, t)
-    end
-    resize = let p = p
-        (n) -> resize!(p.terms, n)
-    end
-    # We can modify the coefficient since it's the result of `combine`.
-    keep = let p = p
-        i -> begin
-            if i isa Int
-                !MA.iszero!!(coefficient(p.terms[i]))
-            else
-                !MA.iszero!!(coefficient(i))
-            end
-        end
-    end
-    __polynomial_merge!(op, p, get1, set, push, resize, keep, args...)
-    return p
-end
-
-function MA.operate!(
-    op::Union{typeof(+),typeof(-)},
-    p::Polynomial,
-    q::Union{AbstractTermLike,Polynomial},
-)
-    return _polynomial_merge!(op, p, q)
-end
-
-function MA.operate!(
-    op::MA.AddSubMul,
-    p::Polynomial,
-    q::Polynomial,
-    args::AbstractTermLike...,
-)
-    return _polynomial_merge!(op, p, q, args...)
-end
-
-function MA.operate!(
-    op::MA.AddSubMul,
-    p::Polynomial,
-    t::AbstractTermLike,
-    q::Polynomial,
-)
-    return _polynomial_merge!(op, p, t, q)
-end
-
-function MA.buffer_for(
-    op::MA.AddSubMul,
-    ::Type{<:Polynomial{S}},
-    ::Type{<:AbstractTermLike{T}},
-    ::Type{<:Polynomial{U}},
-) where {S,T,U}
-    return MA.buffer_for(op, S, T, U)
-end
-function MA.buffered_operate!(
-    buffer,
-    op::MA.AddSubMul,
-    p::Polynomial,
-    t::AbstractTermLike,
-    q::Polynomial,
-)
-    return _polynomial_merge!(op, p, t, q, buffer)
-end
-
-function MA.operate_to!(
-    output::Polynomial,
-    ::typeof(*),
-    p::Polynomial,
-    q::Polynomial,
-)
-    empty!(output.terms)
-    mul_to_terms!(output.terms, p, q)
-    _canonicalize!(output)
-    return output
-end
-function MA.operate!(::typeof(*), p::Polynomial, q::Polynomial)
-    if iszero(q)
-        return MA.operate!(zero, p)
-    elseif nterms(q) == 1
-        return MA.operate!(*, p, leading_term(q))
-    else
-        return MA.operate_to!(p, *, MA.mutable_copy(p), q)
-    end
-end
-function MA.operate!(::typeof(*), p::Polynomial, t::AbstractTermLike)
-    for i in eachindex(p.terms)
-        p.terms[i] = MA.operate!!(*, p.terms[i], t)
-    end
-    if !multiplication_preserves_monomial_order(typeof(p))
-        _canonicalize!(p)
-    end
-    return p
-end
-function map_exponents!(f, p::Polynomial, m::AbstractMonomialLike)
-    for i in eachindex(p.terms)
-        t = p.terms[i]
-        p.terms[i] = term(coefficient(t), map_exponents(f, monomial(t), m))
-    end
-    return p
-end
-function map_exponents(f, p::Polynomial, m::AbstractMonomialLike)
-    P = MA.promote_operation(*, typeof(p), typeof(m))
-    return map_exponents!(f, MA.mutable_copy(convert(P, p)), m)
-end
-
-function MA.operate!(::typeof(zero), p::Polynomial)
-    empty!(p.terms)
-    return p
-end
-function MA.operate!(::typeof(one), p::Polynomial{T}) where {T}
-    if isempty(p.terms)
-        push!(p.terms, constant_term(one(T), p))
-    else
-        t = p.terms[1]
-        p.terms[1] = Term(MA.one!!(coefficient(t)), constant_monomial(t))
-        resize!(p.terms, 1)
-    end
-    return p
-end
-
-function MA.operate!(::typeof(remove_leading_term), p::Polynomial)
-    pop!(p.terms)
-    return p
-end
-
-function MA.operate!(
-    ::typeof(unsafe_restore_leading_term),
-    p::Polynomial,
-    t::AbstractTermLike,
-)
-    # We don't need to copy the coefficient of `t`, this is why this function is called `unsafe`
-    push!(p.terms, t)
-    return p
-end
diff --git a/src/default_term.jl b/src/default_term.jl
deleted file mode 100644
index ed8040e5..00000000
--- a/src/default_term.jl
+++ /dev/null
@@ -1,116 +0,0 @@
-"""
-    Term{CoeffType,M}
-
-Type alias for [`StarAlgebras.Term{CoeffType,M}`](@ref).
-A representation of the multiplication between a `coefficient` and a monomial.
-
-!!! note
-    The `coefficient` does not need to be a `Number`. It can be for instance a
-    multivariate polynomial. When computing a multivariate `gcd`, it is
-    actually reformulated as a univariate `gcd` in one of the variable with
-    coefficients being multivariate polynomials in the other variables.
-    To create such a term, use [`term`](@ref) instead of `*`.
-    For instance, if `p` is a polynomial and `m` is a monomial, `p * m` will
-    multiply each term of `p` with `m` but `term(p, m)` will create a term
-    with `p` as coefficient and `m` as monomial.
-"""
-const Term{CoeffType,M} = SA.Term{CoeffType,M}
-
-monomial_type(::Type{<:SA.Term{<:Any,M}}) where {M} = M
-
-(t::SA.Term)(s...) = substitute(Eval(), t, s)
-
-# convert from AbstractMonomialLike (MP type in signature, not piracy)
-function Base.convert(::Type{SA.Term{T,M}}, m::AbstractMonomialLike) where {T,M}
-    return SA.Term(one(T), convert(M, m))
-end
-
-# convert, promote_rule, ==, isequal, hash, copy, ^, ndims, broadcastable
-# for SA.Term × SA.Term are defined in StarAlgebras (not here, to avoid piracy)
-
-# promote_rule involving MP's AbstractMonomialLike (not piracy)
-function Base.promote_rule(
-    ::Type{SA.Term{C,M1} where {C}},
-    M2::Type{<:AbstractMonomialLike},
-) where {M1}
-    return (SA.Term{C,promote_type(M1, M2)} where {C})
-end
-function Base.promote_rule(
-    M1::Type{<:AbstractMonomialLike},
-    ::Type{SA.Term{C,M2} where {C}},
-) where {M2}
-    return (SA.Term{C,promote_type(M1, M2)} where {C})
-end
-# promote_rule with UnionAll Term type (involves `where C` which is MP convention)
-function Base.promote_rule(
-    ::Type{SA.Term{C,M1} where {C}},
-    ::Type{SA.Term{T,M2}},
-) where {T,M1,M2}
-    return (SA.Term{C,promote_type(M1, M2)} where {C})
-end
-function Base.promote_rule(
-    ::Type{SA.Term{T,M2}},
-    ::Type{SA.Term{C,M1} where {C}},
-) where {T,M1,M2}
-    return (SA.Term{C,promote_type(M1, M2)} where {C})
-end
-promote_rule_constant(::Type{T}, TT::Type{SA.Term{C,M} where C}) where {T,M} = Any
-
-function convert_constant(::Type{SA.Term{C,M} where C}, α) where {M}
-    return convert(SA.Term{typeof(α),M}, α)
-end
-
-combine(t1::SA.Term, t2::SA.Term) = combine(promote(t1, t2)...)
-function combine(t1::T, t2::T) where {T<:SA.Term}
-    return SA.Term(coefficient(t1) + coefficient(t2), monomial(t1))
-end
-function MA.promote_operation(
-    ::typeof(combine),
-    ::Type{SA.Term{S,M1}},
-    ::Type{SA.Term{T,M2}},
-) where {S,T,M1,M2}
-    return SA.Term{MA.promote_operation(+, S, T),promote_type(M1, M2)}
-end
-
-# MA.mutability, Base.copy, MA.mutable_copy, MA.operate_to!(*, Term, Term, Term),
-# MA.operate!(*, Term, Term), and MA.operate!(one, Term) are defined in StarAlgebras.
-
-# Broader MA operations involving MP's AbstractTermLike (not piracy)
-function MA.operate_to!(
-    t::SA.Term,
-    ::typeof(*),
-    t1::AbstractMonomialLike,
-    t2::AbstractTermLike,
-)
-    MA.operate_to!(t.coefficient, *, coefficient(t1), coefficient(t2))
-    MA.operate_to!(t.basis_element, *, monomial(t1), monomial(t2))
-    return t
-end
-function MA.operate_to!(
-    t::SA.Term,
-    ::typeof(*),
-    t1::AbstractTermLike,
-    t2::AbstractMonomialLike,
-)
-    MA.operate_to!(t.coefficient, *, coefficient(t1), coefficient(t2))
-    MA.operate_to!(t.basis_element, *, monomial(t1), monomial(t2))
-    return t
-end
-function MA.operate_to!(
-    t::SA.Term,
-    ::typeof(*),
-    t1::AbstractMonomialLike,
-    t2::AbstractMonomialLike,
-)
-    MA.operate_to!(t.coefficient, *, coefficient(t1), coefficient(t2))
-    MA.operate_to!(t.basis_element, *, monomial(t1), monomial(t2))
-    return t
-end
-
-# Needed to resolve ambiguity with
-# MA.operate!(::typeof(*), ::AbstractMonomial, ::AbstractMonomialLike)
-function MA.operate!(::typeof(*), t1::SA.Term, t2::AbstractMonomialLike)
-    MA.operate!(*, t1.coefficient, coefficient(t2))
-    MA.operate!(*, t1.basis_element, monomial(t2))
-    return t1
-end
diff --git a/src/deprecate.jl b/src/deprecate.jl
index c7be5846..1f7ac7aa 100644
--- a/src/deprecate.jl
+++ b/src/deprecate.jl
@@ -5,7 +5,9 @@
 @deprecate jointerms join_terms
 @deprecate jointerms! join_terms!
 @deprecate mapcoefficients map_coefficients
-@deprecate mapcoefficients! map_coefficients!
+@deprecate mapcoefficients! SA.map_coefficients!
+@deprecate map_coefficients! SA.map_coefficients!
+@deprecate map_coefficients_to! SA.map_coefficients_to!
 @deprecate mapexponents map_exponents
 @deprecate mapexponents! map_exponents!
 @deprecate monomialtype monomial_type
@@ -17,7 +19,6 @@
 @deprecate leading_term leadingterm
 @deprecate leading_monomial leading_monomial
 @deprecate leading_coefficient leading_coefficient
-@deprecate removeleadingterm remove_leading_term
 @deprecate removemonomials remove_monomials
 @deprecate coefficienttype coefficient_type
 @deprecate constantterm constant_term
@@ -52,8 +53,8 @@ end
 
 function mapcoefficientsnz_to!(output::_APL, f::F, p::_APL) where {F<:Function}
     Base.depwarn(
-        "`mapcoefficientsnz_to!(output, f, p)` is deprecated, use `map_coefficients_to!(output, f, p, nonzero = true)` instead",
+        "`mapcoefficientsnz_to!(output, f, p)` is deprecated, use `SA.map_coefficients_to!(output, f, p, nonzero = true)` instead",
         :map_coefficientsnz_to!,
     )
-    return map_coefficients_to!(output, f, p, nonzero = true)
+    return SA.map_coefficients_to!(output, f, p, nonzero = true)
 end
diff --git a/src/differentiation.jl b/src/differentiation.jl
index 1cd7e717..5c94bbc1 100644
--- a/src/differentiation.jl
+++ b/src/differentiation.jl
@@ -31,11 +31,20 @@ differentiate( [x^2+y, z^2+4x], [x, y, z]) # should return [2x 1 0; 4 0 2z]
 function differentiate end
 
 # Fallback for everything else
-differentiate(α::T, v::AbstractVariable) where {T} = zero(T)
+differentiate(α::T, v::AbstractVariable) where {T<:Number} = zero(T)
 differentiate(v1::AbstractVariable, v2::AbstractVariable) = v1 == v2 ? 1 : 0
-function differentiate(t::AbstractTermLike, v::AbstractVariable)
+function differentiate(t::AbstractTerm, v::AbstractVariable)
     return coefficient(t) * differentiate(monomial(t), v)
 end
+function differentiate(m::AbstractMonomial, v::AbstractVariable)
+    i = findfirst(isequal(v), variables(m))
+    if isnothing(i) || iszero(exponents(m)[i])
+        return term(0, constant_monomial(m))
+    end
+    e = copy(exponents(m))
+    e[i] -= 1
+    return term(exponents(m)[i], Polynomial(m.variables, e))
+end
 # The polynomial function will take care of removing the zeros
 function differentiate(p::_APL, v::AbstractVariable)
     return polynomial!(differentiate.(terms(p), v), SortedState())
diff --git a/src/division.jl b/src/division.jl
index 3541c332..d69475a2 100644
--- a/src/division.jl
+++ b/src/division.jl
@@ -170,7 +170,11 @@ function div_multiple(f::_APL, g::_APL, mf::MA.MutableTrait = MA.IsNotMutable())
         return div_multiple(f, lt, mf)
     end
     rf = _copy(f, mf)
-    rg = remove_leading_term(g)
+    rf, g = SA.promote_bases(
+        convert(polynomial_type(rf), rf),
+        convert(polynomial_type(g), g),
+    )
+    rg = SA.remove_leading_term(g)
     q = zero(rf)
     while !iszero(rf)
         ltf = leading_term(rf)
@@ -183,7 +187,7 @@ function div_multiple(f::_APL, g::_APL, mf::MA.MutableTrait = MA.IsNotMutable())
         end
         qt = div_multiple(ltf, lt)
         q = MA.add!!(q, qt)
-        rf = MA.operate!!(remove_leading_term, rf)
+        rf = MA.operate!!(SA.remove_leading_term, rf)
         rf = MA.operate!!(MA.sub_mul, rf, qt, rg)
     end
     return q
@@ -197,7 +201,7 @@ function Base.rem(f::_APL, g::Union{_APL,AbstractVector{<:_APL}}; kwargs...)
 end
 
 """
-    pseudo_divrem(f::_APL{S}, g::_APL{T}, algo) where {S,T}
+    pseudo_divrem(f::_APL, g::_APL, algo)
 
 Return the pseudo divisor and remainder of `f` modulo `g` as defined in [Knu14, Algorithm R, p. 425].
 
@@ -216,7 +220,9 @@ See also [`pseudo_rem`](@ref).
 *Art of computer programming, volume 2: Seminumerical algorithms.*
 Addison-Wesley Professional. Third edition.
 """
-function pseudo_divrem(f::_APL{S}, g::_APL{T}, algo) where {S,T}
+function pseudo_divrem(f::_APL, g::_APL, algo)
+    S = coefficient_type(f)
+    T = coefficient_type(g)
     return _pseudo_divrem(
         algebraic_structure(MA.promote_operation(-, S, T)),
         f,
@@ -232,19 +238,18 @@ end
 
 function _pseudo_divrem(::UFD, f::_APL, g::_APL, algo)
     ltg = leading_term(g)
-    rg = remove_leading_term(g)
+    rg = SA.remove_leading_term(g)
     ltf = leading_term(f)
     if iszero(f) || !divides(monomial(ltg), ltf)
         return one(f), zero(f), zero(f)
     else
         st = constant_term(coefficient(ltg), f)
-        new_f = st * remove_leading_term(f)
+        new_f = st * SA.remove_leading_term(f)
         qt = term(coefficient(ltf), div_multiple(monomial(ltf), monomial(ltg)))
         new_g = qt * rg
+        R = polynomial_type(f)
         # Check with `::` that we don't have any type unstability on this variable.
-        return convert(typeof(f), st),
-        convert(typeof(f), qt),
-        (new_f - new_g)::typeof(f)
+        return convert(R, st), convert(R, qt), (new_f - new_g)::R
     end
 end
 
@@ -260,7 +265,7 @@ See [`pseudo_divrem`](@ref) for more details.
 Addison-Wesley Professional. Third edition.
 """
 function pseudo_rem(f::_APL, g::_APL, algo)
-    return MA.operate!!(pseudo_rem, MA.mutable_copy(f), g, algo)
+    return MA.operate(pseudo_rem, f, g, algo)
 end
 
 function MA.promote_operation(
@@ -268,11 +273,16 @@ function MA.promote_operation(
     ::Type{P},
     ::Type{Q},
     ::Type{A},
-) where {T,S,P<:_APL{T},Q<:_APL{S},A}
+) where {P<:_APL,Q<:_APL,A}
+    T = coefficient_type(P)
+    S = coefficient_type(Q)
     U1 = MA.promote_operation(*, S, T)
     U2 = MA.promote_operation(*, T, S)
-    # `promote_type(P, Q)` is needed for TypedPolynomials in case they use different variables
-    return polynomial_type(promote_type(P, Q), MA.promote_operation(-, U1, U2))
+    return MA.promote_operation(
+        -,
+        polynomial_type(P, U1),
+        polynomial_type(Q, U2),
+    )
 end
 
 function MA.buffer_for(::typeof(pseudo_rem), F::Type, G::Type, ::Type)
@@ -311,14 +321,15 @@ function MA.buffered_operate!(
     g::_APL,
     algo,
 )
+    SA._assert_same_basis(op, f, g)
     ltg = leading_term(g)
     ltf = leading_term(f)
     # This only makes sense in the univariate case but it's only used for univariate gcd anyway
     skipped_divisions = maxdegree(f) - maxdegree(g) + 1
-    MA.operate!(remove_leading_term, g)
+    MA.operate!(SA.remove_leading_term, g)
     while !iszero(f)
         if isapproxzero(ltf) # TODO `, kwargs...)`
-            MA.operate!(remove_leading_term, f)
+            MA.operate!(SA.remove_leading_term, f)
         elseif !divides(monomial(ltg), ltf)
             # Since the monomials are sorted in decreasing order,
             # lm is larger than all of them hence it cannot divide any of them
@@ -327,9 +338,9 @@ function MA.buffered_operate!(
             if monomial(ltg) > monomial(ltf)
                 break
             end
-            MA.operate!(remove_leading_term, f)
+            MA.operate!(SA.remove_leading_term, f)
         else
-            MA.operate!(remove_leading_term, f)
+            MA.operate!(SA.remove_leading_term, f)
             t = _prepare_s_poly!(op, f, ltf, ltg)
             skipped_divisions -= 1
             MA.buffered_operate!(buffer, MA.sub_mul, f, t, g)
@@ -354,18 +365,50 @@ end
 If the coefficient type is a field, return `rem`, otherwise, return [`pseudo_rem`](@ref).
 """
 function rem_or_pseudo_rem(f::_APL, g::_APL, algo)
-    return MA.operate!!(rem_or_pseudo_rem, MA.mutable_copy(f), g, algo)
+    return MA.operate(rem_or_pseudo_rem, f, g, algo)
+end
+
+function MA.operate(
+    op::Union{typeof(rem),typeof(pseudo_rem),typeof(rem_or_pseudo_rem)},
+    f::_APL,
+    g::_APL,
+    algo,
+)
+    R = MA.promote_operation(op, typeof(f), typeof(g), typeof(algo))
+    rf, g = SA.promote_bases(
+        convert(R, MA.copy_if_mutable(f)),
+        convert(polynomial_type(g), g),
+    )
+    return MA.operate!(op, rf, g, algo)
+end
+
+function MA.mutability(
+    output::AbstractPolynomial,
+    op::Union{typeof(rem),typeof(pseudo_rem),typeof(rem_or_pseudo_rem)},
+    f::AbstractPolynomial,
+    g::AbstractPolynomial,
+    algo,
+)
+    # The kernel temporarily removes and restores the divisor's leading term.
+    # Aliasing its storage requires an independent dividend to preserve `g`.
+    if SA.coeffs(output) === SA.coeffs(g)
+        return MA.IsNotMutable()
+    end
+    return MA.mutability(
+        typeof(output),
+        op,
+        typeof(f),
+        typeof(g),
+        typeof(algo),
+    )
 end
 
 _op(::Field) = rem
 _op(::UFD) = pseudo_rem
 
-function MA.operate!(
-    ::typeof(rem_or_pseudo_rem),
-    f::_APL{S},
-    g::_APL{T},
-    algo,
-) where {S,T}
+function MA.operate!(::typeof(rem_or_pseudo_rem), f::_APL, g::_APL, algo)
+    S = coefficient_type(f)
+    T = coefficient_type(g)
     return MA.operate!(
         _op(algebraic_structure(MA.promote_operation(-, S, T))),
         f,
@@ -377,10 +420,12 @@ end
 function MA.buffered_operate!(
     buffer,
     ::typeof(rem_or_pseudo_rem),
-    f::_APL{S},
-    g::_APL{T},
+    f::_APL,
+    g::_APL,
     algo,
-) where {S,T}
+)
+    S = coefficient_type(f)
+    T = coefficient_type(g)
     return MA.buffered_operate!(
         buffer,
         _op(algebraic_structure(MA.promote_operation(-, S, T))),
@@ -392,10 +437,12 @@ end
 
 function MA.buffer_for(
     ::typeof(rem_or_pseudo_rem),
-    F::Type{<:_APL{S}},
-    G::Type{<:_APL{T}},
+    F::Type{<:_APL},
+    G::Type{<:_APL},
     A::Type,
-) where {S,T}
+)
+    S = coefficient_type(F)
+    T = coefficient_type(G)
     return MA.buffer_for(
         _op(algebraic_structure(MA.promote_operation(-, S, T))),
         F,
@@ -409,7 +456,9 @@ function MA.promote_operation(
     ::Type{P},
     ::Type{Q},
     ::Type{A},
-) where {T,S,P<:_APL{T},Q<:_APL{S},A}
+) where {P<:_APL,Q<:_APL,A}
+    T = coefficient_type(P)
+    S = coefficient_type(Q)
     return _promote_operation_rem_or_pseudo_rem(
         algebraic_structure(MA.promote_operation(-, S, T)),
         P,
@@ -446,29 +495,31 @@ function MA.promote_operation(
     ::Union{typeof(div),typeof(rem)},
     ::Type{P},
     ::Type{Q},
-) where {T,S,P<:_APL{T},Q<:_APL{S}}
+) where {P<:_APL,Q<:_APL}
+    T = coefficient_type(P)
+    S = coefficient_type(Q)
     U = MA.promote_operation(/, promote_to_field(T), promote_to_field(S))
-    # `promote_type(P, Q)` is needed for TypedPolynomials in case they use different variables
-    return polynomial_type(promote_type(P, Q), MA.promote_operation(-, U, U))
+    return MA.promote_operation(-, polynomial_type(P, U), polynomial_type(Q, U))
 end
-function Base.divrem(f::_APL{T}, g::_APL{S}; kwargs...) where {T,S}
-    rf = convert(
-        MA.promote_operation(div, typeof(f), typeof(g)),
-        MA.copy_if_mutable(f),
+function Base.divrem(f::_APL, g::_APL; kwargs...)
+    R = MA.promote_operation(div, typeof(f), typeof(g))
+    rf, g = SA.promote_bases(
+        convert(R, MA.copy_if_mutable(f)),
+        convert(polynomial_type(g), g),
     )
     q = zero(rf)
     r = zero(rf)
     lt = leading_term(g)
-    rg = remove_leading_term(g)
+    rg = SA.remove_leading_term(g)
     lm = monomial(lt)
     while !iszero(rf)
         ltf = leading_term(rf)
         if isapproxzero(ltf; kwargs...)
-            rf = MA.operate!!(remove_leading_term, rf)
+            rf = MA.operate!!(SA.remove_leading_term, rf)
         elseif divides(lm, ltf)
             qt = div_multiple(ltf, lt)
             q = MA.add!!(q, qt)
-            rf = MA.operate!!(remove_leading_term, rf)
+            rf = MA.operate!!(SA.remove_leading_term, rf)
             rf = MA.operate!!(MA.sub_mul, rf, qt, rg)
         elseif lm > monomial(ltf)
             # Since the monomials are sorted in decreasing order,
@@ -477,33 +528,32 @@ function Base.divrem(f::_APL{T}, g::_APL{S}; kwargs...) where {T,S}
             break
         else
             r = MA.add!!(r, ltf)
-            rf = MA.operate!!(remove_leading_term, rf)
+            rf = MA.operate!!(SA.remove_leading_term, rf)
         end
     end
     return q, r
 end
-function Base.divrem(
-    f::_APL{T},
-    g::AbstractVector{<:_APL{S}};
-    kwargs...,
-) where {T,S}
-    rf = convert(
-        MA.promote_operation(div, typeof(f), eltype(g)),
-        MA.copy_if_mutable(f),
-    )
+function Base.divrem(f::_APL, g::AbstractVector{<:_APL}; kwargs...)
+    R = MA.promote_operation(div, typeof(f), eltype(g))
+    rf = convert(R, MA.copy_if_mutable(f))
+    g = map(p -> convert(polynomial_type(p), p), g)
+    for p in g
+        rf, _ = SA.promote_bases(rf, p)
+    end
+    g = map(p -> first(SA.promote_bases(p, rf)), g)
     r = zero(rf)
     q = similar(g, typeof(rf))
     for i in eachindex(q)
         q[i] = zero(rf)
     end
     lt = leading_term.(g)
-    rg = remove_leading_term.(g)
+    rg = SA.remove_leading_term.(g)
     lm = monomial.(lt)
     useful = BitSet(eachindex(g))
     while !iszero(rf)
         ltf = leading_term(rf)
         if isapproxzero(ltf; kwargs...)
-            rf = MA.operate!!(remove_leading_term, rf)
+            rf = MA.operate!!(SA.remove_leading_term, rf)
             continue
         end
         divisionoccured = false
@@ -511,7 +561,7 @@ function Base.divrem(
             if divides(lm[i], ltf)
                 qt = div_multiple(ltf, lt[i])
                 q[i] = MA.add!!(q[i], qt)
-                rf = MA.operate!!(remove_leading_term, rf)
+                rf = MA.operate!!(SA.remove_leading_term, rf)
                 rf = MA.operate!!(MA.sub_mul, rf, qt, rg[i])
                 divisionoccured = true
                 break
@@ -527,7 +577,7 @@ function Base.divrem(
                 break
             else
                 r = MA.add!!(r, ltf)
-                rf = MA.operate!!(remove_leading_term, rf)
+                rf = MA.operate!!(SA.remove_leading_term, rf)
             end
         end
     end
diff --git a/src/gcd.jl b/src/gcd.jl
index 33d60206..5799b5b7 100644
--- a/src/gcd.jl
+++ b/src/gcd.jl
@@ -128,8 +128,39 @@ function MA.promote_operation(
     return MA.promote_operation(rem_or_pseudo_rem, P, Q, A)
 end
 
+function MA.promote_operation(
+    ::typeof(gcd),
+    P::Type{<:AbstractMonomialLike},
+    Q::Type{<:AbstractMonomialLike},
+)
+    return promote_type(monomial_type(P), monomial_type(Q))
+end
+
+function MA.promote_operation(
+    ::typeof(gcd),
+    P::Type{<:AbstractTermLike},
+    Q::Type{<:AbstractTermLike},
+    ::Type{<:AbstractUnivariateGCDAlgorithm} = SubresultantAlgorithm,
+)
+    M = MA.promote_operation(gcd, monomial_type(P), monomial_type(Q))
+    T = MA.promote_operation(
+        _coefficient_gcd,
+        coefficient_type(P),
+        coefficient_type(Q),
+    )
+    return term_type(M, T)
+end
+
+function MA.promote_operation(
+    ::typeof(_coefficient_gcd),
+    P::Type{<:_APL},
+    Q::Type{<:_APL},
+)
+    return MA.promote_operation(gcd, P, Q)
+end
+
 """
-    function gcd(p1::_APL{T}, p2::_APL{S}) where {T, S}
+    function gcd(p1::_APL, p2::_APL)
 
 Returns a greatest common divisor of `p1` and `p2`. Note that it does not make
 sense, in general, to speak of "the" greatest common divisor of u and v; there
@@ -174,12 +205,12 @@ This is the [`GeneralizedEuclideanAlgorithm`](@ref).
 Addison-Wesley Professional. Third edition.
 """
 function Base.gcd(
-    p1::_APL{T},
-    p2::_APL{S},
+    p1::_APL,
+    p2::_APL,
     algo::AbstractUnivariateGCDAlgorithm = SubresultantAlgorithm(),
     m1::MA.MutableTrait = MA.IsNotMutable(),
     m2::MA.MutableTrait = MA.IsNotMutable(),
-) where {T,S}
+)
     # If one of these is zero, `shift` should be infinite
     # for this method to work so we exclude these cases.
     if isapproxzero(p1)
@@ -213,12 +244,12 @@ function Base.gcd(
 end
 
 function Base.gcd(
-    t1::AbstractTermLike{T},
-    t2::AbstractTermLike{S},
+    t1::AbstractTermLike,
+    t2::AbstractTermLike,
     ::AbstractUnivariateGCDAlgorithm = SubresultantAlgorithm(),
     m1::MA.MutableTrait = MA.IsNotMutable(),
     m2::MA.MutableTrait = MA.IsNotMutable(),
-) where {T,S}
+)
     return term(
         _coefficient_gcd(coefficient(t1), coefficient(t2)),
         gcd(monomial(t1), monomial(t2)),
@@ -323,12 +354,12 @@ function MA.operate(
 end
 
 function deflated_gcd(
-    p1::_APL{T},
-    p2::_APL{S},
+    p1::_APL,
+    p2::_APL,
     algo,
     m1::MA.MutableTrait,
     m2::MA.MutableTrait,
-) where {T,S}
+)
     i1, i2, num_common = _extracted_variable(p1, p2)
     if iszero(i1)
         if iszero(i2)
@@ -370,12 +401,12 @@ function deflated_gcd(
 end
 
 function Base.gcdx(
-    p1::_APL{T},
-    p2::_APL{S},
+    p1::_APL,
+    p2::_APL,
     algo::AbstractUnivariateGCDAlgorithm = GeneralizedEuclideanAlgorithm(),
-) where {T,S}
+)
     i1, i2, num_common = _extracted_variable(p1, p2)
-    R = MA.promote_operation(gcd, typeof(p1), typeof(p2))
+    R = polynomial_type(MA.promote_operation(gcd, typeof(p1), typeof(p2)))
     if iszero(i1)
         if iszero(i2)
             return univariate_gcdx(p1, p2, algo)
@@ -568,7 +599,7 @@ function primitive_univariate_gcd!(
     if maxdegree(p) < maxdegree(q)
         return primitive_univariate_gcd!(q, p, algo)
     end
-    R = MA.promote_operation(gcd, typeof(p), typeof(q))
+    R = polynomial_type(MA.promote_operation(gcd, typeof(p), typeof(q)))
     u = convert(R, p)
     v = convert(R, q)
     while true
@@ -611,7 +642,7 @@ function primitive_univariate_gcd!(
     if maxdegree(p) < maxdegree(q)
         return primitive_univariate_gcd!(q, p, algo)
     end
-    R = MA.promote_operation(gcd, typeof(p), typeof(q))
+    R = polynomial_type(MA.promote_operation(gcd, typeof(p), typeof(q)))
     u = convert(R, p)
     v = convert(R, q)
     if isapproxzero(v)
@@ -695,7 +726,7 @@ function primitive_univariate_gcdx(
         a, b, g = primitive_univariate_gcdx(v0, u0, algo)
         return b, a, g
     end
-    R = MA.promote_operation(gcd, typeof(u0), typeof(v0))
+    R = polynomial_type(MA.promote_operation(gcd, typeof(u0), typeof(v0)))
     u = convert(R, u0)
     v = convert(R, v0)
     if isapproxzero(v)
@@ -753,12 +784,14 @@ If the coefficients are not `AbstractFloat`, this
 Addison-Wesley Professional. Third edition.
 """
 function univariate_gcd(
-    p1::_APL{S},
-    p2::_APL{T},
+    p1::_APL,
+    p2::_APL,
     algo::AbstractUnivariateGCDAlgorithm,
     m1::MA.MutableTrait,
     m2::MA.MutableTrait,
-) where {S,T}
+)
+    S = coefficient_type(p1)
+    T = coefficient_type(p2)
     return univariate_gcd(
         _field_absorb(algebraic_structure(S), algebraic_structure(T)),
         p1,
@@ -799,10 +832,12 @@ function univariate_gcd(
 end
 
 function univariate_gcdx(
-    p1::_APL{S},
-    p2::_APL{T},
+    p1::_APL,
+    p2::_APL,
     algo::AbstractUnivariateGCDAlgorithm,
-) where {S,T}
+)
+    S = coefficient_type(p1)
+    T = coefficient_type(p2)
     return univariate_gcdx(
         _field_absorb(algebraic_structure(S), algebraic_structure(T)),
         p1,
@@ -865,7 +900,7 @@ function termwise_content(p::_APL, algo, mutability::MA.MutableTrait)
 end
 
 """
-    content(poly::_APL{T}, algo::AbstractUnivariateGCDAlgorithm, mutability::MA.MutableTrait) where {T}
+    content(poly::_APL, algo::AbstractUnivariateGCDAlgorithm, mutability::MA.MutableTrait)
 
 Return the *content* of the polynomial `poly` over a unique factorization
 domain `S` as defined in [Knu14, (3) p. 423].
@@ -880,10 +915,14 @@ The output can be mutated without affecting `poly` if `mutability` is
 Addison-Wesley Professional. Third edition.
 """
 function content(
-    poly::_APL{T},
+    poly::_APL,
     algo::AbstractUnivariateGCDAlgorithm,
     mutability::MA.MutableTrait,
-) where {T}
+)
+    T = coefficient_type(poly)
+    if T <: AbstractFloat
+        return one(T)
+    end
     P = MA.promote_operation(gcd, T, T)
     coefs = coefficients(poly)
     if isempty(coefs)
@@ -926,16 +965,9 @@ function content(
     end
     return g::P
 end
-function content(
-    ::_APL{T},
-    ::AbstractUnivariateGCDAlgorithm,
-    ::MA.MutableTrait,
-) where {T<:AbstractFloat}
-    return one(T)
-end
 
 """
-    primitive_part(poly::_APL{T}, algo::AbstractUnivariateGCDAlgorithm) where {T}
+    primitive_part(poly::_APL, algo::AbstractUnivariateGCDAlgorithm)
 
 Return the *primitive part* of the polynomial `poly` over a unique
 factorization domain `S` as defined in [Knu14, (3) p. 423].
@@ -952,25 +984,14 @@ function primitive_part(
     algo::AbstractUnivariateGCDAlgorithm,
     mutability::MA.MutableTrait,
 )
+    if coefficient_type(p) <: Union{AbstractFloat,Complex{<:AbstractFloat}}
+        return p
+    end
     return primitive_part_content(p, algo, mutability)[1]
 end
-function primitive_part(
-    p::_APL{<:AbstractFloat},
-    ::AbstractUnivariateGCDAlgorithm,
-    ::MA.MutableTrait,
-)
-    return p
-end
-function primitive_part(
-    p::_APL{<:Complex{<:AbstractFloat}},
-    ::AbstractUnivariateGCDAlgorithm,
-    ::MA.MutableTrait,
-)
-    return p
-end
 
 """
-    primitive_part_content(poly::_APL{T}, algo::AbstractUnivariateGCDAlgorithm) where {T}
+    primitive_part_content(poly::_APL, algo::AbstractUnivariateGCDAlgorithm)
 
 Return the *primitive part* and *content* of the polynomial `poly` over a unique
 factorization domain `S` as defined in [Knu14, (3) p. 423]. This is more
diff --git a/src/hash.jl b/src/hash.jl
index 1eb7d127..046b278a 100644
--- a/src/hash.jl
+++ b/src/hash.jl
@@ -15,7 +15,12 @@ function Base.hash(m::AbstractMonomial, u::UInt)
     end
 end
 
-# hash for SA.Term is defined in StarAlgebras
+function Base.hash(t::AbstractTerm, u::UInt)
+    if !iszero(t) && isconstant(t)
+        return hash(coefficient(t), u)
+    end
+    return invoke(hash, Tuple{SA.Term,UInt}, t, u)
+end
 
 function Base.hash(p::AbstractPolynomial, u::UInt)
     if iszero(p)
diff --git a/src/mb_algebra.jl b/src/mb_algebra.jl
new file mode 100644
index 00000000..519a5a4a
--- /dev/null
+++ b/src/mb_algebra.jl
@@ -0,0 +1,35 @@
+# From MultivariateBases/src/MultivariateBases.jl — algebra() function
+
+function algebra(
+    basis::Union{FullBasis{B,V},SubBasis{B,V}},
+) where {B,V}
+    mstruct = if is_commutative(V)
+        MStruct(basis)
+    else
+        if B != Monomial
+            error(
+                "Only the `Monomial` basis is supported with noncommutative variables",
+            )
+        end
+        SA.DiracMStructure(basis, *)
+    end
+    return SA.StarAlgebra(Variables{B}(variables(basis)), mstruct)
+end
+
+function MA.promote_operation(
+    ::typeof(algebra),
+    BT::Type{
+        <:Union{
+            FullBasis{B,V,E},
+            SubBasis{B,V,E},
+        },
+    },
+) where {B,V,E}
+    P = Polynomial{B,V,E}
+    MS = if is_commutative(V)
+        MStruct{B,V,E,SA.key_type(BT),BT}
+    else
+        SA.DiracMStructure{P,SA.key_type(BT),BT,typeof(*)}
+    end
+    return SA.StarAlgebra{Variables{B,V},P,MS}
+end
diff --git a/src/mb_arithmetic.jl b/src/mb_arithmetic.jl
new file mode 100644
index 00000000..bbac18d5
--- /dev/null
+++ b/src/mb_arithmetic.jl
@@ -0,0 +1,58 @@
+# SA.Term constructor for Polynomial{...} basis elements (monomials)
+function SA.Term(coeff, p::Polynomial{B}) where {B}
+    alg = algebra(FullBasis{B}(p))
+    idx = exponents(p)
+    return SA.Term(alg, idx, coeff)
+end
+
+function Base.zero(
+    ::Type{SA.AlgebraElement{T,A,C}},
+) where {T,A<:_PolynomialAlgebra,C}
+    M = monomial_type(A)
+    return zero(algebra_element(term(zero(T), constant_monomial(M))))
+end
+
+# Arithmetic for AlgebraElement polynomials.
+# Since polynomials ARE AlgebraElements now, SA handles most arithmetic.
+
+const _AE = AbstractPolynomial
+
+# polynomial() on an AlgebraElement is the identity
+polynomial(a::AbstractPolynomial) = a
+polynomial(a::AbstractTermLike) = algebra_element(a)
+function polynomial(a::AbstractPolynomialLike, ::Type{T}) where {T}
+    return SA.AlgebraElement{T}(polynomial(a))
+end
+
+# _APL (monomials, terms) + _AE: convert to AE first.
+# Note: SA already defines +(::Term, ::AlgebraElement); we only handle
+# AbstractMonomialLike (not Term) to avoid ambiguity.
+for op in [:+, :-, :*]
+    @eval begin
+        Base.$op(p::AbstractMonomialLike, q::_AE) = $op(algebra_element(p), q)
+        Base.$op(p::_AE, q::AbstractMonomialLike) = $op(p, algebra_element(q))
+        Base.$op(p::Polynomial{<:AbstractMonomialIndexed}, q::_AE) =
+            $op(algebra_element(SA.Term(1, p)), q)
+        Base.$op(p::_AE, q::Polynomial{<:AbstractMonomialIndexed}) =
+            $op(p, algebra_element(SA.Term(1, q)))
+    end
+end
+
+# Scalars use the polynomial's coefficient type.
+for op in [:+, :-]
+    @eval begin
+        function Base.$op(p::Union{T,Number}, q::_AE{T}) where {T}
+            i = implicit(q)
+            return $op(constant_algebra_element(SA.basis(i), convert(T, p)), i)
+        end
+        function Base.$op(p::_AE{T}, q::Union{T,Number}) where {T}
+            i = implicit(p)
+            return $op(i, constant_algebra_element(SA.basis(i), convert(T, q)))
+        end
+    end
+end
+
+function term_element(α, p::Polynomial)
+    return SA.algebra_element(SA.Term(α, p))
+end
+Base.:*(α::Number, p::Polynomial) = term_element(α, p)
diff --git a/src/mb_bases.jl b/src/mb_bases.jl
new file mode 100644
index 00000000..208e3759
--- /dev/null
+++ b/src/mb_bases.jl
@@ -0,0 +1,282 @@
+# From MultivariateBases/src/bases.jl
+
+const FullBasis{B,V,E} = SA.MappedBasis{Polynomial{B,V,E}}
+const SubBasis{B,V,E} = SA.SubBasis{Polynomial{B,V,E}}
+const MonomialIndexedBasis{B,V,E} = Union{SubBasis{B,V,E},FullBasis{B,V,E}}
+
+exponents_index(::FullBasis, exp) = exp
+exponents_index(b::SubBasis, exp) = SA.key_index(b, exp)
+
+function ordering(::Type{MonomialIndexedBasis{B,V,E}}) where {B,V,E}
+    return ordering(E)
+end
+ordering(b::MonomialIndexedBasis) = ordering(typeof(b))
+
+struct _ExponentsWithVariables{V}
+    variables::V
+end
+
+(e::_ExponentsWithVariables)(mono) = exponents(mono, e.variables)
+
+function Base.:(==)(a::_ExponentsWithVariables, b::_ExponentsWithVariables)
+    return unique(sort(a.variables)) == unique(sort(b.variables))
+end
+
+function FullBasis(vars::Variables{B,V}) where {B,V}
+    O = ordering(vars.variables)
+    exps = ExponentsIterator{O}(constant_monomial_exponents(vars))
+    inverse_map = if is_commutative(V)
+        exponents
+    else
+        _ExponentsWithVariables(vars.variables)
+    end
+    return SA.MappedBasis{Polynomial{B,V,eltype(exps)}}(exps, vars, inverse_map)
+end
+
+function FullBasis{B}(vars) where {B}
+    return FullBasis(Variables{B,typeof(vars)}(vars))
+end
+
+function FullBasis{B}(p::_APL) where {B}
+    return FullBasis{B}(variables(p))
+end
+# Also accept Polynomial basis elements (monomials in the new world)
+function FullBasis{B}(p::Polynomial) where {B}
+    return FullBasis{B}(variables(p))
+end
+
+function SA.star(b::MonomialIndexedBasis{B,V,E}, exp::E) where {B,V,E}
+    return b[SA.star(b[exp])]
+end
+
+constant_monomial_exponents(b::FullBasis) = constant_monomial_exponents(b.map)
+
+monomial_type(::Type{<:FullBasis{B,V}}) where {B,V} = monomial_type(V)
+function polynomial_type(basis::FullBasis, ::Type{T}) where {T}
+    return polynomial_type(typeof(basis), T)
+end
+
+nvariables(basis::FullBasis) = nvariables(basis.map)
+nvariables(basis::SubBasis) = nvariables(parent(basis))
+variables(basis::FullBasis) = variables(basis.map)
+variables(basis::SubBasis) = variables(parent(basis))
+
+function explicit_basis_covering(
+    full::FullBasis{B},
+    target::SubBasis{B},
+) where {B}
+    return SA.SubBasis(full, target.keys)
+end
+
+function monomial_type(::Type{<:SA.SubBasis{T,I,K,B}}) where {T,I,K,B}
+    return monomial_type(B)
+end
+
+# FIXME workaround for TP, we should redirect to typeof
+monomial_type(b::FullBasis) = monomial_type(b.map)
+monomial_type(b::SubBasis) = monomial_type(parent(b))
+
+# FIXME type piracy
+SA.comparable(::ExponentsIterator{M}) where {M} = M()
+
+function unsafe_basis(
+    ::Type{B},
+    monomials_vec::AbstractVector{<:AbstractMonomial},
+) where {B<:AbstractMonomialIndexed}
+    return SA.SubBasis(
+        FullBasis{B}(variables(monomials_vec)),
+        exponents.(monomials_vec),
+    )
+end
+
+function Base.getindex(
+    ::FullBasis{B},
+    monomials_vec::AbstractVector{M},
+) where {B,M<:AbstractMonomial}
+    return unsafe_basis(B, monomial_vector(monomials_vec)::AbstractVector{M})
+end
+
+function SubBasis{B}(
+    monomials_vec::AbstractVector{<:AbstractTermLike},
+) where {B<:AbstractMonomialIndexed}
+    return unsafe_basis(
+        B,
+        monomial_vector(monomials_vec)::AbstractVector{<:AbstractMonomial},
+    )
+end
+
+SubBasis{B}(monos::Tuple) where {B} = SubBasis{B}([monos...])
+
+implicit_basis(basis::SubBasis) = parent(basis)
+implicit_basis(basis::FullBasis) = basis
+
+function implicit(a::SA.AlgebraElement)
+    basis = implicit_basis(SA.basis(a))
+    return algebra_element(SA.coeffs(a, basis), basis)
+end
+
+function _comparable_type(
+    ::Type{<:SA.MappedBasis{T,I,<:ExponentsIterator{M}}},
+) where {T,I,M}
+    return M
+end
+function _coeffs_type(::Type{C}, ::Type{B}) where {C,B<:FullBasis}
+    T = eltype(C) # Even works for `NTuple`!
+    E = SA.key_type(B)
+    L = _comparable_type(B)
+    return SA.SparseCoefficients{E,T,_similar_type(C, E),C,L}
+end
+
+function _coeffs_type(::Type{C}, ::Type{B}) where {C,B<:SubBasis}
+    return C
+end
+
+function algebra_element_type(
+    ::Type{C},
+    ::Type{B},
+) where {C,B}
+    A = MA.promote_operation(algebra, B)
+    return SA.AlgebraElement{
+        eltype(C),
+        A,
+        _coeffs_type(C, B),
+    }
+end
+
+function MA.promote_operation(
+    ::typeof(implicit),
+    ::Type{AE},
+) where {AG,T,AE<:SA.AlgebraElement{T,AG}}
+    BT =
+        MA.promote_operation(implicit_basis, MA.promote_operation(SA.basis, AE))
+    return algebra_element_type(Vector{T}, BT)
+end
+
+MA.promote_operation(::typeof(implicit_basis), B::Type{<:SA.ImplicitBasis}) = B
+function MA.promote_operation(
+    ::typeof(implicit_basis),
+    ::Type{<:SA.SubBasis{T,I,K,B}},
+) where {T,I,K,B}
+    return B
+end
+
+function _explicit_basis(coeffs, basis::FullBasis{B}) where {B}
+    return SA.SubBasis(basis, _lazy_collect(SA.keys(coeffs)))
+end
+
+_explicit_basis(_, basis::SubBasis) = basis
+
+function explicit_basis(p::AbstractPolynomialLike)
+    return SubBasis{Monomial}(monomials(p))
+end
+
+function explicit_basis(a::SA.AlgebraElement)
+    return _explicit_basis(SA.coeffs(a), SA.basis(a))
+end
+
+function explicit_basis_type(BT::Type{<:FullBasis{B,V,E}}) where {B,V,E}
+    return SA.SubBasis{eltype(BT),Int,SA.key_type(BT),BT,Vector{E}}
+end
+
+function empty_basis(
+    ::Type{<:SubBasis{B,M}},
+) where {B<:AbstractMonomialIndexed,M}
+    return unsafe_basis(B, empty_monomial_vector(M))
+end
+
+function maxdegree_basis(
+    basis::FullBasis{B},
+    maxdeg::Int,
+) where {B<:AbstractMonomialIndexed}
+    return unsafe_basis(B, monomials(variables(basis), 0:maxdeg))
+end
+
+variables(c::SA.AbstractCoefficients) = variables(SA.keys(c))
+
+_lazy_collect(v::AbstractVector) = collect(v)
+_lazy_collect(v::Tuple) = collect(v)
+_lazy_collect(v::Vector) = v
+
+sparse_coefficients(p::AbstractPolynomial) = SA.coeffs(p)
+sparse_coefficients(t::AbstractTermLike) = SA.coeffs(algebra_element(t))
+algebra_element(p::AbstractPolynomial) = p
+algebra_element(t::AbstractTerm) = SA.algebra_element(t)
+algebra_element(m::AbstractMonomialLike) = SA.algebra_element(term(m))
+
+function algebra_element(f::Function, basis::SubBasis)
+    return algebra_element(map(f, eachindex(basis)), basis)
+end
+
+_similar_type(::Type{<:NTuple{N,Any}}, ::Type{T}) where {N,T} = NTuple{N,T}
+function _similar_type(::Type{V}, ::Type{T}) where {V<:AbstractVector,T}
+    return SA.similar_type(V, T)
+end
+
+function full_basis_type(::Type{B}, ::Type{Polynomial{B,V,E}}) where {B,V,E}
+    inverse_map =
+        is_commutative(V) ? typeof(exponents) : _ExponentsWithVariables{V}
+    return SA.MappedBasis{
+        Polynomial{B,V,E},
+        E,
+        ExponentsIterator{ordering(V),Nothing,E},
+        Variables{B,V},
+        inverse_map,
+    }
+end
+
+function MA.promote_operation(
+    ::typeof(algebra_element),
+    P::Type{<:AbstractPolynomialLike},
+)
+    return polynomial_type(P)
+end
+
+_one_if_type(α) = α
+_one_if_type(::Type{T}) where {T} = one(T)
+
+function constant_algebra_element_type(
+    ::Type{BT},
+    ::Type{T},
+) where {BT<:FullBasis,T}
+    return algebra_element_type(Tuple{T}, BT)
+end
+
+function constant_algebra_element(b::FullBasis, α)
+    return SA.algebra_element(
+        SA.Term(algebra(b), constant_monomial_exponents(b), _one_if_type(α)),
+    )
+end
+
+function constant_algebra_element_type(
+    ::Type{B},
+    ::Type{T},
+) where {B<:SubBasis,T}
+    A = MA.promote_operation(algebra, B)
+    return SA.AlgebraElement{T,A,Vector{T}}
+end
+
+function constant_algebra_element(basis::SubBasis, α)
+    return algebra_element(
+        [_one_if_type(α)],
+        SA.SubBasis(
+            parent(basis),
+            [constant_monomial_exponents(parent(basis))],
+        ),
+    )
+end
+
+function SA.promote_with_map(::Variables{B}, vars, map) where {B}
+    return Variables{B}(vars), map
+end
+
+SA.promote_with_map(::FullBasis, vars, m) = FullBasis(vars), m
+
+function promote_variables_with_maps(a::Variables, b::Variables)
+    _a, _b = promote_variables_with_maps(a.variables, b.variables)
+    return SA.maybe_promote(a, _a...), SA.maybe_promote(b, _b...)
+end
+
+function SA.promote_bases_with_maps(a::FullBasis, b::FullBasis)
+    _a, _b = promote_variables_with_maps(a.map, b.map)
+    return SA.maybe_promote(a, _a...), SA.maybe_promote(b, _b...)
+end
diff --git a/src/mb_interface.jl b/src/mb_interface.jl
new file mode 100644
index 00000000..2e6c4dd8
--- /dev/null
+++ b/src/mb_interface.jl
@@ -0,0 +1,17 @@
+# From MultivariateBases/src/interface.jl
+
+"""
+    maxdegree_basis(basis::SA.AbstractBasis, maxdegree::Int)
+
+Return the explicit version of `basis` generating all polynomials of degree up to
+`maxdegree`.
+"""
+function maxdegree_basis end
+
+"""
+    explicit_basis_covering(basis::SA.AbstractBasis, target::SA.ExplicitBasis)
+
+Return the minimal basis of type `B` that can generate all polynomials generated by
+the basis `target`.
+"""
+function explicit_basis_covering end
diff --git a/src/mb_monomial_basis.jl b/src/mb_monomial_basis.jl
new file mode 100644
index 00000000..004e627a
--- /dev/null
+++ b/src/mb_monomial_basis.jl
@@ -0,0 +1,260 @@
+# Polynomial{Monomial,...} acts as a monomial: support term creation and const mult
+term(coeff, p::Polynomial{<:AbstractMonomialBasis}) = SA.Term(coeff, p)
+term(p::Polynomial{<:AbstractMonomialBasis}) = SA.Term(one(Int), p)
+left_constant_mult(α, p::Polynomial{<:AbstractMonomialBasis}) = SA.Term(α, p)
+right_constant_mult(p::Polynomial{<:AbstractMonomialBasis}, α) = left_constant_mult(α, p)
+
+function term_type(::Type{Polynomial{B,V,E}}, ::Type{T}) where {B,V,E,T}
+    A = MA.promote_operation(algebra, full_basis_type(B, Polynomial{B,V,E}))
+    return SA.Term{T,A,E}
+end
+term_type(::Type{P}) where {P<:Polynomial} = term_type(P, Int)
+
+function explicit_basis_covering(
+    full::FullBasis{B},
+    target::SubBasis{<:AbstractMonomialBasis},
+) where {B<:AbstractMonomialBasis}
+    return SA.SubBasis(full, target.keys)
+end
+
+function explicit_basis_covering(
+    full::FullBasis{B},
+    target::SubBasis{B},
+) where {B<:AbstractMonomialBasis}
+    @assert full == parent(target)
+    return SA.SubBasis(parent(target), target.keys)
+end
+
+function Base.adjoint(p::Polynomial{B}) where {B<:AbstractMonomialIndexed}
+    mono = adjoint(monomial(p))
+    return Polynomial(Variables{B}(variables(mono)), exponents(mono))
+end
+
+degree_one_univariate_polynomial(::Type{Monomial}, value) = value
+
+function recurrence_eval(::Type{Monomial}, previous, value, degree)
+    return previous[degree] * value
+end
+
+# /!\ assumes commutative variables
+function (m::MStruct{Monomial,V,E})(a::E, b::E, ::Type{E}) where {V,E}
+    return SA.SparseCoefficients((a .+ b,), (1,))
+end
+
+function SA.term_product_style(
+    ::MStruct{Monomial,V,E,E,<:FullBasis},
+    ::AbstractMonomialOrdering,
+) where {V,E}
+    return SA.OrderedTermProduct()
+end
+
+# Monomial multiplication: align variables via promote_variables_with_maps,
+# then apply f to the aligned exponent vectors.
+function map_exponents(f, a::Polynomial{Monomial}, b::Polynomial{Monomial})
+    a, b = SA.promote_bases(a, b)
+    return Polynomial(a.variables, f.(exponents(a), exponents(b)))
+end
+
+function Base.:*(a::Polynomial{Monomial}, b::Polynomial{Monomial})
+    return map_exponents(+, a, b)
+end
+
+# Monomial power: multiply exponents by n
+function Base.:^(m::Polynomial{Monomial}, n::Integer)
+    return Polynomial(m.variables, m.exponents .* n)
+end
+
+SA.coeffs(p::Polynomial{Monomial}, ::FullBasis{Monomial}) = p
+
+function polynomial_type(
+    ::Union{SubBasis{B,V,E},Type{<:SubBasis{B,V,E}}},
+    ::Type{T},
+) where {B,V,E,T}
+    return _polynomial_type(B, V, T)
+end
+
+function polynomial_type(::Type{P}, ::Type{T}) where {B,P<:Polynomial{B},T}
+    return algebra_element_type(Vector{T}, full_basis_type(B, P))
+end
+
+function keys_as_monomials(keys, mb::FullBasis)
+    return monomial_vector(monomial.(Ref(mb.map.variables), keys))
+end
+keys_as_monomials(mb::SubBasis) = keys_as_monomials(mb.keys, parent(mb))
+
+function polynomial(f::Function, mb::SubBasis{Monomial})
+    return polynomial(f, keys_as_monomials(mb))
+end
+
+function polynomial(Q::AbstractMatrix, mb::SubBasis{Monomial}, T::Type)
+    return polynomial(Q, keys_as_monomials(mb), T)
+end
+
+function coefficients(
+    p::AbstractPolynomialLike,
+    basis::SubBasis{Monomial},
+)
+    return coefficients(p, keys_as_monomials(basis))
+end
+
+function coefficients(p::AbstractPolynomialLike, ::FullBasis{Monomial})
+    return p
+end
+
+function _assert_constant(α) end
+
+function _assert_constant(
+    x::Union{Polynomial,SA.AlgebraElement,AbstractPolynomialLike},
+)
+    return error("Expected constant element, got type `$(typeof(x))`")
+end
+
+mindegree(basis::SubBasis{Monomial}) = minimum(sum, basis.keys)
+maxdegree(basis::SubBasis) = maximum(sum, basis.keys)
+extdegree(basis::SubBasis{Monomial}) = extrema(sum, basis.keys)
+
+variable_index(basis::FullBasis, v) = variable_index(basis.map, v)
+variable_index(basis::SubBasis, v) = variable_index(basis.parent_basis, v)
+
+function mindegree(basis::SubBasis{Monomial}, v)
+    i = variable_index(basis, v)
+    return minimum(Base.Fix2(getindex, i), basis.keys)
+end
+function maxdegree(basis::SubBasis, v)
+    i = variable_index(basis, v)
+    return maximum(Base.Fix2(getindex, i), basis.keys)
+end
+function extdegree(basis::SubBasis{Monomial}, v)
+    i = variable_index(basis, v)
+    return extrema(Base.Fix2(getindex, i), basis.keys)
+end
+
+_promote_coef(::Type{T}, ::Type{Monomial}) where {T} = T
+
+function polynomial_type(::Type{<:FullBasis{B,V}}, ::Type{T}) where {T,B,V}
+    return _polynomial_type(B, V, T)
+end
+
+function _polynomial_type(::Type{B}, ::Type{V}, ::Type{T}) where {B,V,T}
+    return polynomial_type(monomial_type(V), _promote_coef(T, B))
+end
+
+_vec(v::Vector) = v
+_vec(v::AbstractVector) = collect(v)
+
+function SA.coeffs(
+    cfs,
+    source::MonomialIndexedBasis{B1},
+    target::MonomialIndexedBasis{B2},
+) where {B1,B2}
+    source === target && return cfs
+    source == target && return cfs
+    if B1 === B2 && target isa FullBasis
+        return SA.SparseCoefficients(_vec(source.keys), _vec(cfs))
+    else
+        res = SA.zero_coeffs(
+            _promote_coef(_promote_coef(SA.value_type(cfs), B1), B2),
+            target,
+        )
+        return SA.coeffs!(res, cfs, source, target)
+    end
+end
+
+function SA.adjoint_coeffs(
+    cfs,
+    source::MonomialIndexedBasis{B1},
+    target::MonomialIndexedBasis{B2},
+) where {B1,B2}
+    source === target && return cfs
+    source == target && return cfs
+    res = SA.zero_coeffs(
+        _promote_coef(_promote_coef(SA.value_type(cfs), B1), B2),
+        source,
+    )
+    return SA.adjoint_coeffs!(res, cfs, source, target)
+end
+
+function SA.adjoint_coeffs!(
+    res,
+    cfs,
+    source::MonomialIndexedBasis{B1},
+    target::MonomialIndexedBasis{B2},
+) where {B1,B2}
+    if B1 === B2
+        return @invoke SA.adjoint_coeffs!(
+            res,
+            cfs,
+            source::SA.AbstractBasis,
+            target::SA.AbstractBasis,
+        )
+    end
+    MA.operate!(zero, res)
+    target_full = implicit_basis(target)
+    for (i, src_elem) in enumerate(source)
+        full_col = SA.coeffs(algebra_element(src_elem), target_full)
+        col = SA.coeffs(full_col, target_full, target)
+        for (j, v) in SA.nonzero_pairs(col)
+            res[i] += v * cfs[j]
+        end
+    end
+    return res
+end
+
+# FIXME this assumes that the basis is invariant under adjoint
+SA.star(::SubBasis, coeffs) = SA.star.(coeffs)
+
+function _show_vector(io::IO, mime::MIME, v, map = identity)
+    print(io, '[')
+    first = true
+    for el in v
+        if !first
+            print(io, ", ")
+        end
+        first = false
+        show(io, mime, map(el))
+    end
+    return print(io, ']')
+end
+
+function _show(io::IO, mime::MIME, basis::SubBasis{B}) where {B}
+    print(io, "SubBasis{$(nameof(B))}(")
+    _show_vector(
+        io,
+        mime,
+        basis.keys,
+        Base.Fix1(monomial, variables(basis)),
+    )
+    print(io, ')')
+    return
+end
+
+function _show(io::IO, mime::MIME, basis::FullBasis{B}) where {B}
+    print(io, "FullBasis{$(nameof(B))}(")
+    _show_vector(io, mime, variables(basis))
+    print(io, ")")
+    return
+end
+
+function Base.show(
+    io::IO,
+    mime::MIME"text/plain",
+    basis::Union{Variables,SubBasis,FullBasis},
+)
+    return _show(io, mime, basis)
+end
+
+function Base.show(
+    io::IO,
+    mime::MIME"text/print",
+    basis::Union{Variables,SubBasis,FullBasis},
+)
+    return _show(io, mime, basis)
+end
+
+function Base.print(io::IO, basis::Union{Variables,SubBasis,FullBasis})
+    return show(io, MIME"text/print"(), basis)
+end
+
+function Base.show(io::IO, basis::Union{Variables,SubBasis,FullBasis})
+    return show(io, MIME"text/plain"(), basis)
+end
diff --git a/src/mb_mstructures.jl b/src/mb_mstructures.jl
new file mode 100644
index 00000000..2f898bd5
--- /dev/null
+++ b/src/mb_mstructures.jl
@@ -0,0 +1,58 @@
+# From MultivariateBases/src/mstructures.jl
+
+struct MStruct{
+    B<:AbstractMonomialIndexed,
+    V,
+    E,
+    I,
+    BT<:SA.AbstractBasis{Polynomial{B,V,E},I},
+} <: SA.MultiplicativeStructure{Polynomial{B,V,E},I}
+    basis::BT
+end
+
+function Base.:(==)(a::MStruct, b::MStruct)
+    return a.basis == b.basis
+end
+
+function MA.promote_operation(
+    ::typeof(SA.basis),
+    ::Type{MStruct{B,V,E,I,BT}},
+) where {B,V,E,I,BT}
+    return BT
+end
+
+function (m::MStruct)(a::Polynomial, b::Polynomial, ::Type{U}) where {U}
+    return m(m[a], m[b], U)
+end
+
+function (m::MStruct{B,V,E})(
+    a::E,
+    b::E,
+    ::Type{Polynomial{B,V,E}},
+) where {B,V,E}
+    return SA.map_keys(Base.Fix1(getindex, m), m(a, b, E))
+end
+
+SA.promote_with_map(::MStruct, b, m) = MStruct(b), m
+
+function SA.promote_bases_with_maps(a::MStruct, b::SA.AbstractBasis)
+    _b, _a = SA.promote_bases_with_maps(b, SA.basis(a))
+    return SA.maybe_promote(a, _a...), _b
+end
+
+function SA.promote_bases_with_maps(a::MStruct, b::MStruct)
+    _a, _b = SA.promote_bases_with_maps(SA.basis(a), SA.basis(b))
+    return SA.maybe_promote(a, _a...), SA.maybe_promote(b, _b...)
+end
+
+function SA.promote_object(
+    v::Variables,
+    m::Union{MStruct,SA.DiracMStructure},
+    map,
+)
+    return SA.promote_object(v, SA.basis(m), map)
+end
+function SA.promote_object(v::Variables, m::SA.SubBasis, map)
+    return SA.promote_object(v, parent(m), map)
+end
+SA.promote_object(::Variables, m::FullBasis, _) = m.map
diff --git a/src/mb_polynomial.jl b/src/mb_polynomial.jl
new file mode 100644
index 00000000..a07609d5
--- /dev/null
+++ b/src/mb_polynomial.jl
@@ -0,0 +1,171 @@
+variables(p::AbstractPolynomial) = variables(SA.basis(p))
+
+function Polynomial{B}(v::AbstractVariable) where {B}
+    vars = variables(v)
+    return Polynomial(Variables{B}(vars), map(_ -> 1, vars))
+end
+
+function Polynomial{B}(mono::AbstractMonomial) where {B}
+    vars = variables(mono)
+    return Polynomial(Variables{B}(vars), exponents(mono))
+end
+
+exponents(p::Polynomial) = p.exponents
+function exponents(p::Polynomial, vars)
+    (all_vars, map), _ = promote_variables_with_maps(variables(p), vars)
+    if all_vars != vars
+        throw(
+            ArgumentError("The supplied variables cannot represent this monomial"),
+        )
+    end
+    return map === nothing ? copy(exponents(p)) : map(exponents(p))
+end
+monomial(p::Polynomial) = monomial(variables(p), exponents(p))
+monomial(p::AbstractMonomial) = p
+
+function Base.copy(p::Polynomial{B}) where {B}
+    return Polynomial(Variables{B}(copy(variables(p))), copy(exponents(p)))
+end
+
+function Base.hash(p::Polynomial{B}, u::UInt) where {B}
+    return hash(
+        Tuple(
+            (v, e) for (v, e) in zip(variables(p), exponents(p)) if !iszero(e)
+        ),
+        hash(B, u),
+    )
+end
+function Base.isequal(p::Polynomial{B}, q::Polynomial{B}) where {B}
+    a, b = SA.promote_bases(p, q)
+    return isequal(exponents(a), exponents(b))
+end
+
+Base.isone(p::Polynomial) = all(iszero, p.exponents)
+isconstant(p::Polynomial) = all(iszero, p.exponents)
+# A monomial basis element is never zero
+Base.iszero(p::Polynomial) = false
+# Polynomial{Monomial,...} is its own monomial type
+monomial_type(::Type{PT}) where {PT<:Polynomial{<:AbstractMonomialIndexed}} = PT
+monomial_type(p::Polynomial{<:AbstractMonomialIndexed}) = typeof(p)
+# Ordering: derive from the variable type
+ordering(::Type{Polynomial{B,V,E}}) where {B,V,E} = ordering(V)
+# constant_monomial: return a monomial with all-zero exponents
+constant_monomial(p::Polynomial) = Polynomial(p.variables, zero(p.exponents))
+
+Base.iterate(p::Polynomial) = p, nothing
+Base.iterate(::Polynomial, ::Nothing) = nothing
+
+function Base.:(==)(p::Polynomial{B}, q::Polynomial{B}) where {B}
+    a, b = SA.promote_bases(p, q)
+    return exponents(a) == exponents(b)
+end
+
+variables(p::Polynomial) = variables(p.variables)
+nvariables(p::Polynomial) = nvariables(p.variables)
+
+monomial_type(::Type{<:SA.SparseCoefficients{K}}) where {K} = K
+monomial_type(::Type{<:SA.StarAlgebra{O,P}}) where {O<:Variables,P} = P
+# SA.Term constructor for Polynomial{...} basis elements is in mb_monomial_basis.jl
+
+
+function algebra_element(p, basis::SA.AbstractBasis)
+    return SA.AlgebraElement(p, algebra(basis))
+end
+
+function _algebra_element(p, ::Type{B}) where {B<:AbstractMonomialIndexed}
+    return algebra_element(
+        sparse_coefficients(p),
+        FullBasis{B}(variables(p)),
+    )
+end
+
+algebra_element(p::Polynomial{B}) where {B} = _algebra_element(monomial(p), B)
+polynomial(p::Polynomial) = algebra_element(p)
+
+function _show(io::IO, mime::MIME, p::Polynomial{B}) where {B}
+    if B != Monomial
+        print(io, B)
+        print(io, "(")
+    end
+    # Display monomial from variables and exponents directly
+    vars = variables(p)
+    exps = exponents(p)
+    if all(iszero, exps)
+        print(io, "1")
+    else
+        first = true
+        for (v, e) in zip(vars, exps)
+            iszero(e) && continue
+            if !first
+                print(io, "*")
+            end
+            first = false
+            show(io, mime, v)
+            if e > 1
+                print(io, "^", e)
+            end
+        end
+    end
+    if B != Monomial
+        print(io, ")")
+    end
+    return
+end
+
+function Base.show(io::IO, mime::MIME"text/latex", p::Polynomial)
+    print(io, "\$\$ ")
+    _show(io, mime, p)
+    print(io, " \$\$")
+    return
+end
+
+function Base.show(io::IO, mime::MIME"text/plain", p::Polynomial)
+    return _show(io, mime, p)
+end
+
+function Base.show(io::IO, mime::MIME"text/print", p::Polynomial)
+    return _show(io, mime, p)
+end
+
+Base.show(io::IO, p::Polynomial) = show(io, MIME"text/plain"(), p)
+Base.print(io::IO, p::Polynomial) = show(io, MIME"text/print"(), p)
+
+# zero for a Polynomial type: return an AlgebraElement with no terms
+function Base.zero(::Type{Polynomial{B,V,E}}) where {B,V,E}
+    vars = V()
+    basis = FullBasis{B}(vars)
+    sc = SA.SparseCoefficients(E[], Int[])
+    return SA.algebra_element(sc, algebra(basis))
+end
+
+function Base.zero(p::Polynomial)
+    basis = FullBasis{typeof_basis(p)}(p)
+    sc = SA.SparseCoefficients(typeof(p.exponents)[], Int[])
+    return SA.algebra_element(sc, algebra(basis))
+end
+typeof_basis(::Polynomial{B}) where {B} = B
+
+function convert_basis(basis::SA.AbstractBasis, p::AbstractTermLike)
+    return convert_basis(basis, _algebra_element(p, Monomial))
+end
+
+function convert_basis(basis::SA.AbstractBasis, p::SA.AlgebraElement)
+    return SA.AlgebraElement(SA.coeffs(p, basis), algebra(basis))
+end
+
+function Base.isapprox(a::SA.AlgebraElement, α::Number; kws...)
+    return isapprox(
+        a,
+        α * constant_algebra_element(SA.basis(a), typeof(α));
+        kws...,
+    )
+end
+
+# Type operations for AlgebraElement
+function monomial_type(::Type{<:SA.AlgebraElement{T,A}}) where {T,A}
+    return monomial_type(A)
+end
+monomial_type(::Type{<:SA.StarAlgebra{O}}) where {O} = monomial_type(O)
+function polynomial_type(::Type{A}, ::Type{T}) where {A<:SA.StarAlgebra,T}
+    return polynomial_type(MA.promote_operation(SA.basis, A), T)
+end
diff --git a/src/mb_variables.jl b/src/mb_variables.jl
new file mode 100644
index 00000000..e4980426
--- /dev/null
+++ b/src/mb_variables.jl
@@ -0,0 +1,38 @@
+Variables{B}(vars) where {B} = Variables{B,typeof(vars)}(vars)
+
+function Base.one(v::Variables)
+    return Polynomial(v, constant_monomial_exponents(v))
+end
+
+function variable_index(v::Variables, var)
+    return findfirst(isequal(var), v.variables)
+end
+
+function Base.:(==)(v::Variables{B}, w::Variables{B}) where {B}
+    return v.variables === w.variables || v.variables == w.variables
+end
+
+function Base.hash(v::Variables{B}, u::UInt) where {B}
+    return hash(v.variables, hash(B, u))
+end
+
+monomial_type(::Type{Variables{B,V}}) where {B,V} = monomial_type(V)
+# FIXME workaround for TP
+monomial_type(v::Variables) = monomial_type(prod(v.variables))
+
+constant_monomial_exponents(v::Variables) = map(_ -> 0, v.variables)
+
+function (v::Variables)(exponents)
+    return Polynomial(v, exponents)
+end
+
+variables(v::Variables) = v.variables
+nvariables(v::Variables) = length(v.variables)
+
+function _show(io::IO, mime::MIME, v::Variables{B}) where {B}
+    print(io, "$B polynomials in the variables ")
+    # We don't use the default `show` since we don't want to print the `eltype`
+    # and we want to use the `mime`
+    _show_vector(io, mime, v.variables)
+    return
+end
diff --git a/src/monomial.jl b/src/monomial.jl
index 04fde627..83f12171 100644
--- a/src/monomial.jl
+++ b/src/monomial.jl
@@ -9,9 +9,13 @@ Returns the type of the monomials of a polynomial of type `PT`.
 """
 monomial_type(::Union{M,Type{M}}) where {M<:AbstractMonomial} = M
 # AbstractMonomialLike is its own monomial type (like AbstractMonomial).
-# This breaks the term_type/monomial_type cycle for this bare abstract type.
 monomial_type(::Union{AbstractMonomialLike,Type{AbstractMonomialLike}}) = AbstractMonomialLike
-function monomial_type(::Union{PT,Type{PT}}) where {PT<:_APL}
+# For SA.Term{T,A,I}: derive monomial type from the algebra type A
+function monomial_type(::Type{SA.Term{T,A,I}}) where {T,A<:_PolynomialAlgebra,I}
+    return monomial_type(A)
+end
+# Generic fallback for other _APL types (polynomials etc.)
+function monomial_type(::Union{PT,Type{PT}}) where {PT<:AbstractPolynomialLike}
     return monomial_type(term_type(PT))
 end
 function monomial_type(
@@ -136,6 +140,10 @@ function constant_monomial(::Type{PT}) where {PT<:_APL}
     return constant_monomial(monomial_type(PT))
 end
 constant_monomial(t::AbstractTerm) = constant_monomial(monomial(t))
+constant_monomial(v::AbstractVariable) = constant_monomial(monomial(v))
+function constant_monomial(p::AbstractPolynomial)
+    return one(SA.object(parent(p)))
+end
 
 """
     map_exponents(f, m1::AbstractMonomialLike, m2::AbstractMonomialLike)
@@ -161,11 +169,6 @@ Base.one(t::AbstractMonomialLike) = constant_monomial(t)
 function MA.promote_operation(::typeof(one), MT::Type{<:AbstractMonomialLike})
     return monomial_type(MT)
 end
-# Bridge MA.operate!(one, ...) to constant_monomial for monomials
-# so that SA.Term's generic operate!(one, t) works via operate!(one, t.basis_element)
-function MA.operate!(::typeof(one), m::AbstractMonomial)
-    return MA.operate!(constant_monomial, m)
-end
 # See https://github.com/JuliaAlgebra/MultivariatePolynomials.jl/issues/82
 # By default, Base do oneunit(v::VT) = VT(one(v)).
 # This tries to convert a monomial to a variable which does not work.
diff --git a/src/operators.jl b/src/operators.jl
index 2c150274..2c100c9c 100644
--- a/src/operators.jl
+++ b/src/operators.jl
@@ -1,297 +1,51 @@
-# promoting multiplication is not a good idea
-# For example a polynomial of Float64 * a polynomial of JuMP affine expression
-# is a polynomial of JuMP affine expression but if we promote it would be a
-# polynomial of quadratic expression
-for op in [:+, :-, :(==)]
-    @eval Base.$op(p1::_APL, p2::_APL) = $op(promote(p1, p2)...)
-end
-# Promotion between `I` and `1` is `Any`.
-# Promotion between `I` and `2I` is `UniformScaling`.
-for op in [:+, :-]
-    @eval function Base.$op(p1::_APL, p2::_APL{<:LinearAlgebra.UniformScaling})
-        return $op(p1, map_coefficients(J -> J.λ, p2, nonzero = true))
-    end
-    @eval function Base.$op(p1::_APL{<:LinearAlgebra.UniformScaling}, p2::_APL)
-        return $op(map_coefficients(J -> J.λ, p1, nonzero = true), p2)
-    end
-    @eval function Base.$op(
-        p1::_APL{<:LinearAlgebra.UniformScaling},
-        p2::_APL{<:LinearAlgebra.UniformScaling},
-    )
-        return $op(map_coefficients(J -> J.λ, p1, nonzero = true), p2)
+# Variables and monomials enter the shared term/algebra arithmetic.
+for op in [:+, :-, :*]
+    @eval begin
+        Base.$op(a::AbstractMonomialLike, b::AbstractTerm) = $op(term(a), b)
+        Base.$op(a::AbstractTerm, b::AbstractMonomialLike) = $op(a, term(b))
     end
 end
-
-function Base.isapprox(t1::AbstractTermLike, t2::AbstractTermLike; kwargs...)
-    return isapprox(coefficient(t1), coefficient(t2); kwargs...) &&
-           monomial(t1) == monomial(t2)
-end
-# dot for SA.Term is defined in StarAlgebras
-
-function Base.isapprox(p1::_APL, p2::_APL; kwargs...)
-    return isapprox(promote(p1, p2)...; kwargs...)
-end
-
-# @eval $op(p::_APL, α) = $op(promote(p, α)...) would be less efficient
-for (op, fun) in [
-    (:+, :right_constant_plus),
-    (:-, :right_constant_minus),
-    (:*, :right_constant_mult),
-    (:(==), :right_constant_eq),
-]
-    @eval Base.$op(p::_APL, α) = $fun(p, α)
-end
-for (op, fun) in [
-    (:+, :left_constant_plus),
-    (:-, :left_constant_minus),
-    (:*, :left_constant_mult),
-    (:(==), :left_constant_eq),
-]
-    @eval Base.$op(α, p::_APL) = $fun(α, p)
-end
-## Fix ambiguity between above methods and methods in MA
-Base.:+(::MA.Zero, p::_APL) = MA.copy_if_mutable(p)
-Base.:+(p::_APL, ::MA.Zero) = MA.copy_if_mutable(p)
-Base.:-(::MA.Zero, p::_APL) = MA.operate(-, p)
-Base.:-(p::_APL, ::MA.Zero) = MA.copy_if_mutable(p)
-
-# Special case AbstractArrays of _APLs
-# We add these instead of relying on the broadcasting API since the above method definitions are very wide.
-# In particular, there is support for Matrices as coefficents. In order to avoid issues like #104 we therefore
-# explicitly define this (instead of implictly getting unexpected results).
 for op in [:+, :-]
-    @eval Base.$op(p::_APL, A::AbstractArray{<:_APL}) = map(f -> $op(p, f), A)
-    @eval Base.$op(A::AbstractArray{<:_APL}, p::_APL) = map(f -> $op(f, p), A)
-end
-Base.:*(p::_APL, A::AbstractArray) = map(f -> p * f, A)
-Base.:*(A::AbstractArray, p::_APL) = map(f -> f * p, A)
-Base.:/(A::AbstractArray, p::_APL) = map(f -> f / p, A)
-
-right_constant_function(::typeof(+)) = right_constant_plus
-right_constant_function(::typeof(-)) = right_constant_minus
-right_constant_function(::typeof(*)) = right_constant_mult
-function MA.operate!(op::Union{typeof(+),typeof(-),typeof(*)}, p::_APL, α)
-    return MA.operate!(right_constant_function(op), p, α)
-end
-
-MA.operate!(op::typeof(*), α, p::_APL) = MA.operate!(left_constant_mult, α, p)
-MA.operate!(op::typeof(*), p::_APL, α) = MA.operate!(right_constant_mult, p, α)
-MA.operate!(op::typeof(/), p::_APL, α) = map_coefficients!(Base.Fix2(op, α), p)
-function MA.operate_to!(output::AbstractPolynomial, op::typeof(*), α, p::_APL)
-    return MA.operate_to!(output, left_constant_mult, α, p)
-end
-function MA.operate_to!(output::AbstractPolynomial, op::typeof(*), p::_APL, α)
-    return MA.operate_to!(output, right_constant_mult, p, α)
+    @eval Base.$op(a::AbstractMonomialLike, b::AbstractMonomialLike) =
+        $op(term(a), term(b))
 end
-# Disambiguation: both args are _APL — use polynomial multiplication
-function MA.operate_to!(output::AbstractPolynomial, ::typeof(*), p::_APL, q::_APL)
-    return MA.operate_to!(output, *, polynomial(p), polynomial(q))
-end
-function MA.operate_to!(output::_APL, op::typeof(/), p::_APL, α)
-    return map_coefficients_to!(output, Base.Fix2(op, α), p)
+function Base.:*(a::AbstractMonomialLike, b::AbstractMonomialLike)
+    return monomial(a) * monomial(b)
 end
+Base.:*(a::Number, b::AbstractMonomialLike) = term(a, monomial(b))
+Base.:*(a::AbstractMonomialLike, b::Number) = term(b, monomial(a))
 
-function polynomial_merge!(
-    n1::Int,
-    n2::Int,
-    get1::F1,
-    get2::F2,
-    set::F3,
-    push::F4,
-    compare_monomials::F5,
-    combine::F6,
-    keep::F7,
-    resize::F8,
-    ::Type{term_type} = n1 > 0 ? typeof(get1(1)) :
-                        n2 > 0 ? typeof(get2(1)) : Nothing,
-) where {F1,F2,F3,F4,F5,F6,F7,F8,term_type}
-    buffer = nothing
-    i = j = k = 1
-    # Invariant:
-    # The terms p[0] -> p[k-1] are sorted and are smaller than the remaining terms.
-    # The terms p[k] -> p[i-1] are garbage.
-    # The terms p[i] -> p[end] are sorted and still need to be added.
-    # The terms q[j] -> p[end] are sorted and still need to be added.
-    # If `buffer` is not empty:
-    #   The terms in `buffer` are sorted and still need to be added.
-    #   Moreover, they are smaller than the terms p[i] -> p[end].
-    while i <= n1 && j <= n2
-        @assert buffer === nothing || isempty(buffer)
-        comp = compare_monomials(i, j)
-        if comp < 0
-            if k == i
-                t0 = get1(i)
-                if buffer === nothing
-                    buffer = DataStructures.Queue{term_type}()
-                end
-                push!(buffer, t0)
-                i += 1
-            end
-            set(k, get2(j))
-            j += 1
-            k += 1
-        elseif iszero(comp)
-            combine(i, j)
-            if keep(i)
-                if k != i
-                    @assert k < i
-                    set(k, get1(i))
-                end
-                k += 1
-            end
-            i += 1
-            j += 1
-        else
-            if k != i
-                set(k, get1(i))
-            end
-            i += 1
-            k += 1
-        end
-        while buffer !== nothing && !isempty(buffer) && j <= n2
-            @assert i == k
-            t = first(buffer)
-            comp = compare_monomials(t, j)
-            if comp <= 0
-                if comp < 0
-                    t = get2(j)
-                else
-                    t = combine(t, j)
-                end
-                j += 1
-            end
-            if comp >= 0
-                popfirst!(buffer)
-            end
-            # if `comp` is zero, we called `combine` so `t`
-            # might not be kept. If `comp` is not zero, we
-            # skip the `keep` call that might be costly.
-            if iszero(comp) && !keep(t)
-                continue
-            end
-            if k <= n1
-                push!(buffer, get1(i))
-                set(k, t)
-            else
-                push(t)
-                n1 += 1
-            end
-            i += 1
-            k += 1
-        end
-    end
-    if buffer !== nothing && !isempty(buffer)
-        @assert j == n2 + 1
-        @assert i == k
-        n = length(buffer)
-        resize(n1 + n)
-        for k in n1:-1:i
-            set(k + n, get1(k))
-        end
-        for k in i:(i+n-1)
-            set(k, popfirst!(buffer))
-        end
-        n1 += n
-    else
-        len = (k - 1) + (n2 - (j - 1)) + (n1 - (i - 1))
-        if n1 < len
-            resize(len)
-        end
-        while j <= n2
-            set(k, get2(j))
-            j += 1
-            k += 1
-        end
-        @assert j == n2 + 1
-        while i <= n1
-            set(k, get1(i))
-            i += 1
-            k += 1
-        end
-        if len < n1
-            resize(len)
-        end
-        @assert i == n1 + 1
-        @assert k == len + 1
+for op in [:+, :-]
+    @eval begin
+        Base.$op(p::AbstractMonomialLike, a::Number) =
+            $op(algebra_element(p), a)
+        Base.$op(a::Number, p::AbstractMonomialLike) =
+            $op(a, algebra_element(p))
+        Base.$op(p::AbstractTerm{T}, a::Union{T,Number}) where {T} =
+            $op(algebra_element(p), a)
+        Base.$op(a::Union{T,Number}, p::AbstractTerm{T}) where {T} =
+            $op(a, algebra_element(p))
     end
-    return
-end
-
-#MA.operate!(op::Union{typeof(+), typeof(-)}, p::AbstractPolynomial, q::AbstractPolynomial) = MA.operate_to!(p, op, p, q)
-function MA.operate!(
-    op::Union{typeof(+),typeof(-)},
-    p::AbstractPolynomial,
-    q::AbstractPolynomialLike,
-)
-    return MA.operate!(op, p, polynomial(q))
 end
+_term(α, mono) = term(α, MA.copy_if_mutable(mono))
 
-function mul_to_terms!(ts::Vector{<:AbstractTerm}, p1::_APL, p2::_APL)
-    for t1 in terms(p1)
-        for t2 in terms(p2)
-            push!(ts, t1 * t2)
-        end
-    end
-    return ts
-end
-function Base.:*(p::AbstractPolynomial, q::AbstractPolynomial)
-    return polynomial!(
-        mul_to_terms!(
-            MA.promote_operation(*, term_type(p), term_type(q))[],
-            p,
-            q,
-        ),
-    )
+function Base.isapprox(t1::AbstractTermLike, t2::AbstractTermLike; kwargs...)
+    return isapprox(coefficient(t1), coefficient(t2); kwargs...) &&
+           monomial(t1) == monomial(t2)
 end
-
-Base.isapprox(p::_APL, α; kwargs...) = isapprox(promote(p, α)...; kwargs...)
-Base.isapprox(α, p::_APL; kwargs...) = isapprox(promote(p, α)...; kwargs...)
-
 # `MA.operate(-, p)` redirects to `-p` as it assumes that `-p` can be modified
 # through the MA API without modifying `p`. We should either copy the monomial
 # here or implement a `MA.operate(-, p)` that copies it. We choose the first
 # option.
 Base.:-(m::AbstractMonomialLike) = _term(-1, MA.copy_if_mutable(m))
 Base.:-(t::AbstractTermLike) = _term(MA.operate(-, coefficient(t)), monomial(t))
-Base.:-(p::_APL) = map_coefficients(-, p)
-Base.:+(p::Union{_APL,RationalPoly}) = p
-Base.:*(p::Union{_APL,RationalPoly}) = p
-
-# Avoid adding a zero constant that might artificially increase the Newton polytope
-# Need to add polynomial conversion for type stability
-function right_constant_plus(p::_APL{S}, α::T) where {S,T}
-    return iszero(α) ? polynomial(p, MA.promote_operation(+, S, T)) :
-           p + constant_term(α, p)
-end
-function left_constant_plus(α::S, p::_APL{T}) where {S,T}
-    return iszero(α) ? polynomial(p, MA.promote_operation(+, S, T)) :
-           constant_term(α, p) + p
-end
-function MA.operate!(::typeof(right_constant_plus), p::_APL, α)
-    if !iszero(α)
-        MA.operate!(+, p, constant_term(α, p))
-    end
-    return p
-end
-function right_constant_minus(p::_APL{S}, α::T) where {S,T}
-    return iszero(α) ? polynomial(p, MA.promote_operation(-, S, T)) :
-           p - constant_term(α, p)
-end
-function left_constant_minus(α::S, p::_APL{T}) where {S,T}
-    return iszero(α) ? polynomial(-p, MA.promote_operation(-, S, T)) :
-           constant_term(α, p) - p
-end
-function MA.operate!(::typeof(right_constant_minus), p::_APL, α)
-    if !iszero(α)
-        MA.operate!(-, p, constant_term(α, p))
-    end
-    return p
-end
+Base.:+(p::Union{AbstractMonomialLike,RationalPoly}) = p
+Base.:*(p::Union{AbstractMonomialLike,RationalPoly}) = p
 
 # Coefficients and variables commute
-left_constant_mult(α, v::AbstractMonomial) = term_type(v, typeof(α))(α, v)
-left_constant_mult(α, v::AbstractVariable) = left_constant_mult(α, monomial(v)) # TODO linear term
+left_constant_mult(α, v::AbstractMonomialLike) = SA.Term(α, monomial(v))
 right_constant_mult(m::AbstractMonomialLike, α) = left_constant_mult(α, m)
+# Polynomial{Monomial,...} methods added in mb_monomial_basis.jl
 
 function left_constant_mult(α, t::SA.Term)
     return term(α * coefficient(t), monomial(t))
@@ -307,141 +61,12 @@ function right_constant_mult(p::AbstractPolynomialLike, α)
     return map_coefficients(Base.Fix2(*, α), p)
 end
 
-function MA.operate_to!(output, ::typeof(left_constant_mult), α, p::_APL)
-    return map_coefficients_to!(output, Base.Fix1(*, α), p)
-end
-function MA.operate_to!(output, ::typeof(right_constant_mult), p::_APL, α)
-    return map_coefficients_to!(output, Base.Fix2(*, α), p)
-end
-function MA.operate!(::typeof(left_constant_mult), α, p::_APL)
-    return map_coefficients!(Base.Fix1(*, α), p)
-end
-function MA.operate!(::typeof(right_constant_mult), p::_APL, α)
-    return map_coefficients!(Base.Fix2(MA.mul!!, α), p)
-end
-
-function MA.operate_to!(
-    output::AbstractMonomial,
-    ::typeof(*),
-    m1::AbstractMonomialLike,
-    m2::AbstractMonomialLike,
-)
-    return map_exponents_to!(output, +, m1, m2)
-end
-function MA.operate!(
-    ::typeof(*),
-    m1::AbstractMonomial,
-    m2::AbstractMonomialLike,
-)
-    return map_exponents!(+, m1, m2)
-end
-function Base.:*(m1::AbstractMonomialLike, m2::AbstractMonomialLike)
-    return map_exponents(+, m1, m2)
-end
-#Base.:*(m1::AbstractMonomialLike, m2::AbstractMonomialLike) = *(monomial(m1), monomial(m2))
-
-function Base.:*(m::AbstractMonomialLike, t::AbstractTermLike)
-    return left_constant_mult(coefficient(t), m * monomial(t))
-end
-function Base.:*(t::AbstractTermLike, m::AbstractMonomialLike)
-    return left_constant_mult(coefficient(t), monomial(t) * m)
-end
-function Base.:*(t1::AbstractTermLike, t2::AbstractTermLike)
-    return left_constant_mult(
-        coefficient(t1) * coefficient(t2),
-        monomial(t1) * monomial(t2),
-    )
-end
-
-function MA.operate!(::typeof(*), p::_APL, t::AbstractMonomialLike)
-    return map_exponents!(+, p, t)
-end
-Base.:*(p::_APL, t::AbstractMonomialLike) = map_exponents(+, p, t)
-Base.:*(t::AbstractTermLike, p::_APL) = polynomial!(map(te -> t * te, terms(p)))
-Base.:*(p::_APL, t::AbstractTermLike) = polynomial!(map(te -> te * t, terms(p)))
-Base.:*(p::_APL, q::_APL) = polynomial(p) * polynomial(q)
-
-# guaranteed that monomial(t1) > monomial(t2)
-function _polynomial_2terms(
-    t1::TT,
-    t2::TT,
-    ::Type{T},
-) where {TT<:SA.Term,T}
-    if iszero(t1)
-        polynomial(t2, T)
-    elseif iszero(t2)
-        polynomial(t1, T)
-    else
-        # not `polynomial!` because we `t1` and `t2` cannot be modified
-        polynomial(term_type(TT, T)[t1, t2], SortedUniqState())
-    end
-end
-
-_term(α, mono) = term(α, MA.copy_if_mutable(mono))
-
-for op in [:+, :-]
-    @eval begin
-        function Base.$op(t1::AbstractTermLike, t2::AbstractTermLike)
-            return $op(term(t1), term(t2))
-        end
-        function Base.$op(t1::SA.Term, t2::SA.Term)
-            return $op(_promote_terms(t1, t2)...)
-        end
-        function Base.$op(t1::TT, t2::TT) where {T,TT<:SA.Term{T}}
-            S = MA.promote_operation($op, T, T)
-            # t1 > t2 would compare the coefficient in case the monomials are equal
-            # and it will throw a MethodError in case the coefficients are not comparable
-            if monomial(t1) == monomial(t2)
-                return polynomial(
-                    _term($op(coefficient(t1), coefficient(t2)), monomial(t1)),
-                    S,
-                )
-            elseif monomial(t1) < monomial(t2)
-                return _polynomial_2terms(t1, $op(t2), S)
-            else
-                return _polynomial_2terms($op(t2), t1, S)
-            end
-        end
-    end
-end
-_promote_terms(t1, t2) = promote(t1, t2)
-# Promotion between `I` and `1` is `Any`.
-function _promote_terms(
-    t1::SA.Term,
-    t2::SA.Term{<:LinearAlgebra.UniformScaling},
-)
-    return _promote_terms(t1, coefficient(t2).λ * monomial(t2))
-end
-function _promote_terms(
-    t1::SA.Term{<:LinearAlgebra.UniformScaling},
-    t2::SA.Term,
-)
-    return _promote_terms(coefficient(t1).λ * monomial(t1), t2)
-end
-# Promotion between `I` and `2I` is `UniformScaling`, not `UniformScaling{Int}`.
-function _promote_terms(
-    t1::SA.Term{LinearAlgebra.UniformScaling{S}},
-    t2::SA.Term{LinearAlgebra.UniformScaling{T}},
-) where {S<:Number,T<:Number}
-    U = LinearAlgebra.UniformScaling{promote_type(S, T)}
-    return _promote_terms(
-        MA.scaling_convert(U, coefficient(t1)) * monomial(t1),
-        MA.scaling_convert(U, coefficient(t2)) * monomial(t2),
-    )
-end
-function _promote_terms(
-    t1::SA.Term{LinearAlgebra.UniformScaling{T}},
-    t2::SA.Term{LinearAlgebra.UniformScaling{T}},
-) where {T<:Number}
-    return promote(t1, t2)
-end
-
 LinearAlgebra.adjoint(v::AbstractVariable) = conj(v)
 LinearAlgebra.adjoint(m::AbstractMonomial) = conj(m)
-function LinearAlgebra.adjoint(t::SA.Term)
+function LinearAlgebra.adjoint(t::AbstractTerm)
     return _term(adjoint(coefficient(t)), adjoint(monomial(t)))
 end
-function LinearAlgebra.adjoint(p::AbstractPolynomialLike)
+function LinearAlgebra.adjoint(p::AbstractPolynomial)
     return polynomial(map(LinearAlgebra.adjoint, terms(p)))
 end
 function LinearAlgebra.adjoint(r::RationalPoly)
@@ -459,11 +84,11 @@ LinearAlgebra.ishermitian(p::AbstractPolynomialLike) = p == conj(p)
 
 LinearAlgebra.transpose(v::AbstractVariable) = v
 LinearAlgebra.transpose(m::AbstractMonomial) = m
-function LinearAlgebra.transpose(t::SA.Term)
+function LinearAlgebra.transpose(t::AbstractTerm)
     return _term(LinearAlgebra.transpose(coefficient(t)), monomial(t))
 end
-function LinearAlgebra.transpose(p::AbstractPolynomialLike)
-    return polynomial(map(LinearAlgebra.transpose, terms(p)))
+function LinearAlgebra.transpose(p::AbstractPolynomial)
+    return map_coefficients(LinearAlgebra.transpose, p; nonzero = true)
 end
 function LinearAlgebra.transpose(r::RationalPoly)
     return transpose(numerator(r)) / transpose(denominator(r))
@@ -475,8 +100,104 @@ function LinearAlgebra.dot(
 )
     return p1' * p2
 end
-LinearAlgebra.dot(x, p::AbstractPolynomialLike) = x' * p
-LinearAlgebra.dot(p::AbstractPolynomialLike, x) = p' * x
+LinearAlgebra.dot(x::Number, p::AbstractPolynomialLike) = x' * p
+LinearAlgebra.dot(p::AbstractPolynomialLike, x::Number) = p' * x
+
+_sum_product_operand(p::Union{AbstractPolynomial,AbstractTerm}) = p
+_sum_product_operand(m::AbstractMonomialLike) = term(m)
+
+for (A, B) in (
+    (AbstractPolynomialLike, AbstractPolynomialLike),
+    (Number, AbstractPolynomialLike),
+    (AbstractPolynomialLike, Number),
+)
+    @eval function MA.operate(
+        ::typeof(LinearAlgebra.dot),
+        a::AbstractArray{<:$A},
+        b::AbstractArray{<:$B},
+    )
+        # Conjugation can change the basis, so do it before basis promotion.
+        return LinearAlgebra._dot_nonrecursive(adjoint.(a), b)
+    end
+end
+
+function LinearAlgebra._dot_nonrecursive(
+    a::AbstractArray{<:AbstractPolynomialLike},
+    b::AbstractArray{<:AbstractPolynomialLike},
+)
+    return SA.sum_products(
+        map(_sum_product_operand, a),
+        map(_sum_product_operand, b),
+    )
+end
+
+function LinearAlgebra._dot_nonrecursive(
+    a::AbstractArray{<:Number},
+    b::AbstractArray{<:Union{AbstractPolynomial,AbstractTerm}},
+)
+    return SA.sum_products(a, b)
+end
+function LinearAlgebra._dot_nonrecursive(
+    a::AbstractArray{<:Union{AbstractPolynomial,AbstractTerm}},
+    b::AbstractArray{<:Number},
+)
+    return SA.sum_products(a, b)
+end
+
+# A bare monomial takes its coefficient from the numeric factor. Existing
+# terms and polynomials instead convert numeric factors to their coefficients.
+function _sum_monomial_products(a, b)
+    MA._check_same_length(a, b)
+    return sum(map(term, vec(a), vec(b)))
+end
+function LinearAlgebra._dot_nonrecursive(
+    a::AbstractArray{<:Number},
+    b::AbstractArray{<:AbstractMonomialLike},
+)
+    return _sum_monomial_products(a, b)
+end
+function LinearAlgebra._dot_nonrecursive(
+    a::AbstractArray{<:AbstractMonomialLike},
+    b::AbstractArray{<:Number},
+)
+    return _sum_monomial_products(b, a)
+end
+
+function _matrix_product_array(a, ::Type)
+    return a
+end
+function _matrix_product_array(
+    a::AbstractArray{<:AbstractMonomialLike},
+    ::Type{T},
+) where {T}
+    return map(m -> term(one(T), m), a)
+end
+
+for (A, B) in (
+    (AbstractMonomialLike, AbstractMonomialLike),
+    (AbstractMonomialLike, Union{Number,AbstractPolynomial,AbstractTerm}),
+    (Union{Number,AbstractPolynomial,AbstractTerm}, AbstractMonomialLike),
+)
+    @eval function MA.operate_to!(
+        output::VecOrMat{P},
+        ::typeof(*),
+        A::AbstractMatrix{<:$A},
+        B::AbstractVecOrMat{<:$B},
+        α::Number = true,
+    ) where {P<:AbstractPolynomial}
+        # Bare monomials acquire the numeric product's coefficient type;
+        # algebra factors keep their own coefficient types.
+        a = _matrix_product_array(
+            A,
+            eltype(B) <: Number ? coefficient_type(P) : Int,
+        )
+        b = _matrix_product_array(
+            B,
+            eltype(A) <: Number ? coefficient_type(P) : Int,
+        )
+        return MA.operate_to!(output, *, a, b, α)
+    end
+end
 
 LinearAlgebra.symmetric_type(PT::Type{<:_APL}) = PT
 LinearAlgebra.symmetric(p::_APL, ::Symbol) = p
@@ -494,50 +215,15 @@ Base.vec(vars::Tuple{Vararg{AbstractVariable}}) = [vars...]
 Base.:^(x::AbstractPolynomialLike, p::Integer) = Base.power_by_squaring(x, p)
 # ^(::SA.Term, ::Integer) is defined in StarAlgebras
 
-function MA.operate_to!(
-    output::AbstractPolynomial,
-    op::MA.AddSubMul,
-    x,
-    args::Vararg{Any,N},
-) where {N}
-    return MA.operate_to!(output, MA.add_sub_op(op), x, *(args...))
-end
-function MA.operate!(
-    op::MA.AddSubMul,
-    x::AbstractPolynomial,
-    y,
-    z,
-    args::Vararg{Any,N},
-) where {N}
-    return MA.operate!(MA.add_sub_op(op), x, *(y, z, args...))
+function MA.operate_to!(output, ::typeof(left_constant_mult), α, p::_APL)
+    return SA.map_coefficients_to!(output, Base.Fix1(*, α), p)
 end
-function MA.buffer_for(
-    ::MA.AddSubMul,
-    ::Type{<:AbstractPolynomial},
-    args::Vararg{Type,N},
-) where {N}
-    return zero(MA.promote_operation(*, args...))
+function MA.operate_to!(output, ::typeof(right_constant_mult), p::_APL, α)
+    return SA.map_coefficients_to!(output, Base.Fix2(*, α), p)
 end
-function MA.buffered_operate_to!(
-    buffer::AbstractPolynomial,
-    output::AbstractPolynomial,
-    op::MA.AddSubMul,
-    x::AbstractPolynomial,
-    y,
-    z,
-    args::Vararg{Any,N},
-) where {N}
-    product = MA.operate_to!!(buffer, *, y, z, args...)
-    return MA.operate_to!(output, MA.add_sub_op(op), x, product)
+function MA.operate!(::typeof(left_constant_mult), α, p::_APL)
+    return SA.map_coefficients!(Base.Fix1(*, α), p)
 end
-function MA.buffered_operate!(
-    buffer,
-    op::MA.AddSubMul,
-    x::AbstractPolynomial,
-    y,
-    z,
-    args::Vararg{Any,N},
-) where {N}
-    product = MA.operate_to!!(buffer, *, y, z, args...)
-    return MA.operate!(MA.add_sub_op(op), x, product)
+function MA.operate!(::typeof(right_constant_mult), p::_APL, α)
+    return SA.map_coefficients!(Base.Fix2(MA.mul!!, α), p)
 end
diff --git a/src/polynomial.jl b/src/polynomial.jl
index 009821df..feeb6781 100644
--- a/src/polynomial.jl
+++ b/src/polynomial.jl
@@ -9,17 +9,11 @@ function similar_type(::Type{PT}, ::Type{T}) where {PT<:AbstractPolynomial,T}
     return polynomial_type(PT, T)
 end
 
-function Base.similar(p::PT, ::Type{T}) where {PT<:_APL,T}
-    return convert(similar_type(PT, T), p)
-end
-
 abstract type ListState end
 abstract type UnsortedState <: ListState end
 struct MessyState <: UnsortedState end
 # No duplicates or zeros
 struct UniqState <: UnsortedState end
-sortstate(::MessyState) = SortedState()
-sortstate(::UniqState) = SortedUniqState()
 struct SortedState <: ListState end
 struct SortedUniqState <: ListState end
 
@@ -49,9 +43,6 @@ Creates a polynomial equal to `sum(f(i) * mv[i] for i in 1:length(mv))`.
 Calling `polynomial([2, 4, 1], [x, x^2*y, x*y])` should return ``4x^2y + xy + 2x``.
 """
 function polynomial end
-function polynomial(p::_APL, args::Vararg{Type,N}) where {N}
-    return polynomial!(copy(p), args...)
-end
 function polynomial(Q::AbstractMatrix, mv::AbstractVector)
     return LinearAlgebra.dot(mv, Q * mv)
 end
@@ -80,36 +71,22 @@ end
 
 polynomial(ts::AbstractVector, s::ListState = MessyState()) = sum(ts)
 
-function polynomial!(p::_APL, args::Vararg{Any,N}) where {N}
-    return convert(polynomial_type(p, args...), p)
+function polynomial!(p::_APL, args::Type...)
+    return polynomial(p, args...)
 end
 
 polynomial!(ts::AbstractVector, s::ListState = MessyState()) = sum(ts)
 
-function polynomial!(
-    ts::AbstractVector{TT},
-    s::SortedUniqState,
-) where {TT<:AbstractTerm}
-    return polynomial_type(TT)(ts)
+function Base.sum(a::AbstractArray{<:AbstractMonomialLike}; kwargs...)
+    return sum(term.(a); kwargs...)
 end
 
-function polynomial!(ts::AbstractVector{<:AbstractTerm}, s::SortedState)
-    return polynomial!(uniqterms!(ts), SortedUniqState())
-end
-function polynomial!(
-    ts::AbstractVector{<:AbstractTerm},
-    s::UnsortedState = MessyState(),
+function MA.operate(
+    ::typeof(sum),
+    a::AbstractArray{<:AbstractMonomialLike};
+    kwargs...,
 )
-    return polynomial!(sort!(ts, lt = (<)), sortstate(s))
-end
-
-_collect(v::Vector) = v
-_collect(v::AbstractVector) = collect(v)
-function polynomial(
-    ts::AbstractVector{<:AbstractTerm},
-    args::Vararg{ListState,N},
-) where {N}
-    return polynomial!(MA.mutable_copy(_collect(ts)), args...)
+    return MA.operate(sum, term.(a); kwargs...)
 end
 
 """
@@ -130,24 +107,24 @@ Returns the type that `p` would have if it was converted into a polynomial of co
 Returns the same as `polynomial_type(::PT, ::Type{T})`.
 """
 function polynomial_type end
-function polynomial_type(::Type{T}) where {T<:AbstractTerm}
-    return error("`polynomial_type` not implemented for $T")
+polynomial_type(p::_APL, args::Type...) = polynomial_type(typeof(p), args...)
+function polynomial_type(::Type{M}) where {M<:AbstractMonomialLike}
+    return polynomial_type(M, Int)
 end
-function polynomial_type(::Union{P,Type{P}}) where {P<:_APL}
-    return polynomial_type(term_type(P))
+function polynomial_type(::Type{V}, ::Type{T}) where {V<:AbstractVariable,T}
+    return polynomial_type(monomial_type(V), T)
 end
-polynomial_type(::Union{P,Type{P}}) where {P<:AbstractPolynomial} = P
-function polynomial_type(::Union{M,Type{M}}) where {M<:AbstractMonomialLike}
-    return polynomial_type(term_type(M))
+polynomial_type(::Type{P}) where {P<:AbstractPolynomial} = P
+function polynomial_type(::Type{P}, ::Type{T}) where {P<:AbstractPolynomial,T}
+    return SA.similar_type(P, T)
 end
 function polynomial_type(
-    ::Union{M,Type{M}},
-    ::Type{T},
-) where {M<:AbstractMonomialLike,T}
-    return polynomial_type(term_type(M, T))
+    ::Type{<:SA.Term{T,A,I}},
+) where {T,A<:_PolynomialAlgebra,I}
+    return polynomial_type(monomial_type(A), T)
 end
-function polynomial_type(::Union{P,Type{P}}, ::Type{T}) where {P<:_APL,T}
-    return polynomial_type(polynomial_type(P), T)
+function polynomial_type(::Type{P}, ::Type{T}) where {P<:AbstractTerm,T}
+    return polynomial_type(term_type(P, T))
 end
 function polynomial_type(
     ::Union{AbstractVector{PT},Type{<:AbstractVector{PT}}},
@@ -161,32 +138,6 @@ function polynomial_type(
     return polynomial_type(PT, T)
 end
 
-function uniqterms!(ts::AbstractVector{<:AbstractTerm})
-    i = firstindex(ts)
-    for j in Iterators.drop(eachindex(ts), 1)
-        if !iszero(ts[j])
-            if monomial(ts[i]) == monomial(ts[j])
-                ts[i] = term(
-                    MA.add!!(coefficient(ts[i]), coefficient(ts[j])),
-                    monomial(ts[i]),
-                )
-            else
-                if !iszero(ts[i])
-                    i += 1
-                end
-                ts[i] = MA.copy_if_mutable(ts[j])
-            end
-        end
-    end
-    if i < length(ts)
-        if iszero(ts[i])
-            i -= 1
-        end
-        resize!(ts, i)
-    end
-    return ts
-end
-
 """
     terms(p::AbstractPolynomialLike)
 
@@ -197,7 +148,11 @@ Returns an iterator over the nonzero terms of the polynomial `p` sorted in the d
 Calling `terms` on ``4x^2y + xy + 2x`` should return an iterator of ``[4x^2y, xy, 2x]``.
 """
 terms(t::AbstractTermLike) = OneOrZeroElementVector(iszero(t), term(t))
-terms(p::AbstractPolynomialLike) = terms(polynomial(p))
+function terms(p::AbstractPolynomial)
+    return [
+        SA.Term(parent(p), i, c) for (i, c) in SA.nonzero_pairs(SA.coeffs(p))
+    ]
+end
 
 """
     nterms(p::AbstractPolynomialLike)
@@ -227,8 +182,14 @@ Returns an iterator over the coefficients of the monomials of `X` in `p` where `
 Calling `coefficients` on ``4x^2y + xy + 2x`` should return an iterator of ``[4, 1, 2]``.
 Calling `coefficients(4x^2*y + x*y + 2x + 3, [x, 1, x*y, y])` should return an iterator of ``[2, 3, 1, 0]``.
 """
-coefficients(p::_APL{T}) where {T} = LazyMap{T}(coefficient, terms(p))
-function coefficients(p::_APL{T}, X::AbstractVector) where {T}
+function coefficients(p::AbstractPolynomial)
+    return LazyMap{coefficient_type(p)}(coefficient, terms(p))
+end
+function coefficients(p::AbstractTermLike)
+    return OneOrZeroElementVector(iszero(p), coefficient(p))
+end
+function coefficients(p::_APL, X::AbstractVector)
+    T = coefficient_type(p)
     σ, mv = sort_monomial_vector(X)
     @assert length(mv) == length(X) # no duplicate in X
     c = zeros(T, length(mv))
@@ -264,7 +225,7 @@ Calling `monomials` on ``4x^2y + xy + 2x`` should return an iterator of ``[x^2y,
 
 Calling `monomials((x, y), [1, 3], m -> degree(m, y) != 1)` should return `[x^3, x*y^2, y^3, x]` where `x^2*y` and `y` have been excluded by the filter.
 """
-monomials(p::_APL) = monomial_vector(monomial.(terms(p)))
+monomials(p::AbstractPolynomial) = monomial.(terms(p))
 monomials(t::AbstractTermLike) = OneOrZeroElementVector(iszero(t), monomial(t))
 
 function isconstant(p::_APL)
@@ -385,9 +346,9 @@ Returns the coefficient of the leading term of `p`, i.e. `coefficient(leading_te
 
 Calling `leading_coefficient` on ``4x^2y + xy + 2x`` should return ``4`` and calling it on ``0`` should return ``0``.
 """
-function leading_coefficient(p::AbstractPolynomialLike{T}) where {T}
+function leading_coefficient(p::AbstractPolynomialLike)
     if iszero(p)
-        zero(T)
+        zero(coefficient_type(p))
     else
         last(coefficients(p))
     end
@@ -413,39 +374,22 @@ function leading_monomial(p::AbstractPolynomialLike)
 end
 leading_monomial(t::AbstractTermLike) = monomial(t)
 
-#$(SIGNATURES)
-"""
-    remove_leading_term(p::AbstractPolynomialLike)
-
-Returns a polynomial with the leading term removed in the polynomial `p`.
-
-### Examples
-
-Calling `remove_leading_term` on ``4x^2y + xy + 2x`` should return ``xy + 2x``.
-"""
-function remove_leading_term(p::AbstractPolynomialLike)
-    # Iterators.drop returns an Interators.Drop which is not an AbstractVector
-    return polynomial(terms(p)[1:(end-1)], SortedUniqState())
-end
-function MA.promote_operation(
-    ::typeof(remove_leading_term),
-    ::Type{PT},
-) where {PT<:AbstractPolynomial}
-    return PT
-end
-function MA.operate(::typeof(remove_leading_term), t::AbstractTermLike)
-    return remove_leading_term(t)
+SA.remove_leading_term(m::AbstractMonomialLike) = zero(m)
+function MA.operate(::typeof(SA.remove_leading_term), m::AbstractMonomialLike)
+    return SA.remove_leading_term(m)
 end
-remove_leading_term(t::AbstractTermLike) = zero(t)
 
 function unsafe_restore_leading_term end
 function MA.operate!(
     ::typeof(unsafe_restore_leading_term),
     p::AbstractPolynomial,
-    t::AbstractTermLike,
+    t::AbstractTerm,
 )
-    # `MA.add!` will copy the coefficient of `t` so `Polynomial` redefines this
-    return MA.add!!(p, t)
+    @assert parent(p) == parent(t)
+    if !iszero(t)
+        SA.coeffs(p)[t.index] = coefficient(t)
+    end
+    return p
 end
 
 #$(SIGNATURES)
@@ -532,10 +476,10 @@ function monic(p::_APL)
     return polynomial!(_div_to_one.(terms(p), α))
 end
 monic(m::AbstractMonomialLike) = m
-monic(t::AbstractTermLike{T}) where {T} = term(one(T), monomial(t))
+monic(t::AbstractTermLike) = term(one(coefficient_type(t)), monomial(t))
 
-function _div_to_one(t::AbstractTermLike{T}, α::S) where {T,S}
-    U = Base.promote_op(/, T, S)
+function _div_to_one(t::AbstractTermLike, α::S) where {S}
+    U = Base.promote_op(/, coefficient_type(t), S)
     β = coefficient(t)
     if β == α
         term(one(U), monomial(t))
@@ -552,7 +496,7 @@ The function may return zero in which case the term is dropped.
 If the function is known to never return zero for a nonzero input, `nonzero`
 can be set to `true` to get a small speedup.
 
-See also [`map_coefficients!`](@ref) and [`map_coefficients_to!`](@ref).
+See also [`SA.map_coefficients!`](@ref) and [`SA.map_coefficients_to!`](@ref).
 
 ### Examples
 
@@ -561,16 +505,11 @@ Calling `map_coefficients(α -> mod(3α, 6), 2x*y + 3x + 1)` should return `3x +
 function map_coefficients end
 function map_coefficients(
     f::F,
-    p::AbstractPolynomialLike;
+    p::AbstractPolynomial;
     nonzero = false,
-) where {F<:Function} # 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
-    return polynomial!(
-        map_coefficients.(f, terms(p)),
-        nonzero ? SortedUniqState() : SortedState(),
-    )
+) where {F<:Function}
+    output = similar(p, Base.promote_op(f, coefficient_type(p)))
+    return SA.map_coefficients_to!(output, f, p; nonzero)
 end
 function map_coefficients(
     f::F,
@@ -580,24 +519,6 @@ function map_coefficients(
     return term(f(coefficient(t)), monomial(t))
 end
 
-"""
-    map_coefficients!(f::Function, p::AbstractPolynomialLike, nonzero = false)
-
-Mutate `p` by replacing each coefficient `α` by `f(α)`.
-The function may return zero in which case the term is dropped.
-If the function is known to never return zero for a nonzero input, `nonzero`
-can be set to `true` to get a small speedup.
-The function returns `p`, which is identically equal to the second argument.
-
-See also [`map_coefficients`](@ref) and [`map_coefficients_to!`](@ref).
-
-### Examples
-
-Let `p = 2x*y + 3x + 1`, after `map_coefficients!(α -> mod(3α, 6), p)`, `p` is
-equal to `3x + 3`.
-"""
-function map_coefficients! end
-
 function map_coefficients(
     f::F,
     p::_APL,
@@ -612,22 +533,9 @@ function map_coefficients(
     ::MA.IsMutable;
     nonzero = false,
 ) where {F<:Function}
-    return map_coefficients!(f, p; nonzero = nonzero)
+    return SA.map_coefficients!(f, p; nonzero = nonzero)
 end
 
-"""
-    map_coefficients_to!(output::AbstractPolynomialLike, f::Function, p::AbstractPolynomialLike, nonzero = false)
-
-Mutate `output` by replacing each coefficient `α` of `p` by `f(α)`.
-The function may return zero in which case the term is dropped.
-If the function is known to never returns zero for a nonzero input, `nonzero`
-can be set to `true` to get a small speedup.
-The function returns `output`, which is identically equal to the first argument.
-
-See also [`map_coefficients!`](@ref) and [`map_coefficients`](@ref).
-"""
-function map_coefficients_to! end
-
 """
     deg_num_leading_terms(p::AbstractPolynomialLike, var)
 
diff --git a/src/promote.jl b/src/promote.jl
index 6f43ea11..37747cd6 100644
--- a/src/promote.jl
+++ b/src/promote.jl
@@ -1,150 +1,61 @@
-"""
-    promote_variables(p::AbstractPolynomialLike, q::AbstractPolynomialLike)
-
-Return two polynomials over the same variables.
-"""
 function promote_variables end
 
-# MonomialLike
 Base.promote_rule(::Type{M}, ::Type{M}) where {M<:AbstractMonomialLike} = M
 function Base.promote_rule(
-    M1::Type{<:AbstractMonomialLike},
-    M2::Type{<:AbstractMonomialLike},
-)
-    return promote_type(monomial_type(M1), monomial_type(M2))
-end
-
-# SA.Term
-Base.promote_rule(::Type{T}, ::Type{T}) where {T<:SA.Term} = T
-function Base.promote_rule(
-    TS::Type{<:SA.Term{S}},
-    TT::Type{<:SA.Term{T}},
-) where {S,T}
-    U = promote_type(S, T)
-    M = promote_type(monomial_type(TS), monomial_type(TT))
-    return term_type(M, U)
-end
-# MonomialLike-Term promote
-function Base.promote_rule(
-    TS::Type{<:AbstractMonomialLike},
-    TT::Type{<:SA.Term{T}},
-) where {T}
-    U = promote_type(Int, T)
-    M = promote_type(monomial_type(TS), monomial_type(TT))
-    return term_type(M, U)
-end
-function Base.promote_rule(
-    TT::Type{<:SA.Term{T}},
-    TS::Type{<:AbstractMonomialLike},
-) where {T}
-    U = promote_type(T, Int)
-    M = promote_type(monomial_type(TT), monomial_type(TS))
-    return term_type(M, U)
-end
-# Bare AbstractMonomialLike: we can't compute monomial_type/term_type on it,
-# so return the abstract union type directly.
-function Base.promote_rule(
-    ::Type{AbstractMonomialLike},
-    ::Type{<:SA.Term{T}},
-) where {T}
-    return AbstractTermLike{promote_type(Int, T)}
-end
-function Base.promote_rule(
-    ::Type{<:SA.Term{T}},
-    ::Type{AbstractMonomialLike},
-) where {T}
-    return AbstractTermLike{promote_type(T, Int)}
-end
-function promote_rule_constant(
-    ::Type{S},
-    TT::Type{<:SA.Term{T}},
-) where {S,T}
-    return term_type(TT, promote_type(S, T))
-end
-
-# PolynomialLike
-Base.promote_rule(::Type{PT}, ::Type{PT}) where {PT<:_APL} = PT
-function Base.promote_rule(PS::Type{<:_APL}, PT::Type{<:_APL})
-    # For abstract or UnionAll types (e.g. SA.Term{T,M} where T from
-    # promote_typejoin), term_type/polynomial_type may not work.
-    # Return Union{} so Julia falls back to typejoin.
-    TS = try
-        term_type(PS)
-    catch
-        return Union{}
-    end
-    TT = try
-        term_type(PT)
-    catch
-        return Union{}
-    end
-    return polynomial_type(promote_type(TS, TT))
+    ::Type{M},
+    ::Type{N},
+) where {M<:AbstractMonomialLike,N<:AbstractMonomialLike}
+    return promote_type(monomial_type(M), monomial_type(N))
 end
-# When one side is the bare AbstractMonomialLike type, avoid calling term_type on it
-function Base.promote_rule(::Type{AbstractMonomialLike}, PT::Type{<:_APL{T}}) where {T}
-    return _apl(Int, T)
+function MA.promote_operation(
+    ::typeof(*),
+    ::Type{M},
+    ::Type{N},
+) where {M<:AbstractMonomialLike,N<:AbstractMonomialLike}
+    return promote_type(monomial_type(M), monomial_type(N))
 end
-function Base.promote_rule(PT::Type{<:_APL{T}}, ::Type{AbstractMonomialLike}) where {T}
-    return _apl(Int, T)
+function MA.promote_operation(
+    op::Union{typeof(+),typeof(-)},
+    ::Type{M},
+    ::Type{N},
+) where {M<:AbstractMonomialLike,N<:AbstractMonomialLike}
+    return polynomial_type(MA.promote_operation(*, M, N))
 end
-
-function promote_rule_constant(::Type{S}, PT::Type{<:_APL{T}}) where {S,T}
-    try
-        return polynomial_type(PT, promote_type(S, T))
-    catch
-        return Any
-    end
+function MA.promote_operation(
+    ::typeof(*),
+    ::Type{T},
+    ::Type{M},
+) where {T<:Number,M<:AbstractMonomialLike}
+    return term_type(M, T)
 end
-function Base.promote_rule(::Type{PT}, ::Type{T}) where {T,PT<:_APL}
-    return promote_rule_constant(T, PT)
+function MA.promote_operation(
+    ::typeof(*),
+    ::Type{M},
+    ::Type{T},
+) where {M<:AbstractMonomialLike,T<:Number}
+    return term_type(M, T)
 end
 
-# We don't have any information on the MultivariatePolynomials implementation,
-# so we won't be able to convert the constant to `_APL`.
-# These exact-type matches must be listed explicitly because the generic
-# promote_rule_constant(::Type{S}, ::Type{<:_APL{T}}) above would try to call
-# polynomial_type/term_type on these abstract/union types and fail.
-promote_rule_constant(::Type, ::Type{AbstractMonomialLike}) = Any
-promote_rule_constant(::Type, ::Type{AbstractPolynomialLike{T}}) where {T} = Any
-promote_rule_constant(::Type, ::Type{AbstractPolynomialLike}) = Any
-promote_rule_constant(::Type, ::Type{_APL{T}}) where {T} = Any
-promote_rule_constant(::Type, ::Type{_APL}) = Any
-
-# AbstractMonomialLike{T}
-function Base.promote_rule(
-    ::Type{AbstractMonomialLike},
-    ::Type{<:AbstractMonomialLike},
-)
-    return AbstractMonomialLike
-end
-function Base.promote_rule(
-    ::Type{<:AbstractMonomialLike},
-    ::Type{AbstractMonomialLike},
-)
-    return AbstractMonomialLike
+function MA.promote_operation(
+    op::Union{typeof(+),typeof(-),typeof(*)},
+    ::Type{T},
+    ::Type{S},
+) where {T<:AbstractTerm,S<:AbstractTerm}
+    return MA.promote_operation(op, polynomial_type(T), polynomial_type(S))
 end
-# AbstractMonomialLike vs _APL and AbstractTermLike promote rules are
-# handled above (lines 69-75 and the AbstractMonomialLike-SA.Term rules).
-
-# _APL{T}
-_apl(::Type{T}, ::Type{T}) where {T} = _APL{T}
-_apl(::Type, ::Type) = _APL
-Base.promote_rule(::Type{_APL{T}}, ::Type{<:_APL{S}}) where {S,T} = _apl(S, T)
-Base.promote_rule(::Type{<:_APL{S}}, ::Type{_APL{T}}) where {S,T} = _apl(S, T)
-
-# _APL
-Base.promote_rule(::Type{_APL}, ::Type{<:_APL}) = _APL
-Base.promote_rule(::Type{<:_APL}, ::Type{_APL}) = _APL
 
 # Rational
 function promote_rule_constant(
     ::Type{T},
     ::Type{RationalPoly{NT,DT}},
-) where {T,NT,DT}
+) where {T<:Number,NT,DT}
     return RationalPoly{promote_type(T, NT),promote_type(DT, term_type(DT))}
 end
 
-function Base.promote_rule(::Type{RT}, ::Type{T}) where {T,RT<:RationalPoly}
+function Base.promote_rule(
+    ::Type{RT},
+    ::Type{T},
+) where {T<:Number,RT<:RationalPoly}
     return promote_rule_constant(T, RT)
 end
 
@@ -173,116 +84,6 @@ function Base.promote_rule(
 ) where {PT<:_APL,RT<:RationalPoly}
     return promote_rule_rational(PT, RT)
 end
-function Base.promote_rule(
-    ::Type{RT},
-    ::Type{PT},
-) where {PT<:_APL,RT<:RationalPoly}
-    return promote_rule_rational(PT, RT)
-end
-
-# MutableArithmetics
-function MA.promote_operation(
-    op::Union{typeof(+),typeof(-)},
-    PT::Type{<:_APL{S}},
-    QT::Type{<:_APL{T}},
-) where {S,T}
-    U = MA.promote_operation(op, S, T)
-    return polynomial_type(
-        promote_type(monomial_type(PT), monomial_type(QT)),
-        U,
-    )
-end
-function MA.promote_operation(
-    ::typeof(*),
-    MT1::Type{<:AbstractMonomialLike},
-    MT2::Type{<:AbstractMonomialLike},
-)
-    return promote_type(monomial_type(MT1), monomial_type(MT2))
-end
-function MA.promote_operation(
-    ::typeof(*),
-    TT::Type{<:SA.Term{S}},
-    ST::Type{<:SA.Term{T}},
-) where {S,T}
-    UT = MA.promote_operation(*, monomial_type(TT), monomial_type(ST))
-    U = MA.promote_operation(*, S, T)
-    return promote_operation_left_constant(*, U, UT)
-end
-function MA.promote_operation(
-    ::typeof(*),
-    TT::Type{<:AbstractMonomialLike},
-    ST::Type{<:SA.Term{T}},
-) where {T}
-    UT = MA.promote_operation(*, monomial_type(TT), monomial_type(ST))
-    U = MA.promote_operation(*, Int, T)
-    return promote_operation_left_constant(*, U, UT)
-end
-function MA.promote_operation(
-    ::typeof(*),
-    TT::Type{<:SA.Term{S}},
-    ST::Type{<:AbstractMonomialLike},
-) where {S}
-    UT = MA.promote_operation(*, monomial_type(TT), monomial_type(ST))
-    U = MA.promote_operation(*, S, Int)
-    return promote_operation_left_constant(*, U, UT)
-end
-function MA.promote_operation(
-    ::typeof(*),
-    PT::Type{<:_APL{S}},
-    QT::Type{<:_APL{T}},
-) where {S,T}
-    UP = MA.promote_operation(*, monomial_type(PT), monomial_type(QT))
-    U = MA.promote_sum_mul(S, T)
-    return polynomial_type(promote_operation_left_constant(*, U, UP))
-end
-
-function promote_operation_left_constant(
-    ::typeof(*),
-    ::Type{T},
-    ::Type{M},
-) where {T,M<:AbstractMonomialLike}
-    return term_type(M, T)
-end
-
-function promote_operation_right_constant(
-    ::typeof(*),
-    ::Type{M},
-    ::Type{T},
-) where {T,M<:AbstractMonomialLike}
-    return term_type(M, T)
-end
-
-function promote_operation_left_constant(
-    ::typeof(*),
-    ::Type{T},
-    ::Type{P},
-) where {T,U,P<:_APL{U}}
-    return similar_type(P, MA.promote_operation(*, T, U))
-end
-
-function promote_operation_right_constant(
-    ::typeof(*),
-    ::Type{P},
-    ::Type{T},
-) where {T,U,P<:_APL{U}}
-    return similar_type(P, MA.promote_operation(*, U, T))
-end
-
-function MA.promote_operation(
-    ::typeof(*),
-    ::Type{T},
-    ::Type{P},
-) where {T,P<:_APL}
-    return promote_operation_left_constant(*, T, P)
-end
-
-function MA.promote_operation(
-    ::typeof(*),
-    ::Type{P},
-    ::Type{T},
-) where {T,P<:_APL}
-    return promote_operation_right_constant(*, P, T)
-end
 
 function MA.promote_operation(
     ::typeof(*),
@@ -335,9 +136,6 @@ end
 
 function (map::ExponentMap{NTuple{N,Int}})(exp::NTuple{N,Int}) where {N}
     return ntuple(map.length::Val) do i
-        # This does not have the best complexity since `findfirst`
-        # search through the whole list each time but since we're using
-        # tuples, we're probably not having a large list if indices anyway
         j = findfirst(isequal(i), map.indices)
         if isnothing(j)
             return 0
@@ -384,22 +182,16 @@ function promote_variables_with_maps(a, b)
     return (all_vars, _map(a, all_vars)), (all_vars, _map(b, all_vars))
 end
 
-function SA.promote_with_map(t::AbstractTerm, all_vars, map::ExponentMap)
-    mono, _ = SA.promote_with_map(monomial(t), all_vars, map)
-    return term(coefficient(t), mono), map
+function SA.promote_with_map(p::Polynomial{B}, vars, map::ExponentMap) where {B}
+    return Polynomial(Variables{B}(vars), map(exponents(p))), map
 end
-
-function SA.promote_with_map(p::AbstractPolynomial, all_vars, map::ExponentMap)
-    new_terms = [
-        term(
-            coefficient(t),
-            first(SA.promote_with_map(monomial(t), all_vars, map)),
-        ) for t in terms(p)
-    ]
-    return polynomial(new_terms, SortedUniqState()), map
+function SA.promote_bases_with_maps(
+    p::Polynomial{B},
+    q::Polynomial{B},
+) where {B}
+    a, b = promote_variables_with_maps(variables(p), variables(q))
+    return SA.maybe_promote(p, a...), SA.maybe_promote(q, b...)
 end
-
-function SA.promote_bases_with_maps(p::_APL, q::_APL)
-    _p, _q = promote_variables_with_maps(variables(p), variables(q))
-    return SA.maybe_promote(p, _p...), SA.maybe_promote(q, _q...)
+function promote_variables(p::AbstractMonomial, q::AbstractMonomial)
+    return SA.promote_bases(p, q)
 end
diff --git a/src/rational.jl b/src/rational.jl
index aa9235ae..a28a930a 100644
--- a/src/rational.jl
+++ b/src/rational.jl
@@ -30,7 +30,7 @@ end
 function Base.convert(::Type{RationalPoly{NT,DT}}, p::_APL) where {NT,DT}
     return convert(NT, p) / one(DT)
 end
-function Base.convert(::Type{RationalPoly{NT,DT}}, α) where {NT,DT}
+function Base.convert(::Type{RationalPoly{NT,DT}}, α::Number) where {NT,DT}
     return convert(NT, α) / one(DT)
     #convert(RationalPoly{NT, DT}, convert(NT, α))
 end
@@ -39,19 +39,24 @@ end
 RationalPoly{NT,DT}(r::RationalPoly{NT,DT}) where {NT,DT} = r
 
 Base.inv(r::RationalPoly) = r.den / r.num
-Base.inv(p::_APL{T}) where {T} = one(T) / p
-Base.:/(r::RationalPoly, p) = r.num / (r.den * p)
+Base.inv(p::_APL) = one(coefficient_type(p)) / p
+Base.:/(r::RationalPoly, p::Number) = r.num / (r.den * p)
 Base.:/(r::RationalPoly, p::_APL) = r.num / (r.den * p)
 Base.:/(r::RationalPoly, s::RationalPoly) = (r.num * s.den) / (s.num * r.den)
 function Base.:/(num::NT, den::DT) where {NT<:_APL,DT<:_APL}
     return RationalPoly{NT,DT}(num, den)
 end
-function Base.:/(num, den::_APL)
+function Base.:/(num::Number, den::_APL)
     return constant_term(num, den) / den
 end
 # Polynomial divided by coefficient is a polynomial not a rational polynomial
 # (1/den) * num would not be correct in case of noncommutative coefficients
-Base.:/(num::_APL, den) = map_coefficients(α -> α / den, num, nonzero = true)
+function Base.:/(num::AbstractMonomialLike, den::Number)
+    return map_coefficients(α -> α / den, num, nonzero = true)
+end
+function Base.:/(num::AbstractTerm{T}, den::Union{T,Number}) where {T}
+    return algebra_element(num) / den
+end
 
 function Base.:+(r::RationalPoly, s::RationalPoly)
     return (r.num * s.den + r.den * s.num) / (r.den * s.den)
@@ -61,8 +66,8 @@ function _plus(r::RationalPoly, p)
 end
 Base.:+(p::_APL, r::RationalPoly) = _plus(r, p)
 Base.:+(r::RationalPoly, p::_APL) = _plus(r, p)
-Base.:+(r::RationalPoly, α) = _plus(r, α)
-Base.:+(α, r::RationalPoly) = _plus(r, α)
+Base.:+(r::RationalPoly, α::Number) = _plus(r, α)
+Base.:+(α::Number, r::RationalPoly) = _plus(r, α)
 function Base.:-(r::RationalPoly, s::RationalPoly)
     return (r.num * s.den - r.den * s.num) / (r.den * s.den)
 end
@@ -70,15 +75,15 @@ _minus(p, s::RationalPoly) = (p * s.den - s.num) / s.den
 _minus(s::RationalPoly, p) = (s.num - p * s.den) / s.den
 Base.:-(p::_APL, r::RationalPoly) = _minus(p, r)
 Base.:-(r::RationalPoly, p::_APL) = _minus(r, p)
-Base.:-(r::RationalPoly, α) = _minus(r, α)
-Base.:-(α, r::RationalPoly) = _minus(α, r)
+Base.:-(r::RationalPoly, α::Number) = _minus(r, α)
+Base.:-(α::Number, r::RationalPoly) = _minus(α, r)
 Base.:-(r::RationalPoly) = (-r.num) / r.den
 
 Base.:*(r::RationalPoly, s::RationalPoly) = (r.num * s.num) / (r.den * s.den)
 Base.:*(p::_APL, r::RationalPoly) = (p * r.num) / r.den
 Base.:*(r::RationalPoly, p::_APL) = (r.num * p) / r.den
-Base.:*(α, r::RationalPoly) = (α * r.num) / r.den
-Base.:*(r::RationalPoly, α) = (r.num * α) / r.den
+Base.:*(α::Number, r::RationalPoly) = (α * r.num) / r.den
+Base.:*(r::RationalPoly, α::Number) = (r.num * α) / r.den
 
 Base.zero(r::RationalPoly) = zero(typeof(r))
 Base.zero(::Type{RationalPoly{NT,DT}}) where {NT,DT} = zero(NT) / one(DT)
diff --git a/src/sequences.jl b/src/sequences.jl
deleted file mode 100644
index deee5c12..00000000
--- a/src/sequences.jl
+++ /dev/null
@@ -1,72 +0,0 @@
-# Copied from TypedPolynomials
-module Sequences
-
-import MutableArithmetics
-const MA = MutableArithmetics
-
-function merge_sorted!(
-    result,
-    v1::AbstractArray,
-    v2::AbstractArray,
-    compare,
-    combine,
-    filter = x -> !iszero(x),
-)
-    i = firstindex(result)
-    i1 = firstindex(v1)
-    i2 = firstindex(v2)
-    while i1 <= lastindex(v1) && i2 <= lastindex(v2)
-        x1 = v1[i1]
-        x2 = v2[i2]
-        Δ = compare(x1, x2)
-        if Δ < 0
-            if filter(x1)
-                result[i] = x1
-                i += 1
-            end
-            i1 += 1
-        elseif Δ > 0
-            if filter(x2)
-                result[i] = x2
-                i += 1
-            end
-            i2 += 1
-        else
-            c = combine(x1, x2)
-            if filter(c)
-                result[i] = combine(x1, x2)
-                i += 1
-            end
-            i1 += 1
-            i2 += 1
-        end
-    end
-    for j in i1:lastindex(v1)
-        if filter(v1[j])
-            result[i] = v1[j]
-            i += 1
-        end
-    end
-    for j in i2:lastindex(v2)
-        if filter(v2[j])
-            result[i] = v2[j]
-            i += 1
-        end
-    end
-    resize!(result, i - 1)
-    return result
-end
-
-function merge_sorted(
-    v1::AbstractArray,
-    v2::AbstractArray,
-    compare,
-    combine = Base.:+,
-    filter = x -> !iszero(x),
-)
-    T = MA.promote_operation(combine, eltype(v1), eltype(v2))
-    result = Vector{T}(undef, length(v1) + length(v2))
-    return merge_sorted!(result, v1, v2, compare, combine, filter)
-end
-
-end
diff --git a/src/show.jl b/src/show.jl
index bae999ca..a7eef5d0 100644
--- a/src/show.jl
+++ b/src/show.jl
@@ -59,7 +59,7 @@ const unicode_subscripts = ("₀", "₁", "₂", "₃", "₄", "₅", "₆", "
 unicode_subscript(i) = join(unicode_subscripts[d+1] for d in reverse(digits(i)))
 
 # MONOMIALS
-function _show(io::IO, mime, m::AbstractMonomial)
+function _show(io::IO, mime::MIME, m::AbstractMonomial)
     if isconstant(m)
         print(io, '1')
     else
diff --git a/src/substitution.jl b/src/substitution.jl
index 8dfa0b74..0fdad4c3 100644
--- a/src/substitution.jl
+++ b/src/substitution.jl
@@ -48,6 +48,15 @@ end
 function substitute(st::_AST, v::AbstractVariable, s::Substitutions)
     return substitute(st, v, s...)
 end
+function substitute(st::_AST, v::AbstractVariable, s::Substitution...)
+    for rule in s
+        if v == first(rule)
+            return last(rule)
+        end
+    end
+    st isa Eval && throw(ArgumentError("No value provided for variable $v"))
+    return v
+end
 
 ## Monomials
 function powersubstitute(
@@ -129,9 +138,13 @@ end
 
 function substitute(st::_AST, m::AbstractMonomial, s::Substitutions)
     if isconstant(m)
-        return one(power_promote(typeof(st), variables(m), s))
+        return 1
     else
-        return powersubstitute(st, s, powers(m)...)
+        return powersubstitute(
+            st,
+            s,
+            filter(p -> !iszero(last(p)), collect(powers(m)))...,
+        )
     end
 end
 
@@ -167,14 +180,7 @@ function substitute(st::_AST, p::AbstractPolynomial, s::Substitutions)
     if iszero(p)
         _polynomial(substitute(st, zero_term(p), s))
     else
-        ts = terms(p)
-        r1 = substitute(st, ts[1], s)
-        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)
-        end
-        result
+        return sum(t -> substitute(st, t, s), terms(p))
     end
 end
 
@@ -237,3 +243,7 @@ p(y => 2) # Do not do that, this works fine with TypedPolynomials but it will no
 subs(p, s::AbstractSubstitution...) = substitute(Subs(), p, s)
 
 (p::RationalPoly)(s::AbstractSubstitution...) = p.num(s...) / p.den(s...)
+
+function (p::Union{AbstractTerm,AbstractPolynomial})(s::AbstractSubstitution...)
+    return substitute(Eval(), p, s)
+end
diff --git a/src/term.jl b/src/term.jl
index f8b07c43..cc52b4bf 100644
--- a/src/term.jl
+++ b/src/term.jl
@@ -22,9 +22,11 @@ When applied to a monomial, it create a term of type `AbstractTerm{Int}`.
 function term end
 term(coef, var::AbstractVariable) = term(coef, monomial(var))
 function term(coef, mono::AbstractMonomialLike)
-    return term_type(mono, typeof(coef))(coef, mono)
+    return SA.Term(coef, mono)
 end
-term(p::_APL) = convert(term_type(typeof(p)), p)
+# Convert any polynomial-like to a term
+term(m::AbstractMonomialLike) = term(one(Int), monomial(m))
+term(t::AbstractTerm) = t
 
 """
     term_type(p::AbstractPolynomialLike)
@@ -43,22 +45,25 @@ Returns the type of the terms of `p` but with coefficient type `T`.
 
 Returns the type of the terms of a polynomial of type `PT` but with coefficient type `T`.
 """
-term_type(::Type{<:SA.Term{T,B}}) where {T,B} = SA.Term{T,B}
-# Handle UnionAll types like SA.Term{T,M} where T (from promote_typejoin)
-term_type(::Type{<:SA.Term{<:Any,B}}) where {B} = SA.Term{<:Any,B}
-term_type(::Type{<:SA.Term{<:Any,B}}, ::Type{T}) where {T,B} = SA.Term{T,B}
-term_type(p::Type{<:_APL}, ::Type{T}) where {T} = term_type(term_type(p), T)
-term_type(::Type{M}) where {M<:AbstractMonomialLike} = term_type(M, Int)
-# Break the term_type/monomial_type cycle for bare AbstractMonomialLike:
-# term_type(AbstractMonomialLike, T) should return SA.Term{T, AbstractMonomialLike}
-# rather than recursing through the generic _APL 2-arg path.
-function term_type(::Type{AbstractMonomialLike}, ::Type{T}) where {T}
-    return SA.Term{T,AbstractMonomialLike}
+# SA.Term{T,A,I}: identity for term_type
+term_type(::Type{TT}) where {TT<:SA.Term} = TT
+# With new coefficient type: change T, keep A and I
+function term_type(::Type{<:SA.Term{<:Any,A,I}}, ::Type{T}) where {T,A,I}
+    return SA.Term{T,A,I}
 end
+term_type(::Type{M}) where {M<:AbstractMonomialLike} = term_type(M, Int)
 term_type(v::Type{<:AbstractVariable}) = term_type(monomial_type(v))
 function term_type(v::Type{<:AbstractVariable}, ::Type{T}) where {T}
     return term_type(monomial_type(v), T)
 end
+function term_type(
+    ::Type{<:SA.AlgebraElement{T,A,C}},
+) where {T,A<:_PolynomialAlgebra,C}
+    return SA.Term{T,A,SA.key_type(C)}
+end
+function term_type(::Type{P}, ::Type{T}) where {P<:AbstractPolynomial,T}
+    return term_type(term_type(P), T)
+end
 term_type(p::_APL, ::Type{T}) where {T} = term_type(typeof(p), T)
 term_type(p::_APL) = term_type(typeof(p))
 function term_type(
@@ -91,16 +96,13 @@ Calling `coefficient(2x + 4y^2 + 3, x^2)` should return ``0``.
 function coefficient end
 coefficient(t::SA.Term) = SA.coefficient(t)
 coefficient(m::AbstractMonomialLike) = 1
-function coefficient(
-    p::AbstractPolynomialLike{T},
-    m::AbstractMonomialLike,
-) where {T}
+function coefficient(p::AbstractPolynomialLike, m::AbstractMonomialLike)
     for t in terms(p)
         if monomial(t) == m
             return coefficient(t)
         end
     end
-    return zero(T)
+    return zero(coefficient_type(p))
 end
 
 """
@@ -146,13 +148,19 @@ calling `coefficient_type` on ``1.0x^2y + 2.0x`` should return `Float64` and
 calling `coefficient_type` on ``xy`` should return `Int`.
 """
 function coefficient_type(
-    ::Union{PT,Type{PT},AbstractVector{PT},Type{<:AbstractVector{PT}}},
-) where {T,PT<:AbstractPolynomialLike{T}}
-    return T
+    ::Union{AbstractVector{P},Type{<:AbstractVector{P}}},
+) where {P<:AbstractPolynomialLike}
+    return coefficient_type(P)
 end
 function coefficient_type(
-    ::Union{SA.Term{T,B},Type{SA.Term{T,B}}},
-) where {T,B}
+    ::Union{AbstractMonomialLike,Type{<:AbstractMonomialLike}},
+)
+    return Int
+end
+coefficient_type(::Union{AbstractTerm{T},Type{<:AbstractTerm{T}}}) where {T} = T
+function coefficient_type(
+    ::Union{AbstractPolynomial{T},Type{<:AbstractPolynomial{T}}},
+) where {T}
     return T
 end
 
@@ -177,7 +185,8 @@ In order to create `x^2 * y`,
 * with TypedPolynomials, use `monomial((x, y), (2, 1))`.
 """
 function monomial end
-monomial(t::SA.Term) = SA.basis_element(t)
+# basis_element returns Polynomial{Monomial,...}, monomial() converts to actual monomial
+monomial(t::AbstractTerm) = SA.basis_element(t)
 monomial(m::AbstractMonomial) = m
 
 """
@@ -194,21 +203,20 @@ constant_term(α, p) = term(α, constant_monomial(p))
 # zero should return a polynomial since it is often used to keep the result of a summation of terms.
 # For example, Base.vecdot(x::Vector{<:AbstractTerm}, y:Vector{Int}) starts with `s = zero(dot(first(x), first(y)))` and then adds terms.
 # We want `s` to start as a polynomial for this operation to be efficient.
-#Base.zero(::Type{TT}) where {T, TT<:AbstractTermLike{T}} = zero(T) * constant_monomial(TT)
-#Base.zero(t::AbstractTermLike{T}) where {T} = zero(T) * constant_monomial(t)
 """
-    zero_term(p::AbstractPolynomialLike{T}) where T
+    zero_term(p::AbstractPolynomialLike)
 
-Equivalent to `constant_term(zero(T), p)`.
+Equivalent to `constant_term(zero(coefficient_type(p)), p)`.
 
-    zero_term(α, ::Type{PT} where {T, PT<:AbstractPolynomialLike{T}}
+    zero_term(::Type{PT}) where {PT<:AbstractPolynomialLike}
 
-Equivalent to `constant_term(zero(T), PT)`.
+Equivalent to `constant_term(zero(coefficient_type(PT)), PT)`.
 """
-zero_term(::Type{PT}) where {T,PT<:AbstractPolynomialLike{T}} = constant_term(zero(T), PT)
-zero_term(p::AbstractPolynomialLike{T}) where {T} = constant_term(zero(T), p)
-zero_term(::Type{SA.Term{T,B}}) where {T,B} = constant_term(zero(T), SA.Term{T,B})
-zero_term(t::SA.Term) = constant_term(zero(coefficient(t)), t)
+function zero_term(
+    p::Union{AbstractPolynomialLike,Type{<:AbstractPolynomialLike}},
+)
+    return constant_term(zero(coefficient_type(p)), p)
+end
 
 function Base.zero(::Type{TT}) where {TT<:AbstractMonomialLike}
     return zero(polynomial_type(TT))
@@ -219,19 +227,16 @@ end
 function MA.promote_operation(::typeof(zero), PT::Type{<:AbstractMonomialLike})
     return polynomial_type(PT)
 end
-# SA.Term already defines `zero(t::Term)` and `one(t::Term)` in StarAlgebras.
-# For type-level zero, return a polynomial (for summation efficiency):
-function Base.zero(::Type{SA.Term{T,M}}) where {T,M}
-    return zero(polynomial_type(SA.Term{T,M}))
+function Base.zero(::Type{T}) where {T<:AbstractTerm}
+    return zero(polynomial_type(T))
 end
-function MA.promote_operation(::typeof(zero), PT::Type{<:SA.Term})
-    return polynomial_type(PT)
+function MA.promote_operation(::typeof(zero), ::Type{T}) where {T<:AbstractTerm}
+    return polynomial_type(T)
 end
-# `one` and `MA.promote_operation(one, ...)` for monomials is defined in monomial.jl
-# `one(t::SA.Term)` and `one(::Type{SA.Term{T,M}})` are defined in StarAlgebras
-function MA.promote_operation(::typeof(one), ::Type{SA.Term{T,M}}) where {T,M}
-    return SA.Term{T,M}
-end
-function MA.promote_operation(::typeof(one), PT::Type{<:AbstractPolynomialLike})
-    return polynomial_type(PT)
+Base.one(t::AbstractTerm) = one(algebra_element(t))
+function MA.promote_operation(
+    ::typeof(one),
+    ::Type{T},
+) where {T<:Union{AbstractTerm,AbstractPolynomial}}
+    return polynomial_type(T)
 end
diff --git a/src/types.jl b/src/types.jl
new file mode 100644
index 00000000..60c5471e
--- /dev/null
+++ b/src/types.jl
@@ -0,0 +1,33 @@
+"""A variable supplied by a polynomial frontend."""
+abstract type AbstractVariable <: MA.AbstractMutable end
+
+struct Variables{B,V}
+    variables::V
+end
+
+abstract type AbstractMonomialIndexed end
+abstract type AbstractMonomialBasis <: AbstractMonomialIndexed end
+
+"""The monomial basis."""
+struct Monomial <: AbstractMonomialBasis end
+
+"""
+    Polynomial{B,V,E}
+
+A single element of basis `B`, identified by its variables and exponents.
+Sums of basis elements are stored in `StarAlgebras.AlgebraElement`.
+"""
+struct Polynomial{B<:AbstractMonomialIndexed,V,E} <: MA.AbstractMutable
+    variables::Variables{B,V}
+    exponents::E
+end
+
+const AbstractMonomial = Polynomial{Monomial}
+const AbstractMonomialLike = Union{AbstractVariable,AbstractMonomial}
+const _PolynomialAlgebra = SA.StarAlgebra{<:Variables}
+const AbstractTerm{T} = SA.Term{T,<:_PolynomialAlgebra}
+const AbstractPolynomial{T} = SA.AlgebraElement{T,<:_PolynomialAlgebra}
+const AbstractTermLike{T} = Union{AbstractMonomialLike,AbstractTerm{T}}
+const AbstractPolynomialLike{T} =
+    Union{AbstractTermLike{T},AbstractPolynomial{T}}
+const _APL = AbstractPolynomialLike
diff --git a/test/Project.toml b/test/Project.toml
index cba9cceb..d7e66e17 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -6,6 +6,7 @@ DynamicPolynomials = "7c1d4256-1411-5781-91ec-d7bc3513ac07"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MultivariatePolynomials = "102ac46a-7ee4-5c85-9060-abc95bfdeaa3"
 MutableArithmetics = "d8a4904e-b15c-11e9-3269-09a3773c0cb0"
+SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 StarAlgebras = "0c0c59c1-dc5f-42e9-9a8b-b5dc384a6cd1"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 TypedPolynomials = "afbbf031-7a57-5f58-a1b9-b774a0fad08d"
diff --git a/test/commutative/polynomial.jl b/test/commutative/polynomial.jl
index 792b2153..77f549ed 100644
--- a/test/commutative/polynomial.jl
+++ b/test/commutative/polynomial.jl
@@ -1,6 +1,7 @@
 using Test
 using LinearAlgebra
 import MutableArithmetics as MA
+import StarAlgebras as SA
 using MultivariatePolynomials
 const MP = MultivariatePolynomials
 @testset "Polynomial" begin
@@ -242,16 +243,16 @@ const MP = MultivariatePolynomials
         p = 2x + 1
         @test map_coefficients(x -> x / 2, p, nonzero = nz) == 1.0x + 0.5
         @test map_coefficients(x -> x / 2, 3x, nonzero = nz) == 1.5x
-        @test p === map_coefficients!(x -> x + 1, p, nonzero = nz)
+        @test p === SA.map_coefficients!(x -> x + 1, p, nonzero = nz)
         @test p == 3x + 2
         q = zero(p)
-        @test q === map_coefficients_to!(q, x -> 2x, p, nonzero = nz)
+        @test q === SA.map_coefficients_to!(q, x -> 2x, p, nonzero = nz)
         @test q == 6x + 4
-        @test q === map_coefficients_to!(q, x -> 2x, 3x, nonzero = nz)
+        @test q === SA.map_coefficients_to!(q, x -> 2x, 3x, nonzero = nz)
         @test q == 6x
-        @test q === map_coefficients_to!(q, x -> 2x, x, nonzero = nz)
+        @test q === SA.map_coefficients_to!(q, x -> 2x, x, nonzero = nz)
         @test q == 2x
-        @test q === map_coefficients_to!(q, x -> 2x, x^2, nonzero = nz)
+        @test q === SA.map_coefficients_to!(q, x -> 2x, x^2, nonzero = nz)
         @test q == 2x^2
     end
 end
diff --git a/test/runtests.jl b/test/runtests.jl
index 41dce304..5708da0c 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -21,6 +21,7 @@ end
 
 if try_import(:DynamicPolynomials)
     Mod = DynamicPolynomials
+    include("sa_polynomial.jl")
     include("commutativetests.jl")
     include("noncommutativetests.jl")
 end
diff --git a/test/sa_polynomial.jl b/test/sa_polynomial.jl
new file mode 100644
index 00000000..abb35f90
--- /dev/null
+++ b/test/sa_polynomial.jl
@@ -0,0 +1,1235 @@
+module TestSAPolynomial
+
+using Test
+import DynamicPolynomials as DP
+import MultivariatePolynomials as MP
+import StarAlgebras as SA
+import LinearAlgebra
+import SparseArrays
+
+@testset "Shared polynomial representation" begin
+    DP.@polyvar x y
+    p = 2x + y
+    for a in (x, x * y)
+        @test +a === a
+        @test *(a) === a
+    end
+    @test x isa MP.AbstractPolynomialLike
+    @test x^2 isa MP.AbstractTermLike
+    @test 2x isa MP.AbstractTerm
+    @test p isa MP.AbstractPolynomial
+    @test typeof(p) === DP.Polynomial{Int}
+    @test MP.term_type(x) === DP.Term{Int}
+    @test MP.polynomial_type(x) === typeof(p)
+    @test MP.polynomial_type(p, Float64) === DP.Polynomial{Float64}
+    @test MP.coefficient_type(x^2) === Int
+    @test MP.coefficients(p) == [1, 2]
+    @test MP.coefficient.(MP.terms(p)) == [1, 2]
+    @test length(MP.monomials(p)) == 2
+    @test MP.polynomial(p) === p
+    @test MP.polynomial!(p) === p
+    @test MP.polynomial!(2x, Float64) isa DP.Polynomial{Float64}
+    @test sprint(show, typeof(p)) == "DynamicPolynomials.Polynomial{$Int}"
+    @test sprint(show, typeof(2x)) == "DynamicPolynomials.Term{$Int}"
+    @test promote_type(typeof(x), String) === Any
+    @test_throws MethodError p * "x"
+    @test_throws InexactError p * 0.5
+    @test MP.coefficients(p / 2) == [0.5, 1.0]
+    @test MP.coefficients(p // 2) == [1 // 2, 1 // 1]
+    @test div(p, 2) == MP.polynomial(x)
+    @test_throws InexactError p / 2.5
+    @test iszero(p - (x + y) - x)
+    @test iszero((x + y)^2 - (x^2 + 2x*y + y^2))
+    @test iszero(sum([x, y]) - (x + y))
+    @test iszero(sum([p, p]) - 2p)
+    @test MP.coefficients(p) == [1, 2]
+    @test iszero(sum(DP.Polynomial{Int}[]))
+    @test !isone(zero(p))
+    @test SA.coeffs(zero(p)) != SA.coeffs(p)
+    @test MP.coefficients(MP.differentiate(p, x)) == [2]
+
+    @test MP.coefficient(p, x) == 2
+    @test MP.coefficient(x, y) == 0
+    @test MP.coefficients(x, [x, y]) == [1, 0]
+    @test MP.MA.scaling(one(x)) == 1
+    @test p - y == MP.polynomial(2x)
+    @test MP.subs(p, x => 2, y => 3) == 7
+    @test iszero(MP.subs(p, x => 2) - (y + 4))
+    @test MP.monomial([x, y], [1, 0]) == MP.monomial(x)
+    @test hash(MP.monomial([x, y], [1, 0])) == hash(MP.monomial(x))
+
+    @test p(x => 2, y => 3) == 7
+    @test (2x)(x => 3) == 6
+    @test MP.coefficient_type([x, y]) === Int
+
+    DP.@polyvar u v monomial_order=MP.LexOrder
+    q = u + v
+    @test !(q isa DP.Polynomial)
+    @test MP.polynomial_type(u) === typeof(q)
+    @test iszero(sum([u, v]) - q)
+end
+
+@testset "Ordered noncommutative monomials" begin
+    DP.@ncpolyvar x y z
+    xy, yx = @inferred(x * y), @inferred(y * x)
+    @test xy != yx
+    @test MP.variables(xy) == [x, y]
+    @test MP.variables(yx) == [y, x]
+    @test MP.exponents(xy) == MP.exponents(yx) == [1, 1]
+    @test x * x == x^2
+    for (word, vars, exps) in (
+        (x * (y * x), [x, y, x], [1, 1, 1]),
+        ((x * y^2) * (y^3 * z), [x, y, z], [1, 5, 1]),
+        ((x * x * y) * (y * x), [x, y, x], [2, 2, 1]),
+        ((x * y)^2, [x, y, x, y], [1, 1, 1, 1]),
+        ((x * y)^3, [x, y, x, y, x, y], ones(Int, 6)),
+        (x^1000000 * y, [x, y], [1000000, 1]),
+    )
+        @test MP.variables(word) == vars
+        @test MP.exponents(word) == exps
+    end
+    @test x^0 * y == y
+    @test y * x^0 == y
+    @test isone((x * y)^0)
+    @test (x * y)^1 == xy
+    copied = xy^1
+    @test MP.variables(copied) !== MP.variables(xy)
+    @test MP.exponents(copied) !== MP.exponents(xy)
+    MP.variables(copied)[1] = z
+    MP.exponents(copied)[2] = 3
+    @test MP.variables(xy) == [x, y]
+    @test MP.exponents(xy) == [1, 1]
+    @test (@inferred MP.degree(x * y^4 * x^2, x)) == 3
+
+    a, b = @inferred SA.promote_bases(xy, yx)
+    @test MP.variables(a) == MP.variables(b) == [x, y, x]
+    @test MP.exponents(a) == [1, 1, 0]
+    @test MP.exponents(b) == [0, 1, 1]
+    @test a == xy && b == yx
+    @test isequal(a, xy) && hash(a) == hash(xy)
+    @test MP.exponents(xy) == MP.exponents(yx) == [1, 1]
+    @test (@inferred MP.exponents(x^2 * y * x, [x, y, x, z])) == [2, 1, 1, 0]
+    @test MP.exponents(a, MP.variables(a)) == MP.exponents(a)
+    @test MP.exponents(a, MP.variables(a)) !== MP.exponents(a)
+    @test_throws ArgumentError MP.exponents(yx, [x, y])
+
+    p = @inferred x + y
+    @test typeof(2x) === MP.term_type(x)
+    @test typeof(p) === MP.polynomial_type(x)
+    @test p == y + x
+    @test (@inferred sum([x, y])) == p
+    @test (@inferred LinearAlgebra.dot([x, y], [x, y])) == x^2 + y^2
+    @test (x + 1) * (x + 2) == x^2 + 3x + 2
+    p = MP.term(2, xy) + MP.term(3, yx)
+    @test MP.coefficient(p, xy) == 2
+    @test MP.coefficient(p, yx) == 3
+
+    DP.@polyvar u v w
+    @test (@inferred MP.exponents(u^2 * v, [u, v, w])) == [2, 1, 0]
+    @test_throws ArgumentError MP.exponents(u * v, [u])
+end
+
+@testset "Parent-aware sums" begin
+    DP.@polyvar x y z
+    left, right = @inferred SA.promote_bases(2x, 3y)
+    @test parent(left) == parent(right)
+    @test left == 2x
+    @test right == 3y
+    for a in ([x, y], [2x, 3y], [x + 1, y + 2])
+        originals = deepcopy(a)
+        expected = a[1] + a[2]
+        @test (@inferred sum(a)) == expected
+        @test (@inferred MP.MA.operate(sum, a)) == expected
+        @test sum(a; init = z + 1) == expected + z + 1
+        @test sum(a; init = 0.5) ==
+              convert(DP.Polynomial{Float64}, expected) + 0.5
+        @test a == originals
+        @test parent(sum(a)) == parent(expected)
+    end
+
+    p = x + y
+    a = [p, p]
+    original = deepcopy(p)
+    @test sum(a; init = p) == 3p
+    @test p == original
+    @test iszero(sum(DP.Polynomial{Int}[]))
+    @test iszero(sum(DP.Term{Int}[]))
+    @test iszero(sum(typeof(x)[]))
+    @test sum(DP.Polynomial{Int}[]; init = p) == p
+    @test parent(sum(DP.Polynomial{Int}[]; init = p)) === parent(p)
+
+    a = [x y; z x]
+    @test sum(a; dims = 1) == [x + z y + x]
+    @test sum(a; dims = 2) == reshape([x + y, z + x], 2, 1)
+    @test sum(a; dims = (1, 2)) == fill(2x + y + z, 1, 1)
+
+    @test MP.polynomial([2x, 3y, -2x]) == 3y
+    @test MP.differentiate(x^2 + 2x*y + y, x) == 2x + 2y
+    DP.@complex_polyvar w
+    @test adjoint((1 + im) * w + x) == (1 - im) * conj(w) + x
+end
+
+@testset "Parent-aware dot products" begin
+    DP.@polyvar x y z
+    for (a, b) in (
+        ([x, y], [y, z]),
+        ([2x, 3y], [4y, 5z]),
+        ([x + 1, y + 2], [y - 1, z + 3]),
+        ([2x, 3y], [y - 1, z + 3]),
+        ([x + 1, y + 2], [4y, 5z]),
+        ([x + 1, y + 2], [big(2) * y, big(3) * z]),
+        (
+            [big(2) * x + 1, big(3) * y + 2],
+            [big(4) * y + 1, big(5) * z + 2],
+        ),
+    )
+        originals = MP.MA.mutable_copy.((a, b))
+        expected =
+            LinearAlgebra.dot(a[1], b[1]) + LinearAlgebra.dot(a[2], b[2])
+        @test (@inferred LinearAlgebra.dot(a, b)) == expected
+        @test (@inferred MP.MA.operate(LinearAlgebra.dot, a, b)) == expected
+        @test a' * b == expected
+        @test LinearAlgebra.dot(reshape(a, 1, 2), reshape(b, 2, 1)) == expected
+        @test (a, b) == originals
+        @test_throws DimensionMismatch LinearAlgebra.dot(a, b[1:1])
+    end
+    for T in (Int, BigInt, Complex{Int}), P in (DP.Polynomial{T}, DP.Term{T})
+        @test (@inferred LinearAlgebra.dot(P[], P[])) isa DP.Polynomial{T}
+        result = LinearAlgebra.dot(P[], P[])
+        @test iszero(result)
+        @test MP.coefficient_type(result) === T
+    end
+    @test iszero(LinearAlgebra.dot(typeof(x)[], typeof(x)[]))
+
+    DP.@complex_polyvar w
+    a = DP.Polynomial{Complex{Int}}[(1 + 2im) * w + x, w + y]
+    b = DP.Polynomial{Complex{Int}}[z + w, (2 - im) * w + x]
+    originals = MP.MA.mutable_copy.((a, b))
+    @test (@inferred LinearAlgebra.dot(a, b)) ==
+          ((1 - 2im) * conj(w) + x) * b[1] + (conj(w) + y) * b[2]
+    @test a' * b == LinearAlgebra.dot(a, b)
+    @test transpose(a) * b == a[1] * b[1] + a[2] * b[2]
+    @test (a, b) == originals
+    @test LinearAlgebra.dot([(1 + im) * w, 2x], [3y, 1z]) ==
+          3 * (1 - im) * conj(w) * y + 2x * z
+end
+
+@testset "Mixed numeric and polynomial dot products" begin
+    DP.@polyvar x y
+    for polynomials in (
+        [x, y],
+        [x^2, y^2],
+        [2x, 3y],
+        [x + 1, y + 2],
+        [big(2) * x, big(3) * y],
+        [big(2) * x + 1, big(3) * y + 2],
+    ), numbers in ([2, 3], [2.0, 3.0]), left in (false, true)
+
+        a, b = left ? (numbers, polynomials) : (polynomials, numbers)
+        originals = MP.MA.mutable_copy.((a, b))
+        expected =
+            LinearAlgebra.dot(a[1], b[1]) + LinearAlgebra.dot(a[2], b[2])
+        @test (@inferred LinearAlgebra.dot(a, b)) == expected
+        @test (@inferred MP.MA.operate(LinearAlgebra.dot, a, b)) == expected
+        @test a' * b == expected
+        @test transpose(a) * b == a[1] * b[1] + a[2] * b[2]
+        @test LinearAlgebra.dot(reshape(a, 1, 2), reshape(b, 2, 1)) == expected
+        @test (a, b) == originals
+        @test_throws DimensionMismatch LinearAlgebra.dot(a, b[1:1])
+        @test_throws DimensionMismatch LinearAlgebra.dot(a[1:0], b)
+    end
+
+    for monomials in ([x, y], [x^2, y^2])
+        expected = 0.5 * monomials[1] + 1.5 * monomials[2]
+        @test (@inferred LinearAlgebra.dot([0.5, 1.5], monomials)) == expected
+        @test (@inferred LinearAlgebra.dot(monomials, [0.5, 1.5])) == expected
+    end
+    for polynomials in ([2x, 3y], [x + 1, y + 2])
+        @test_throws InexactError LinearAlgebra.dot([0.5, 1.5], polynomials)
+        @test_throws InexactError LinearAlgebra.dot(polynomials, [0.5, 1.5])
+        @test MP.coefficient_type(LinearAlgebra.dot([2.0, 3.0], polynomials)) ===
+              Int
+    end
+    numbers = BigInt[2, 3]
+    result = @inferred LinearAlgebra.dot(numbers, [x, y])
+    MP.MA.operate!(+, first(MP.coefficients(result)), big(1))
+    @test numbers == [2, 3]
+
+    for T in (Int, BigInt, Complex{Int}), P in (DP.Polynomial{T}, DP.Term{T})
+        for (a, b) in ((Float64[], P[]), (P[], Float64[]))
+            result = @inferred LinearAlgebra.dot(a, b)
+            @test result isa DP.Polynomial{T}
+            @test iszero(result)
+        end
+    end
+    for T in (Float64, BigInt), P in (typeof(x), typeof(x^2))
+        for (a, b) in ((T[], P[]), (P[], T[]))
+            result = @inferred LinearAlgebra.dot(a, b)
+            @test result isa DP.Polynomial{T}
+            @test iszero(result)
+        end
+    end
+
+    DP.@complex_polyvar w
+    numbers = [1 + 2im, 3 - im]
+    for polynomials in (
+        [w, x],
+        [w^2, x^2],
+        [(2 + im) * w, (3 - im) * x],
+        [(2 + im) * w + x, (3 - im) * x + 1],
+    ), left in (false, true)
+
+        a, b = left ? (numbers, polynomials) : (polynomials, numbers)
+        originals = MP.MA.mutable_copy.((a, b))
+        expected =
+            LinearAlgebra.dot(a[1], b[1]) + LinearAlgebra.dot(a[2], b[2])
+        @test (@inferred LinearAlgebra.dot(a, b)) == expected
+        @test a' * b == expected
+        @test transpose(a) * b == a[1] * b[1] + a[2] * b[2]
+        @test (a, b) == originals
+    end
+end
+
+@testset "Matrix products in a common basis" begin
+    DP.@polyvar x y z
+    for (A, b) in (
+        ([1 2; 3 4], [x, y]),
+        ([0.5 1.5; 2.5 3.5], [x^2, y^2]),
+        ([1 2; 3 4], [2x, 3y]),
+        ([2.0 3.0; 4.0 5.0], [x + 1, y + 2]),
+        ([1 2; 3 4], [big(2) * x + 1, big(3) * y + 2]),
+        ([x y; z x], [0.5, 1.5]),
+        ([x y; z x], [y, z]),
+        ([x^2 y^2; z^2 x*y], [y, z]),
+        ([2x 3y; 4z 1x], [x, z]),
+        ([x y; z x], [2y, 3z]),
+        ([x + 1 y + 2; z + 3 x + y], [x, y]),
+        ([x y; z x], [y + 1, z + 2]),
+        ([x + 1 y + 2; z + 3 x + y], [2, 3]),
+        ([x + 1 y + 2; z + 3 x + y], [y + 1, z + 2]),
+    )
+        for B in (b, hcat(b, reverse(b), b))
+            originals = MP.MA.mutable_copy.((A, B))
+            expected = if B isa AbstractVector
+                [A[i, 1] * B[1] + A[i, 2] * B[2] for i in 1:2]
+            else
+                [A[i, 1] * B[1, j] + A[i, 2] * B[2, j] for i in 1:2, j in 1:3]
+            end
+            result = @inferred A * B
+            @test result == expected
+            @test typeof(result) === typeof(expected)
+            @test (@inferred MP.MA.operate(*, A, B)) == expected
+            @test all(parent(r) == parent(first(result)) for r in result)
+            output = similar(result)
+            @test (@inferred LinearAlgebra.mul!(output, A, B)) === output
+            @test output == expected
+            for (α, β) in ((2, 0), (2, 3), (0, 1))
+                initial = fill(convert(eltype(result), z + 1), size(result))
+                output = iszero(β) ? similar(result) : copy(initial)
+                @test (@inferred LinearAlgebra.mul!(output, A, B, α, β)) ===
+                      output
+                @test output == expected .* α .+ initial .* β
+                @test all(c == z + 1 for c in initial)
+            end
+            @test (A, B) == originals
+            bad = B isa AbstractVector ? B[1:1] : B[1:1, :]
+            @test_throws DimensionMismatch A * bad
+        end
+    end
+
+    DP.@complex_polyvar w
+    A = DP.Polynomial{Complex{Int}}[
+        (1 + im) * w + x w + y
+        w + z (2 - im) * w + x
+    ]
+    b = [1 + im, 2 - im]
+    for matrix in (A, transpose(A), adjoint(A), view(A, :, :))
+        expected = [matrix[i, 1] * b[1] + matrix[i, 2] * b[2] for i in 1:2]
+        @test (@inferred matrix * b) == expected
+        @test (@inferred MP.MA.operate(*, matrix, b)) == expected
+    end
+    N = [1 + im 2 - im; 2 + im 1 - im]
+    for (left, right) in (
+        (transpose(A), N),
+        (adjoint(A), N),
+        (view(A, :, :), N),
+        (N, transpose(A)),
+        (N, adjoint(A)),
+        (N, view(A, :, :)),
+        (transpose(A), transpose(A)),
+        (adjoint(A), adjoint(A)),
+        (transpose(A), adjoint(A)),
+        (adjoint(A), transpose(A)),
+    )
+        expected = [
+            left[i, 1] * right[1, j] + left[i, 2] * right[2, j] for
+            i in 1:2, j in 1:2
+        ]
+        @test (@inferred left * right) == expected
+        @test (@inferred MP.MA.operate(*, left, right)) == expected
+        output = similar(expected)
+        @test (@inferred LinearAlgebra.mul!(output, left, right)) === output
+        @test output == expected
+        initial = fill(convert(eltype(expected), z + 1), size(expected))
+        output = copy(initial)
+        @test (@inferred LinearAlgebra.mul!(
+            output,
+            left,
+            right,
+            1 + im,
+            2 - im,
+        )) === output
+        @test output == expected .* (1 + im) .+ initial .* (2 - im)
+    end
+
+    for P in (typeof(x), DP.Term{Int}, DP.Polynomial{Int})
+        result = @inferred zeros(Float64, 2, 0) * P[]
+        @test length(result) == 2
+        @test all(iszero, result)
+        @test eltype(result) === MP.polynomial_type(
+            P,
+            P <: MP.AbstractMonomialLike ? Float64 : Int,
+        )
+        result = @inferred zeros(Float64, 2, 0) * Matrix{P}(undef, 0, 3)
+        @test size(result) == (2, 3)
+        @test all(iszero, result)
+        @test eltype(result) === MP.polynomial_type(
+            P,
+            P <: MP.AbstractMonomialLike ? Float64 : Int,
+        )
+        @test isempty(@inferred zeros(Int, 0, 2) * Matrix{P}(undef, 2, 0))
+    end
+    @test isempty(@inferred zeros(Int, 0, 2) * [x, y])
+    @test_throws InexactError [0.5 1.5] * [2x, 3y]
+    @test_throws InexactError [0.5 1.5] * [x + 1, y + 2]
+    @test_throws InexactError [0.5 1.5] * [2x 3y; 3y 2x]
+    @test_throws InexactError [x + 1 y + 2] * [0.5 1.5; 0.5 1.5]
+    @test MP.polynomial([1 2; 3 4], [x, y]) == x^2 + 5x*y + 4y^2
+    @test (@inferred MP.polynomial([0.5 1.5; 0.5 3.0], [x, y])) ==
+          0.5x^2 + 2x*y + 3y^2
+
+    A = [x + 1 y + 2; x + y z + 3]
+    originals = MP.MA.mutable_copy(A)
+    for B in (A, view(A, :, :))
+        @test_throws ArgumentError LinearAlgebra.mul!(A, B, [1 2; 3 4])
+        @test_throws ArgumentError LinearAlgebra.mul!(A, [1 2; 3 4], B)
+        @test A == originals
+        for (α, β) in ((0, 0), (2, 3))
+            @test_throws ArgumentError LinearAlgebra.mul!(
+                A,
+                B,
+                [1 2; 3 4],
+                α,
+                β,
+            )
+            @test_throws ArgumentError LinearAlgebra.mul!(
+                A,
+                [1 2; 3 4],
+                B,
+                α,
+                β,
+            )
+            @test A == originals
+        end
+    end
+    output = [x + 1 y + 2]
+    @test_throws DimensionMismatch LinearAlgebra.mul!(output, A, A)
+    @test output == [x + 1 y + 2]
+    @test_throws DimensionMismatch LinearAlgebra.mul!(output, [1 2], A[1:1, :])
+    @test output == [x + 1 y + 2]
+
+    @test_throws DimensionMismatch LinearAlgebra.mul!(output, A, A, 2, 3)
+    @test_throws DimensionMismatch LinearAlgebra.mul!(
+        output,
+        [1 2],
+        A[1:1, :],
+        2,
+        3,
+    )
+    @test output == [x + 1 y + 2]
+
+    output = Vector{DP.Polynomial{Int}}(undef, 1)
+    @test LinearAlgebra.mul!(output, [NaN Inf], [x + 1, y + 2], 0, 0) === output
+    @test iszero(only(output))
+    output = [z + 1]
+    @test LinearAlgebra.mul!(output, [NaN Inf], [x, y], 0, 2) === output
+    @test output == [2z + 2]
+    @test_throws InexactError LinearAlgebra.mul!(output, [2;;], [x + 1], 0.5, 0)
+    @test_throws InexactError LinearAlgebra.mul!(output, [2;;], [x + 1], 1, 0.5)
+
+    for output in ([x + 1, z + 2], [x + 1 z + 2; z + 2 x + 1])
+        initial = copy(output)
+        B =
+            output isa AbstractVector ? DP.Polynomial{Int}[] :
+            Matrix{DP.Polynomial{Int}}(undef, 0, 2)
+        @test (@inferred LinearAlgebra.mul!(
+            output,
+            zeros(Int, 2, 0),
+            B,
+            2,
+            3,
+        )) === output
+        @test output == initial .* 3
+    end
+    output = Matrix{DP.Polynomial{Int}}(undef, 0, 2)
+    @test (@inferred LinearAlgebra.mul!(output, zeros(Int, 0, 2), A, 2, 3)) ===
+          output
+
+    result =
+        @inferred zeros(Int, 2, 0) * Matrix{DP.Polynomial{BigInt}}(undef, 0, 2)
+    MP.MA.operate!(+, result[1], one(result[1]))
+    @test result == [1 0; 0 0]
+end
+
+@testset "Mixed sparse and dense matrix products" begin
+    DP.@polyvar x y z w
+    N = SparseArrays.sparse([0 2; 3 0])
+    P = SparseArrays.sparse([2, 1, 2], [1, 2, 2], [x + 1, y + 2, z + 3], 2, 2)
+    D = [x + 1 y + 2; y + 1 z + 2]
+    originals = MP.MA.mutable_copy(SparseArrays.nonzeros(P))
+    for (A, B) in (
+        (N, [x, y]),
+        (N, [2x, 3y]),
+        (N, D),
+        (D, N),
+        (transpose(N), D),
+        (D, adjoint(N)),
+        (P, [2, 3]),
+        (P, [x + 1, y + 2]),
+        (P, [1 2; 3 4]),
+        ([1 2; 3 4], P),
+        (P, D),
+        (D, P),
+        (N, P),
+        (P, N),
+        (P, P),
+    )
+        expected = if B isa AbstractVector
+            [A[i, 1] * B[1] + A[i, 2] * B[2] for i in 1:2]
+        else
+            [A[i, 1] * B[1, j] + A[i, 2] * B[2, j] for i in 1:2, j in 1:2]
+        end
+        @test (@inferred A * B) == expected
+        @test (@inferred MP.MA.operate(*, A, B)) == expected
+        output = similar(expected)
+        @test (@inferred LinearAlgebra.mul!(output, A, B)) === output
+        @test output == expected
+        if B isa AbstractMatrix
+            output = SparseArrays.spzeros(eltype(expected), size(expected)...)
+            @test (@inferred LinearAlgebra.mul!(output, A, B)) === output
+            @test output == expected
+            @test all(
+                p -> MP.variables(p) == MP.variables(sum(expected)),
+                SparseArrays.nonzeros(output),
+            )
+        end
+        for (α, β) in ((2, 0), (2, 3), (0, 1))
+            initial = fill(w + 1, size(expected))
+            output = iszero(β) ? similar(expected) : copy(initial)
+            @test (@inferred LinearAlgebra.mul!(output, A, B, α, β)) === output
+            @test output == expected .* α .+ initial .* β
+            @test initial == fill(w + 1, size(expected))
+        end
+        @test SparseArrays.nonzeros(P) == originals
+    end
+
+    DP.@complex_polyvar u
+    P = SparseArrays.sparse(
+        [2, 1, 2],
+        [1, 2, 2],
+        DP.Polynomial{Complex{Int}}[(1+im)*u+x, y+u, z+2u],
+        2,
+        2,
+    )
+    D = Complex{Int}[1 2; 3 4]
+    for wrapper in (identity, transpose, adjoint),
+        (A, B) in (
+            (wrapper(P), D),
+            (D, wrapper(P)),
+            (wrapper(P), SparseArrays.sparse(D)),
+            (SparseArrays.sparse(D), wrapper(P)),
+        )
+
+        expected =
+            [A[i, 1] * B[1, j] + A[i, 2] * B[2, j] for i in 1:2, j in 1:2]
+        @test (@inferred A * B) == expected
+        output = similar(expected)
+        @test (@inferred LinearAlgebra.mul!(output, A, B, 2, 0)) === output
+        @test output == expected .* 2
+    end
+
+    empty = SparseArrays.spzeros(DP.Polynomial{Int}, 2, 2)
+    @test (@inferred empty * [2, 3]) == [0, 0]
+    @test (@inferred [1 2; 3 4] * empty) == zeros(Int, 2, 2)
+    @test isempty(@inferred zeros(Int, 0, 2) * empty)
+    @test SparseArrays.nnz(@inferred empty * empty) == 0
+end
+
+@testset "Equality and adjoints" begin
+    DP.@polyvar x y
+    p = x + y
+    for comp in (==, isequal)
+        @test comp(zero(p), 0) && comp(0, zero(p))
+        @test comp(one(p), 1) && comp(1, one(p))
+        @test comp(MP.polynomial(x), x) && comp(x, MP.polynomial(x))
+        @test comp(MP.term(2, one(x)), 2) && comp(2, MP.term(2, one(x)))
+        @test comp(one(x), 1) && comp(1, one(x))
+        @test !comp(p, 0) && !comp(0, p)
+    end
+    @test adjoint(p) == transpose(p) == p
+    constant = MP.term(2, one(x))
+    @test hash(constant) == hash(MP.polynomial(constant)) == hash(2)
+    @test haskey(Dict{Any,Bool}(constant => true), 2)
+    t = (1 + 2im) * x
+    @test adjoint(t) == (1 - 2im) * x
+    @test transpose(t) == t
+    c = t + (3 - im) * y
+    @test adjoint(c) == (1 - 2im) * x + (3 + im) * y
+    @test transpose(c) == c
+    @test LinearAlgebra.dot(p, p) == p * p
+    @test LinearAlgebra.dot(t, t) == 5 * x^2
+    @test LinearAlgebra.dot([p, p], [p, p]) == 2 * p * p
+    @test LinearAlgebra.dot([t, t], [t, t]) == 10 * x^2
+end
+
+function same_basis_promotion_allocations(output, x)
+    SA._promote_operand(+, output, x)
+    return @allocated SA._promote_operand(+, output, x)
+end
+
+@testset "Prepared destinations promote operands" begin
+    DP.@polyvar x y z
+    for T in (Int, BigInt)
+        a = convert(DP.Polynomial{T}, 2x + 2)
+        b = convert(DP.Polynomial{T}, y + 1)
+        prototype = convert(DP.Polynomial{T}, x + y + z)
+        alg = parent(prototype)
+        originals = MP.MA.mutable_copy.((a, b))
+        same_basis = first(SA.promote_bases(a, alg))
+        @test (@inferred SA._promote_operand(+, prototype, same_basis)) ===
+              same_basis
+        @test same_basis_promotion_allocations(prototype, same_basis) == 0
+        for op in (+, -, *)
+            output = zero(prototype)
+            @test (@inferred MP.MA.operate_to!(output, op, a, b)) === output
+            @test output == op(a, b)
+            @test parent(output) === alg
+        end
+        output = zero(prototype)
+        @test (@inferred MP.MA.operate_to!(output, -, a)) === output
+        @test output == -a
+        @test (@inferred MP.MA.operate_to!(output, *, 2, a)) === output
+        @test output == 2a
+        for op in (*, /, //, div)
+            @test (@inferred MP.MA.operate_to!(output, op, a, 2)) === output
+            @test output == op(a, 2)
+        end
+        for op in (MP.MA.add_mul, MP.MA.sub_mul),
+            (left, right) in
+            ((a, b), (2x, b), (a, 3y), (2x, 3y), (2, b), (b, 2))
+
+            output = zero(prototype)
+            @test (@inferred MP.MA.operate_to!(output, op, a, left, right)) ===
+                  output
+            @test output == op(a, left, right)
+            @test parent(output) === alg
+        end
+        @test (a, b) == originals
+        @test parent(a) != alg && parent(b) != alg
+    end
+end
+
+@testset "Fused polynomial products" begin
+    DP.@polyvar x y
+    for T in (Int, BigInt), op in (MP.MA.add_mul, MP.MA.sub_mul)
+        a, b = SA.promote_bases(
+            convert(DP.Polynomial{T}, x + 1),
+            convert(DP.Polynomial{T}, y + 2),
+        )
+        for (left, right) in ((a, b), (last(MP.terms(a)), last(MP.terms(b))))
+            f = convert(DP.Polynomial{T}, x + y)
+            originals = MP.MA.copy_if_mutable.((f, left, right))
+            expected = op(f, left, right)
+            out = zero(f)
+            @test (@inferred MP.MA.operate_to!!(out, op, f, left, right)) === out
+            @test out == expected
+            @test (f, left, right) == originals
+            @test (@inferred MP.MA.operate!!(op, f, left, right)) === f
+            @test f == expected
+            @test (left, right) == originals[2:3]
+        end
+        f = MP.MA.mutable_copy(a)
+        expected = op(f, f, b)
+        @test (@inferred MP.MA.operate!!(op, f, f, b)) === f
+        @test f == expected
+    end
+end
+
+@testset "Fused scalar products" begin
+    DP.@polyvar x y
+    for T in (Int, BigInt), op in (MP.MA.add_mul, MP.MA.sub_mul)
+        g = convert(DP.Polynomial{T}, x + 2y)
+        for p in (g, last(MP.terms(g))), (a, b) in ((2.0, p), (p, 2.0))
+            f = one(g)
+            expected = op(f, a, b)
+            original = MP.MA.mutable_copy(g)
+            out = zero(f)
+            @test (@inferred MP.MA.operate_to!!(out, op, f, a, b)) === out
+            @test out == expected
+            @test (@inferred MP.MA.operate!!(op, f, a, b)) === f
+            @test f == expected
+            @test g == original
+            # Mutating an inserted exponent vector must not modify the source.
+            fill!(last(keys(SA.coeffs(f))), 0)
+            @test g == original
+        end
+    end
+end
+
+@testset "Matching coefficient type" begin
+    DP.@polyvar x
+    a = [1 2; 3 4]
+    t = MP.term(a, x)
+    c = one(x)
+    for p in (t + a, a + t)
+        @test MP.coefficient(p, x) == a
+        @test MP.coefficient(p, c) == a
+    end
+    @test MP.coefficient(t - a, c) == -a
+    @test MP.coefficient(a - t, x) == -a
+    @test MP.term(a, c) == a
+    @test MP.polynomial(MP.term(a, c)) == a
+end
+
+@testset "Coefficient types in division and GCD" begin
+    DP.@polyvar x y
+    algo = MP.SubresultantAlgorithm()
+    immutable = MP.MA.IsNotMutable()
+    for (p, T, field) in (
+        (x, Int, Rational{Int}),
+        (x^2, Int, Rational{Int}),
+        (2x, Int, Rational{Int}),
+        (2x + 1, Int, Rational{Int}),
+        ((2 // 3) * x, Rational{Int}, Rational{Int}),
+        ((2 // 3) * x + 1, Rational{Int}, Rational{Int}),
+        (2.0x, Float64, Float64),
+        (2.0x + 1, Float64, Float64),
+    )
+        P = typeof(p)
+        for op in (div, rem)
+            @test MP.MA.promote_operation(op, P, P) === DP.Polynomial{field}
+        end
+        for op in (MP.pseudo_rem, MP.rem_or_pseudo_rem)
+            @test MP.MA.promote_operation(op, P, P, typeof(algo)) ===
+                  DP.Polynomial{T}
+        end
+    end
+    for p in (x, x^2)
+        @test (@inferred MP.content(p, algo, immutable)) === 1
+        @test MP.primitive_part(p, algo, immutable) == p
+    end
+    @test MP.content(6x, algo, immutable) == 6
+    @test MP.content(6x + 9, algo, immutable) == 3
+    @test MP.content(zero(x + y), algo, immutable) == 0
+    @test MP.primitive_part(6x + 9, algo, immutable) == 2x + 3
+    for p in (2.5x, 2.5x + 1)
+        @test MP.content(p, algo, immutable) === 1.0
+        @test MP.primitive_part(p, algo, immutable) === p
+    end
+    complex_term = (2.0 + 3.0im) * x
+    @test MP.primitive_part(complex_term, algo, immutable) === complex_term
+    @test gcd(x^2, x^3) == x^2
+    @test gcd(6x, 9x) == 3x
+    @test isempty(Test.detect_unbound_args(MP))
+end
+
+@testset "Monomial divisibility" begin
+    DP.@polyvar x y z
+    for (a, b, expected) in (
+        (x, x, true),
+        (x, y, false),
+        (x, x^2, true),
+        (x^2, x, false),
+        (x^2, x^2, true),
+        (x^2, y^2, false),
+        (x * y, x^2 * y^3, true),
+        (x * y^3, x^2 * y, false),
+        (y, x * y * z, true),
+        (x * z, y^2, false),
+        (x * y, x^2, false),
+        (y * z, y^2, false),
+        (one(x), x^2, true),
+        (x^2, one(x), false),
+        (one(x), one(y), true),
+        (MP.monomial([x, y, z], [0, 2, 0]), y^3, true),
+        (y, MP.monomial([x, y, z], [0, 2, 0]), true),
+        (MP.term(3, x * y), MP.term(2, x^2 * y), true),
+        (MP.term(2, x^2 * y), MP.term(3, x * y), false),
+    )
+        originals = deepcopy((a, b))
+        @test MP.divides(a, b) === expected
+        @test (a, b) == originals
+    end
+    @test (@inferred MP.divides(x, x^2))
+    DP.@ncpolyvar u v
+    @test MP.divides(u, u)
+    @test !MP.divides(u, v)
+    @test_throws ErrorException MP.divides(u, u^2)
+    @test_throws ErrorException MP.divides(u^2, v^2)
+end
+
+@testset "Polynomial conversion" begin
+    DP.@polyvar x
+    R = DP.Polynomial{Rational{Int}}
+    for p in (x, x^2, 2x, 2x + 1, 3, 0)
+        converted = convert(R, p)
+        @test converted isa R
+        @test converted == p
+        @test convert(R, converted) === converted
+    end
+    p = 2x + 1
+    @test parent(convert(R, p)) === parent(p)
+    @test_throws InexactError convert(DP.Polynomial{Int}, (1 // 2) * x + 1)
+    @test_throws InexactError convert(DP.Polynomial{Int}, (1 // 2) * x)
+    @test_throws InexactError convert(DP.Polynomial{Int}, 1 // 2)
+end
+
+@testset "Polynomial conversion in division and GCD" begin
+    DP.@polyvar x y
+    for p in (x, x^2, 2x, 2x + 1)
+        original = deepcopy(p)
+        q, r = divrem(p, typeof(p)[])
+        @test isempty(q)
+        @test r == p
+        @test r isa DP.Polynomial{Rational{Int}}
+        MP.MA.operate!(zero, r)
+        @test p == original
+    end
+    q, r = divrem(zero(2x), 2x)
+    @test q isa DP.Polynomial{Rational{Int}}
+    @test iszero(q) && iszero(r)
+
+    algo = MP.GeneralizedEuclideanAlgorithm()
+    g = MP.primitive_univariate_gcd!(2x, zero(2x), algo)
+    @test g isa DP.Polynomial{Int}
+    @test g == 2x
+    @test MP.inflate(3, one(x), one(x)) == 3
+    nested = MP.term(2x, MP.monomial(y))
+    c = MP.content(nested, algo, MP.MA.IsNotMutable())
+    @test c isa DP.Term{Int}
+    @test c == 2x
+end
+
+@testset "Leading-term removal and restoration" begin
+    DP.@polyvar x
+    p = big(2) * x + 3
+    owned = MP.MA.copy_if_mutable(p)
+    @test parent(owned) === parent(p)
+    @test owned !== p
+    @test last(keys(SA.coeffs(owned))) !== last(keys(SA.coeffs(p)))
+    MP.MA.operate!(+, MP.leading_coefficient(owned), 1)
+    last(keys(SA.coeffs(owned)))[1] = 2
+    @test p == big(2) * x + 3
+    @test MP.MA.operate!!(SA.remove_leading_term, owned) === owned
+    @test owned == 3
+    for p in (x, x^2, 2x)
+        remainder = MP.MA.operate(SA.remove_leading_term, p)
+        @test typeof(remainder) === typeof(zero(p))
+        @test iszero(remainder)
+    end
+    for (p, expected) in (
+        (zero(x + 1), 0),
+        (MP.polynomial(2x), 0),
+        (x^2 + x + 1, x + 1),
+        (big(2) * x + 3, 3),
+    )
+        original = deepcopy(p)
+        leading = MP.leading_term(p)
+        remainder = SA.remove_leading_term(p)
+        @test remainder == expected
+        @test typeof(remainder) ===
+              MP.MA.promote_operation(SA.remove_leading_term, typeof(p)) ===
+              typeof(p)
+        @test parent(remainder) === parent(p)
+        @test p == original
+        @test MP.MA.operate!(SA.remove_leading_term, p) === p
+        @test p == expected
+        @test all(c -> !iszero(c), values(SA.coeffs(p)))
+        @test MP.MA.operate!(MP.unsafe_restore_leading_term, p, leading) === p
+        @test p == original
+        if !iszero(leading)
+            @test MP.leading_coefficient(p) === MP.coefficient(leading)
+        end
+    end
+
+    basis = SA.SubBasis(MP.FullBasis{MP.Monomial}([x]), [[0], [1], [2]])
+    for c in ([1, 0, 3], SA.SparseArrays.sparsevec([1, 3], [1, 3], 3))
+        p = SA.AlgebraElement(c, MP.algebra(basis))
+        leading = MP.leading_term(p)
+        @test SA.remove_leading_term(p) == 1
+        @test MP.MA.operate!(SA.remove_leading_term, p) === p
+        @test p == 1
+        @test MP.MA.operate!(MP.unsafe_restore_leading_term, p, leading) === p
+        @test p == 3x^2 + 1
+    end
+
+    coefficient = [1 2; 3 4]
+    p = MP.polynomial(MP.term(coefficient, x))
+    leading = MP.leading_term(p)
+    MP.MA.operate!(SA.remove_leading_term, p)
+    @test iszero(p)
+    @test MP.coefficient(leading) == [1 2; 3 4]
+    MP.MA.operate!(MP.unsafe_restore_leading_term, p, leading)
+    @test MP.leading_coefficient(p) === coefficient
+end
+
+@testset "Shared coefficient mapping" begin
+    DP.@polyvar x
+    for T in (Int, BigInt)
+        p = convert(DP.Polynomial{T}, 2x^2 + 3x + 2)
+        original = deepcopy(p)
+        mapped = @inferred MP.map_coefficients(c -> c // 2, p)
+        @test mapped == x^2 + (3 // 2) * x + 1
+        @test MP.coefficient_type(mapped) === Rational{T}
+        @test parent(mapped) === parent(p)
+        @test SA.coeffs(mapped) !== SA.coeffs(p)
+        @test MP.map_coefficients(c -> c - 2, p) == x
+        @test MP.map_coefficients(c -> 2c, p; nonzero = true) == 2p
+        @test p == original
+
+        z = zero(p)
+        mapped = @inferred MP.map_coefficients(c -> c // 2, z)
+        @test iszero(mapped)
+        @test MP.coefficient_type(mapped) === Rational{T}
+        @test parent(mapped) === parent(z)
+
+        output = zero(p)
+        @test SA.map_coefficients_to!(output, c -> c - 2, p) === output
+        @test output == x
+        @test p == original
+        @test MP.MA.operate!(MP.right_constant_mult, p, T(2)) === p
+        @test p == 2original
+        @test MP.right_constant_div_multiple(p, T(2), MP.MA.IsMutable()) === p
+        @test p == original
+        @test SA.map_coefficients!(zero, p) === p
+        @test iszero(p)
+
+        f = convert(DP.Polynomial{T}, x^3 + 2x + 1)
+        g = convert(DP.Polynomial{T}, 2x^2 + 1)
+        g_original = deepcopy(g)
+        @test MP.MA.operate!(
+            MP.pseudo_rem,
+            f,
+            g,
+            MP.GeneralizedEuclideanAlgorithm(),
+        ) === f
+        @test f == 3x + 2
+        @test g == g_original
+    end
+
+    basis = SA.SubBasis(MP.FullBasis{MP.Monomial}([x]), [[0], [1], [2]])
+    for c in ([2, 0, 4], SA.SparseArrays.sparsevec([1, 3], [2, 4], 3))
+        p = SA.AlgebraElement(c, MP.algebra(basis))
+        mapped = @inferred MP.map_coefficients(c -> c / 2 + 1, p)
+        @test parent(mapped) === parent(p)
+        @test SA.coeffs(mapped) == [2.0, 0.0, 3.0]
+        @test SA.coeffs(p) == [2, 0, 4]
+    end
+
+    p = x + 1
+    c = SA.SparseCoefficients(([0], [1]), (2, 3), SA.coeffs(p).isless)
+    p = SA.AlgebraElement(c, parent(p))
+    mapped = @inferred MP.map_coefficients(c -> c / 2, p)
+    @test mapped == 1 + 1.5x
+    @test mapped isa DP.Polynomial{Float64}
+    @test parent(mapped) === parent(p)
+    @test values(c) == (2, 3)
+
+    a = [1 2; 3 4]
+    p = MP.polynomial(MP.term(a, x))
+    q = @inferred transpose(p)
+    @test MP.coefficient(q, x) == transpose(a)
+    @test parent(q) === parent(p)
+    @test MP.coefficient_type(q) === typeof(transpose(a))
+    @test MP.coefficient(p, x) == a
+end
+
+@testset "Coefficient mapping deprecations" begin
+    DP.@polyvar x
+    p = x + 1
+    @test_deprecated r"SA\.map_coefficients!" MP.map_coefficients!(
+        identity,
+        p;
+        nonzero = true,
+    )
+    @test_deprecated r"SA\.map_coefficients_to!" MP.map_coefficients_to!(
+        p,
+        identity,
+        p;
+        nonzero = true,
+    )
+end
+
+@testset "Remainder operation routing" begin
+    DP.@polyvar x y
+    algo = MP.GeneralizedEuclideanAlgorithm()
+    for T in (Int, BigInt, Rational{Int}, Float64)
+        f = convert(DP.Polynomial{T}, x^3 + 2x + 1)
+        g = convert(DP.Polynomial{T}, 2x^2 + 1)
+        originals = deepcopy((f, g))
+        for op in (MP.pseudo_rem, MP.rem_or_pseudo_rem)
+            expected =
+                op === MP.pseudo_rem || T <: Integer ? 3x + 2 : (3 // 2) * x + 1
+            r = @inferred op(f, g, algo)
+            @test r == expected
+            @test typeof(r) === MP.MA.promote_operation(
+                op,
+                typeof(f),
+                typeof(g),
+                typeof(algo),
+            )
+            @test (f, g) == originals
+
+            owned = MP.MA.mutable_copy(f)
+            @test (@inferred MP.MA.operate!!(op, owned, g, algo)) === owned
+            @test owned == expected
+            @test g == originals[2]
+        end
+    end
+
+    for op in (MP.pseudo_rem, MP.rem_or_pseudo_rem)
+        for (f, g, expected) in (
+            (x, 2x + 1, -1),
+            (x^2, 2x, 0),
+            (2x, x^2 + 1, 2x),
+            (x^3 + 2x + 1, convert(DP.Polynomial{BigInt}, 2x^2 + 1), 3x + 2),
+        )
+            originals = deepcopy((f, g))
+            r = @inferred MP.MA.operate!!(op, f, g, algo)
+            @test r == expected
+            @test (f, g) == originals
+        end
+
+        f = x^3 + 2x + 1
+        for g in (f, SA.AlgebraElement(SA.coeffs(f), parent(f)))
+            original = deepcopy(f)
+            @test iszero(MP.MA.operate!!(op, f, g, algo))
+            @test f == original
+        end
+
+        f = x^3 + 2x + 1
+        g = 2x^2 + 1
+        buffer = MP.MA.buffer_for(op, typeof(f), typeof(g), typeof(algo))
+        @test MP.MA.buffered_operate!!(buffer, op, f, g, algo) === f
+        @test f == 3x + 2
+        @test g == 2x^2 + 1
+
+        f = x^3 + 2x + 1
+        originals = deepcopy((f, g))
+        # A missing operate_to! method must not silently allocate.
+        @test_throws ErrorException MP.MA.operate_to!!(zero(f), op, f, g, algo)
+        @test (f, g) == originals
+
+        g = 2x^2 + 1 + zero(x + y)
+        originals = deepcopy((f, g))
+        @test_throws ArgumentError MP.MA.operate!!(op, f, g, algo)
+        @test (f, g) == originals
+        @test op(f, g, algo) == 3x + 2
+        @test (f, g) == originals
+    end
+
+    f = x^3 + 2x + 1
+    g = 2x^2 + 1
+    originals = deepcopy((f, g))
+    r = @inferred MP.MA.operate!!(rem, f, g, algo)
+    @test r == (3 // 2) * x + 1
+    @test (f, g) == originals
+    @test MP.MA.operate!!(rem, r, g, algo) === r
+end
+
+@testset "In-place quotient term accumulation" begin
+    DP.@polyvar x y
+    for T in (Int, BigInt)
+        q = convert(DP.Polynomial{T}, x^2 + 1)
+        t = MP.term(T(2), x)
+        original = deepcopy(t)
+        @test (@inferred MP.MA.add!!(q, t)) === q
+        @test q == x^2 + 2x + 1
+        @test MP.MA.add!!(q, -t) === q
+        @test q == x^2 + 1
+        @test t == original
+
+        q = zero(q)
+        @test MP.MA.add!!(q, t) === q
+        stored = only(keys(SA.coeffs(q)))
+        @test stored !== t.index
+        @test stored == t.index
+        stored[1] += 1
+        @test t == original
+        if T === BigInt
+            MP.MA.operate!(+, MP.leading_coefficient(q), 1)
+            @test t == original
+        end
+    end
+
+    q = x + 1
+    t = 0.5x
+    original = deepcopy((q, t))
+    result = @inferred MP.MA.add!!(q, t)
+    @test result == 1.5x + 1
+    @test result isa DP.Polynomial{Float64}
+    @test (q, t) == original
+
+    q = x + 1
+    t = 2y
+    original = deepcopy((q, t))
+    @test_throws ArgumentError MP.MA.add!!(q, t)
+    @test (q, t) == original
+    @test_throws ArgumentError MP.MA.operate!(+, q, t)
+    @test (q, t) == original
+end
+
+@testset "In-place polynomial addition" begin
+    DP.@polyvar x y
+    for T in (Int, BigInt)
+        f = convert(DP.Polynomial{T}, x^2 + 1)
+        g = convert(DP.Polynomial{T}, 2x - 1)
+        original = deepcopy(g)
+        @test (@inferred MP.MA.add!!(f, g)) === f
+        @test f == x^2 + 2x
+        @test g == original
+        alias = SA.AlgebraElement(SA.coeffs(f), parent(f))
+        @test (@inferred MP.MA.add!!(f, alias)) === f
+        @test f == 2x^2 + 4x
+        @test g == original
+    end
+
+    f, g = x + 1, 0.5x - 1
+    originals = deepcopy((f, g))
+    result = @inferred MP.MA.add!!(f, g)
+    @test result == 1.5x
+    @test result isa DP.Polynomial{Float64}
+    @test (f, g) == originals
+
+    g = y + 1
+    originals = deepcopy((f, g))
+    @test_throws ArgumentError MP.MA.add!!(f, g)
+    @test (f, g) == originals
+end
+
+@testset "Polynomial division" begin
+    DP.@polyvar x y
+    for T in (Int, Rational{Int}, Float64, BigInt)
+        for (f, g, q, r) in
+            ((x^2 - 1, x - 1, x + 1, 0), (x^3 + 2x + 1, x^2 + 1, x, x + 1))
+            f = convert(DP.Polynomial{T}, f)
+            g = convert(DP.Polynomial{T}, g)
+            originals = deepcopy((f, g))
+            quotient, remainder = divrem(f, g)
+            @test quotient == q
+            @test remainder == r
+            @test quotient * g + remainder == f
+            @test (f, g) == originals
+        end
+    end
+    f = x^2 * y + x + 1
+    divisors = [MP.polynomial(x * y)]
+    originals = deepcopy((f, divisors))
+    quotients, remainder = divrem(f, divisors)
+    @test quotients == [x]
+    @test remainder == x + 1
+    @test (f, divisors) == originals
+
+    g = convert(DP.Polynomial{Rational{Int}}, 3x^2 + 1)
+    f = 2 * one(g)
+    original = deepcopy(g)
+    @test MP.MA.operate!(rem, f, g, MP.GeneralizedEuclideanAlgorithm()) === f
+    @test f == 2
+    @test g == original
+end
+
+@testset "Ordered term products in division" begin
+    DP.@polyvar x y z
+    f, g, t = x^2 + x + 1, x + 1, 2x
+    @test SA.term_product_style(SA.mstructure(f), SA.coeffs(f).isless) isa
+          SA.OrderedTermProduct
+    buffer = MP.MA.buffer_for(MP.MA.sub_mul, typeof(f), typeof(t), typeof(g))
+    @test buffer === nothing
+    @test MP.MA.buffered_operate!!(buffer, MP.MA.sub_mul, f, t, g) === f
+    @test f == 1 - x - x^2
+    @test g == x + 1
+
+    # Coefficient promotion still uses allocating arithmetic.
+    f, g, t = x + 1, x + 1, (1 // 2) * x
+    original = deepcopy(f)
+    result = MP.MA.operate!!(MP.MA.sub_mul, f, t, g)
+    @test result == f - t * g
+    @test result isa DP.Polynomial{Rational{Int}}
+    @test f == original
+    @test_throws ArgumentError MP.MA.operate!!(MP.MA.sub_mul, f, 2y, y + 1)
+    @test f == original
+    @test_throws ArgumentError MP.MA.operate_to!!(
+        f,
+        MP.MA.sub_mul,
+        y + 1,
+        2y,
+        y + 1,
+    )
+    @test f == original
+
+    f, g = x^2 + y, x + 1
+    originals = deepcopy((f, g))
+    @test divrem(f, g) == (x - 1, y + 1)
+    @test (f, g) == originals
+    f, divisors = x^2 + y^2 + z, [x + 1, y + 1]
+    originals = deepcopy((f, divisors))
+    q, r = divrem(f, divisors)
+    @test q == [x - 1, y - 1]
+    @test r == z + 2
+    @test (f, divisors) == originals
+    @test MP.div_multiple((x + y) * (x + 1), x + 1) == x + y
+
+    f = convert(DP.Polynomial{Rational{Int}}, x^3 + 2x + 1)
+    g = convert(DP.Polynomial{Rational{Int}}, x^2 + 1)
+    original = deepcopy(g)
+    @test MP.MA.operate!(rem, f, g, MP.GeneralizedEuclideanAlgorithm()) === f
+    @test f == x + 1
+    @test g == original
+end
+
+@testset "GCD result type promotion" begin
+    DP.@polyvar x y
+    for (a, b) in (
+        (x, y),
+        (x^2, x),
+        (6x^2, 9x),
+        (2x, y),
+        (2.0x, 3y),
+        ((2 // 3) * x, 3y),
+        ((2 + 3im) * x, 3y),
+        ((2 + 3im) * x, 3.0y),
+        (MP.term(2x, MP.monomial(y)), MP.term(3x, MP.monomial(y))),
+    )
+        for (p, q) in ((a, b), (b, a))
+            @test MP.MA.promote_operation(gcd, typeof(p), typeof(q)) ===
+                  typeof(gcd(p, q))
+            for algo in
+                (MP.GeneralizedEuclideanAlgorithm(), MP.SubresultantAlgorithm())
+                @test MP.MA.promote_operation(
+                    gcd,
+                    typeof(p),
+                    typeof(q),
+                    typeof(algo),
+                ) === typeof(gcd(p, q, algo))
+            end
+        end
+    end
+    @test MP.MA.promote_operation(gcd, typeof(x + 1), typeof(x + 1)) ===
+          DP.Polynomial{Int}
+    nested = MP.term(x, MP.monomial(y))
+    c = MP.content(nested, MP.SubresultantAlgorithm(), MP.MA.IsNotMutable())
+    @test c isa MP.AbstractMonomial
+    @test c == x
+end
+
+end # module