diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index 1b7c314..b0d2123 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -47,7 +47,7 @@ jobs:
         run: |
           using Pkg
           Pkg.add([
-              PackageSpec(name="StarAlgebras", rev="mk/non_monomial_basis"),
+              PackageSpec(name="StarAlgebras", rev="main"),
           ])
       - uses: julia-actions/julia-buildpkg@v1
       - uses: julia-actions/julia-runtest@v1
diff --git a/.github/workflows/documentation.yml b/.github/workflows/documentation.yml
index 7832ef0..eabbdf7 100644
--- a/.github/workflows/documentation.yml
+++ b/.github/workflows/documentation.yml
@@ -19,7 +19,7 @@ jobs:
         run: |
           using Pkg
           Pkg.add([
-              PackageSpec(name="StarAlgebras", rev="mk/non_monomial_basis"),
+              PackageSpec(name="StarAlgebras", rev="main"),
               PackageSpec(path=pwd()),
           ])
           Pkg.instantiate()
diff --git a/docs/src/index.md b/docs/src/index.md
index c887c5b..e901ba5 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -13,7 +13,7 @@ based on the [MultivariatePolynomials](https://github.com/JuliaAlgebra/Multivari
 ```@docs
 AbstractPolynomialBasis
 maxdegree_basis
-basis_covering_monomials
+explicit_basis_covering
 ```
 
 ## Monomial basis
diff --git a/src/MultivariateBases.jl b/src/MultivariateBases.jl
index 8e6f94b..79f692d 100644
--- a/src/MultivariateBases.jl
+++ b/src/MultivariateBases.jl
@@ -5,11 +5,33 @@ import StarAlgebras as SA
 import MultivariatePolynomials as MP
 
 export AbstractPolynomialBasis, FullBasis, SubBasis
-export maxdegree_basis, basis_covering_monomials, empty_basis
+export maxdegree_basis, explicit_basis_covering, empty_basis, monomial_index
 include("interface.jl")
 
 export AbstractMonomialIndexed, Monomial, ScaledMonomial
 include("polynomial.jl")
+MP.monomial_type(::Type{<:SA.AlgebraElement{A}}) where {A} = MP.monomial_type(A)
+function MP.polynomial_type(::Type{<:SA.AlgebraElement{A,T}}) where {A,T}
+    return MP.polynomial_type(A, T)
+end
+struct Algebra{BT,B,M} <:
+       SA.AbstractStarAlgebra{Polynomial{B,M},Polynomial{B,M}}
+    basis::BT
+end
+MP.monomial_type(::Type{<:Algebra{B}}) where {B} = MP.monomial_type(B)
+function MP.polynomial_type(::Type{<:Algebra{B}}, ::Type{T}) where {B,T}
+    return MP.polynomial_type(B, T)
+end
+SA.basis(a::Algebra) = a.basis
+
+#Base.:(==)(::Algebra{BT1,B1,M}, ::Algebra{BT2,B2,M}) where {BT1,B1,BT2,B2,M} = true
+#Base.:(==)(::Algebra, ::Algebra) = false
+
+function Base.show(io::IO, ::Algebra{BT,B}) where {BT,B}
+    ioc = IOContext(io, :limit => true, :compact => true)
+    return print(ioc, "Polynomial algebra of $B basis")
+end
+
 include("monomial.jl")
 include("scaled.jl")
 
@@ -17,7 +39,8 @@ export AbstractMultipleOrthogonal,
     ProbabilistsHermite, PhysicistsHermite, Laguerre
 export AbstractGegenbauer,
     Legendre, Chebyshev, ChebyshevFirstKind, ChebyshevSecondKind
-export generators,
+export algebra_element,
+    sparse_coefficients,
     univariate_orthogonal_basis,
     reccurence_first_coef,
     reccurence_second_coef,
@@ -34,4 +57,19 @@ include("legendre.jl")
 include("chebyshev.jl")
 include("quotient.jl")
 
+function algebra(
+    basis::Union{QuotientBasis{Polynomial{B,M}},FullBasis{B,M},SubBasis{B,M}},
+) where {B,M}
+    return Algebra{typeof(basis),B,M}(basis)
+end
+
+function MA.promote_operation(
+    ::typeof(algebra),
+    BT::Type{
+        <:Union{QuotientBasis{Polynomial{B,M}},FullBasis{B,M},SubBasis{B,M}},
+    },
+) where {B,M}
+    return Algebra{BT,B,M}
+end
+
 end # module
diff --git a/src/chebyshev.jl b/src/chebyshev.jl
index 27126fd..fbd4285 100644
--- a/src/chebyshev.jl
+++ b/src/chebyshev.jl
@@ -1,13 +1,13 @@
 abstract type AbstractChebyshev <: AbstractGegenbauer end
 
 _promote_div(::Type{I}) where {I<:Integer} = Rational{I}
-_promote_div(::Type{T}) where {T} = MA.promote_operation(/, T, Int)
+_promote_div(::Type{T}) where {T<:Number} = MA.promote_operation(/, T, Int)
+# Could be for instance `MathOptInterface.ScalarAffineFunction{Float64}`
+# which does not support division with `Int`
+_promote_div(::Type{F}) where {F} = F
 
-function MP.polynomial_type(
-    ::Type{Polynomial{B,M}},
-    ::Type{T},
-) where {B<:AbstractChebyshev,M,T}
-    return MP.polynomial_type(M, _promote_div(T))
+function _promote_coef(::Type{T}, ::Type{<:AbstractChebyshev}) where {T}
+    return _promote_div(T)
 end
 
 reccurence_first_coef(::Type{<:AbstractChebyshev}, degree) = 2
@@ -26,23 +26,28 @@ const Chebyshev = ChebyshevFirstKind
 # T_n * T_m = T_{n + m} / 2 + T_{|n - m|} / 2
 function (::Mul{Chebyshev})(a::MP.AbstractMonomial, b::MP.AbstractMonomial)
     terms = [MP.term(1 // 1, MP.constant_monomial(a * b))]
+    vars_a = MP.variables(a)
+    var_state_a = iterate(vars_a)
     vars_b = MP.variables(b)
     var_state_b = iterate(vars_b)
-    for var_a in MP.variables(a)
-        if isnothing(var_state_b)
-            break
-        end
-        var_b, state_b = var_state_b
-        if var_b > var_a
-            var_state_b = iterate(vars_b, state_b)
+    while !isnothing(var_state_a) || !isnothing(var_state_b)
+        if isnothing(var_state_a) ||
+           (!isnothing(var_state_b) && var_state_b[1] > var_state_a[1])
+            var_b, state_b = var_state_b
             for i in eachindex(terms)
                 terms[i] = MA.mul!!(terms[i], var_b^MP.degree(b, var_b))
             end
-        elseif var_a > var_b
+            var_state_b = iterate(vars_b, state_b)
+        elseif isnothing(var_state_b) ||
+               (!isnothing(var_state_a) && var_state_a[1] > var_state_b[1])
+            var_a, state_a = var_state_a
             for i in eachindex(terms)
                 terms[i] = MA.mul!!(terms[i], var_a^MP.degree(a, var_a))
             end
+            var_state_a = iterate(vars_a, state_a)
         else
+            var_a, state_a = var_state_a
+            var_b, state_b = var_state_b
             d_a = MP.degree(a, var_a)
             d_b = MP.degree(b, var_b)
             I = eachindex(terms)
@@ -53,9 +58,62 @@ function (::Mul{Chebyshev})(a::MP.AbstractMonomial, b::MP.AbstractMonomial)
                 α = MA.copy_if_mutable(MP.coefficient(terms[i]))
                 push!(terms, MP.term(α, mono))
             end
+            var_state_a = iterate(vars_a, state_a)
+            var_state_b = iterate(vars_b, state_b)
         end
     end
-    return MP.polynomial!(terms)
+    return sparse_coefficients(MP.polynomial!(terms))
+end
+
+function _add_mul_scalar_vector!(res, ::SubBasis, scalar, vector)
+    return MA.operate!(MA.add_mul, res, scalar, vector)
+end
+
+function _add_mul_scalar_vector!(res, ::FullBasis, scalar, vector)
+    for (k, v) in SA.nonzero_pairs(vector)
+        SA.unsafe_push!(res, k, scalar * v)
+    end
+end
+
+function SA.coeffs!(
+    res,
+    cfs,
+    source::MonomialIndexedBasis{<:AbstractMultipleOrthogonal},
+    target::MonomialIndexedBasis{Monomial},
+)
+    MA.operate!(zero, res)
+    for (k, v) in SA.nonzero_pairs(cfs)
+        _add_mul_scalar_vector!(res, target, v, SA.coeffs(source[k], target))
+    end
+    MA.operate!(SA.canonical, res)
+    return res
+end
+
+function SA.coeffs(
+    cfs,
+    source::SubBasis{Monomial},
+    target::FullBasis{Chebyshev},
+)
+    sub = explicit_basis_covering(target, source)
+    # Need to make A square so that it's UpperTriangular
+    extended = SubBasis{Monomial}(sub.monomials)
+    A = zeros(Float64, length(extended), length(sub))
+    for (i, cheby) in enumerate(sub)
+        A[:, i] = SA.coeffs(algebra_element(cheby), extended)
+    end
+    ext = SA.coeffs(algebra_element(cfs, source), extended)
+    return SA.SparseCoefficients(
+        sub.monomials,
+        LinearAlgebra.UpperTriangular(A) \ ext,
+    )
+end
+
+function SA.coeffs(cfs, ::FullBasis{Monomial}, target::FullBasis{Chebyshev})
+    return SA.coeffs(
+        SA.values(cfs),
+        SubBasis{Monomial}(collect(SA.keys(cfs))),
+        target,
+    )
 end
 
 function degree_one_univariate_polynomial(
diff --git a/src/hermite.jl b/src/hermite.jl
index 47c9f71..a95353e 100644
--- a/src/hermite.jl
+++ b/src/hermite.jl
@@ -1,11 +1,6 @@
 abstract type AbstractHermite <: AbstractMultipleOrthogonal end
 
-function MP.polynomial_type(
-    ::Type{Polynomial{B,M}},
-    ::Type{T},
-) where {B<:AbstractHermite,M,T}
-    return MP.polynomial_type(M, float(T))
-end
+_promote_coef(::Type{T}, ::Type{<:AbstractHermite}) where {T} = _float(T)
 
 even_odd_separated(::Type{<:AbstractHermite}) = true
 
diff --git a/src/interface.jl b/src/interface.jl
index 4626960..2fe5c04 100644
--- a/src/interface.jl
+++ b/src/interface.jl
@@ -27,11 +27,12 @@ Return the explicit version of `basis`generating all polynomials of degree up to
 """
 function maxdegree_basis end
 
+# TODO remove, not needed anymore
 """
-    basis_covering_monomials(basis::StarAlgebras.AbstractBasis, monos::AbstractVector{<:AbstractMonomial})
+    explicit_basis_covering(basis::StarAlgebras.AbstractBasis, target::StarAlgebras.ExplicitBasis)
 
-Return the minimal basis of type `B` that can generate all polynomials of the
-monomial basis generated by `monos`.
+Return the minimal basis of type `B` that can generate all polynomials generated by
+the basis `target`.
 
 ## Examples
 
@@ -44,8 +45,8 @@ julia> using DynamicPolynomials
 julia> @polyvar x
 (x,)
 
-julia> basis_covering_monomials(FullBasis{Chebyshev,typeof(x^2)}(), [x^2, x^4])
+julia> explicit_basis_covering(FullBasis{Chebyshev,typeof(x^2)}(), SubBasis{Monomial}([x^2, x^4]))
 SubBasis{ChebyshevFirstKind}([1, x², x⁴])
 ```
 """
-function basis_covering_monomials end
+function explicit_basis_covering end
diff --git a/src/laguerre.jl b/src/laguerre.jl
index 605bb78..c279409 100644
--- a/src/laguerre.jl
+++ b/src/laguerre.jl
@@ -7,12 +7,7 @@ struct Laguerre <: AbstractMultipleOrthogonal end
 
 # TODO implement multiplication with https://www.jstor.org/stable/2002985
 
-function MP.polynomial_type(
-    ::Type{Polynomial{Laguerre,M}},
-    ::Type{T},
-) where {M,T}
-    return MP.polynomial_type(M, float(T))
-end
+_promote_coef(::Type{T}, ::Type{<:Laguerre}) where {T} = _float(T)
 
 even_odd_separated(::Type{Laguerre}) = false
 
diff --git a/src/legendre.jl b/src/legendre.jl
index 1301e8c..0595125 100644
--- a/src/legendre.jl
+++ b/src/legendre.jl
@@ -15,12 +15,7 @@ Orthogonal polynomial with respect to the univariate weight function ``w(x) = 1`
 """
 struct Legendre <: AbstractGegenbauer end
 
-function MP.polynomial_type(
-    ::Type{Polynomial{Legendre,M}},
-    ::Type{T},
-) where {M,T}
-    return MP.polynomial_type(M, float(T))
-end
+_promote_coef(::Type{T}, ::Type{<:Legendre}) where {T} = _float(T)
 
 reccurence_first_coef(::Type{Legendre}, degree) = (2degree - 1)
 reccurence_third_coef(::Type{Legendre}, degree) = -(degree - 1)
diff --git a/src/monomial.jl b/src/monomial.jl
index 6abc7bc..42b92a6 100644
--- a/src/monomial.jl
+++ b/src/monomial.jl
@@ -64,6 +64,8 @@ function Base.getindex(basis::SubBasis{B,M}, value::Polynomial{B,M}) where {B,M}
     return mono
 end
 
+const MonomialIndexedBasis{B,M} = Union{SubBasis{B,M},FullBasis{B,M}}
+
 MP.monomial_type(::Type{<:SubBasis{B,M}}) where {B,M} = M
 
 # The `i`th index of output is the index of occurence of `x[i]` in `y`,
@@ -116,8 +118,28 @@ end
 
 Base.:(==)(a::SubBasis, b::SubBasis) = a.monomials == b.monomials
 
-function _algebra(basis::Union{FullBasis{B,M},SubBasis{B,M}}) where {B,M}
-    return SA.StarAlgebra(Polynomial{B}(MP.constant_monomial(M)), basis)
+function algebra_type(::Type{BT}) where {B,M,BT<:MonomialIndexedBasis{B,M}}
+    return Algebra{BT,B,M}
+end
+
+implicit_basis(::SubBasis{B,M}) where {B,M} = FullBasis{B,M}()
+implicit_basis(basis::FullBasis) = basis
+
+function MA.promote_operation(
+    ::typeof(implicit_basis),
+    ::Type{<:Union{FullBasis{B,M},SubBasis{B,M}}},
+) where {B,M}
+    return FullBasis{B,M}
+end
+
+function _explicit_basis(coeffs, ::FullBasis{B}) where {B}
+    return SubBasis{B}(SA.keys(coeffs))
+end
+
+_explicit_basis(_, basis::SubBasis) = basis
+
+function explicit_basis(a::SA.AlgebraElement)
+    return _explicit_basis(SA.coeffs(a), SA.basis(a))
 end
 
 function explicit_basis_type(::Type{FullBasis{B,M}}) where {B,M}
@@ -138,47 +160,53 @@ function maxdegree_basis(
     return unsafe_basis(B, MP.monomials(variables, 0:maxdegree))
 end
 
-function MP.polynomial(coefs::Vector, basis::SubBasis)
-    return MP.polynomial(Base.Fix1(getindex, coefs), basis)
+MP.variables(c::SA.AbstractCoefficients) = MP.variables(SA.keys(c))
+
+function sparse_coefficients(p::MP.AbstractPolynomial)
+    return SA.SparseCoefficients(MP.monomials(p), MP.coefficients(p))
 end
 
-function MP.polynomial(f::Function, basis::SubBasis)
-    if isempty(basis)
-        return zero(MP.polynomial_type(basis))
-    else
-        return MP.polynomial(
-            mapreduce(
-                ip -> f(ip[1]) * MP.polynomial(ip[2]),
-                MA.add!!,
-                enumerate(basis),
-            ),
-        )
-    end
+function sparse_coefficients(t::MP.AbstractTermLike)
+    return SA.SparseCoefficients((MP.monomial(t),), (MP.coefficient(t),))
 end
 
-_convert(::Type{P}, p) where {P} = convert(P, p)
-_convert(::Type{P}, ::MA.Zero) where {P} = zero(P)
-
-function MP.polynomial(Q::AbstractMatrix, basis::SubBasis, ::Type{T}) where {T}
-    n = length(basis)
-    @assert size(Q) == (n, n)
-    PT = MP.polynomial_type(eltype(basis), T)
-    return _convert(
-        PT,
-        mapreduce(
-            row -> begin
-                adjoint(MP.polynomial(basis[row])) * mapreduce(
-                    col -> Q[row, col] * MP.polynomial(basis[col]),
-                    MA.add!!,
-                    1:n;
-                    init = MA.Zero(),
-                )
-            end,
-            MA.add!!,
-            1:n;
-            init = MA.Zero(),
+function MA.promote_operation(
+    ::typeof(sparse_coefficients),
+    ::Type{P},
+) where {P<:MP.AbstractPolynomialLike}
+    M = MP.monomial_type(P)
+    T = MP.coefficient_type(P)
+    return SA.SparseCoefficients{M,T,MP.monomial_vector_type(M),Vector{T}}
+end
+
+function algebra_element(p::MP.AbstractPolynomialLike)
+    return algebra_element(
+        sparse_coefficients(p),
+        FullBasis{Monomial,MP.monomial_type(p)}(),
+    )
+end
+
+function algebra_element(f::Function, basis::SubBasis)
+    return algebra_element(map(f, eachindex(basis)), basis)
+end
+
+function constant_algebra_element(
+    ::Type{FullBasis{B,M}},
+    ::Type{T},
+) where {B,M,T}
+    return algebra_element(
+        sparse_coefficients(
+            MP.polynomial(MP.term(one(T), MP.constant_monomial(M))),
         ),
-    )::PT
+        FullBasis{B,M}(),
+    )
+end
+
+function constant_algebra_element(
+    ::Type{<:SubBasis{B,M}},
+    ::Type{T},
+) where {B,M,T}
+    return algebra_element([one(T)], SubBasis{B}([MP.constant_monomial(M)]))
 end
 
 function _show(io::IO, mime::MIME, basis::SubBasis{B}) where {B}
@@ -215,11 +243,15 @@ end
 
 abstract type AbstractMonomial <: AbstractMonomialIndexed end
 
-function basis_covering_monomials(
+function explicit_basis_covering(
     ::FullBasis{B},
-    monos::AbstractVector,
+    target::SubBasis{<:AbstractMonomial},
 ) where {B<:AbstractMonomial}
-    return SubBasis{B}(monos)
+    return SubBasis{B}(target.monomials)
+end
+
+function Base.adjoint(p::Polynomial{B}) where {B<:AbstractMonomialIndexed}
+    return Polynomial{B}(adjoint(p.monomial))
 end
 
 """
@@ -235,9 +267,11 @@ one get ths [`ScaledMonomial`](@ref).
 """
 struct Monomial <: AbstractMonomial end
 
-(::Mul{Monomial})(a::MP.AbstractMonomial, b::MP.AbstractMonomial) = a * b
+function (::Mul{Monomial})(a::MP.AbstractMonomial, b::MP.AbstractMonomial)
+    return sparse_coefficients(a * b)
+end
 
-MP.polynomial(p::Polynomial{Monomial}) = MP.polynomial(p.monomial)
+SA.coeffs(p::Polynomial{Monomial}, ::FullBasis{Monomial}) = p.monomial
 
 function MP.polynomial_type(
     ::Union{SubBasis{B,M},Type{<:SubBasis{B,M}}},
@@ -246,15 +280,8 @@ function MP.polynomial_type(
     return MP.polynomial_type(FullBasis{B,M}, T)
 end
 
-function MP.polynomial_type(::Type{FullBasis{B,M}}, ::Type{T}) where {B,M,T}
-    return MP.polynomial_type(Polynomial{B,M}, T)
-end
-
-function MP.polynomial_type(
-    ::Type{Polynomial{Monomial,M}},
-    ::Type{T},
-) where {M,T}
-    return MP.polynomial_type(M, T)
+function MP.polynomial_type(::Type{Polynomial{B,M}}, ::Type{T}) where {B,M,T}
+    return MP.polynomial_type(FullBasis{B,M}, T)
 end
 
 function MP.polynomial(f::Function, mb::SubBasis{Monomial})
@@ -276,6 +303,52 @@ function MP.coefficients(p::MP.AbstractPolynomialLike, ::FullBasis{Monomial})
     return p
 end
 
+function _assert_constant(α) end
+
+function _assert_constant(
+    x::Union{Polynomial,SA.AlgebraElement,MP.AbstractPolynomialLike},
+)
+    return error("Expected constant element, got type `$(typeof(x))`")
+end
+
+#function MA.operate!(::SA.UnsafeAddMul{<:Mul{Monomial}}, p::MP.AbstractPolynomial, args::Vararg{Any,N}) where {N}
+#    return MA.operate!(MA.add_mul, p, args...)
+#end
+
+#function MA.operate!(
+#    ::SA.UnsafeAddMul{Mul{B}},
+#    res::SA.AbstractCoefficients,
+#    v::SA.AbstractCoefficients,
+#    w::SA.AbstractCoefficients,
+#) where {B<:MB.AbstractMonomial}
+#    for (kv, a) in nonzero_pairs(v)
+#        for (kw, b) in nonzero_pairs(w)
+#            SA.unsafe_push!(res, kv * kw, a * b)
+#            c = ms.structure(kv, kw)
+#            for (k, v) in nonzero_pairs(c)
+#            end
+#        end
+#    end
+#    return res
+#end
+
+#function MA.operate!(op::SA.UnsafeAddMul{Mul{B}}, a::SA.AbstractCoefficients, α, b::SA.AbstractCoefficients) where {B<:MB.AbstractMonomial}
+#    for
+#    MA.operate!(op, a, α, Polynomial{Monomial}(x.monomial * y.monomial * z.monomial))
+#    return a
+#end
+
+function MA.operate!(
+    ::SA.UnsafeAddMul{typeof(*)},
+    a::SA.AlgebraElement{<:Algebra{<:MonomialIndexedBasis,Monomial}},
+    α,
+    x::Polynomial{Monomial},
+)
+    _assert_constant(α)
+    SA.unsafe_push!(a, x, α)
+    return a
+end
+
 # Overload some of the `MP` interface for convenience
 MP.mindegree(basis::SubBasis{Monomial}) = MP.mindegree(basis.monomials)
 MP.maxdegree(basis::SubBasis) = MP.maxdegree(basis.monomials)
@@ -289,3 +362,21 @@ end
 function MP.extdegree(basis::SubBasis{Monomial}, v)
     return MP.extdegree(basis.monomials, v)
 end
+
+_promote_coef(::Type{T}, ::Type{Monomial}) where {T} = T
+
+function MP.polynomial_type(::Type{FullBasis{B,M}}, ::Type{T}) where {T,B,M}
+    return MP.polynomial_type(M, _promote_coef(T, B))
+end
+
+# Adapted from SA to incorporate `_promote_coef`
+function SA.coeffs(
+    cfs,
+    source::MonomialIndexedBasis{B},
+    target::MonomialIndexedBasis{Monomial},
+) where {B}
+    source === target && return cfs
+    source == target && return cfs
+    res = SA.zero_coeffs(_promote_coef(valtype(cfs), B), target)
+    return SA.coeffs!(res, cfs, source, target)
+end
diff --git a/src/orthogonal.jl b/src/orthogonal.jl
index bc44905..9307c0e 100644
--- a/src/orthogonal.jl
+++ b/src/orthogonal.jl
@@ -105,12 +105,12 @@ function univariate_orthogonal_basis(
     end
 end
 
-function basis_covering_monomials(
+function explicit_basis_covering(
     ::FullBasis{B,M},
-    monos::AbstractVector,
+    monos::SubBasis{<:AbstractMonomial,M},
 ) where {B<:AbstractMultipleOrthogonal,M}
-    to_add = collect(monos)
-    m = Set{M}(monos)
+    to_add = collect(monos.monomials)
+    m = Set{M}(to_add)
     while !isempty(to_add)
         mono = pop!(to_add)
         for v in MP.variables(mono)
@@ -172,15 +172,31 @@ function MP.coefficients(
     p,
     basis::SubBasis{B,M},
 ) where {B<:AbstractMultipleOrthogonal,M}
+    poly_p = MP.polynomial(p)
     return map(basis) do el
-        q = MP.polynomial(el)
-        return LinearAlgebra.dot(p, q, B) / LinearAlgebra.dot(q, q, B)
+        q = SA.coeffs(el, FullBasis{Monomial,M}())
+        poly_q = MP.polynomial(q)
+        return LinearAlgebra.dot(poly_p, poly_q, B) /
+               LinearAlgebra.dot(poly_q, poly_q, B)
     end
 end
 
-function MP.polynomial(p::Polynomial{B}) where {B<:AbstractMultipleOrthogonal}
-    return prod(
-        univariate_orthogonal_basis(B, var, deg)[deg+1] for
-        (var, deg) in MP.powers(p.monomial)
+function SA.coeffs(
+    p::Polynomial{B},
+    basis::SubBasis{Monomial},
+) where {B<:AbstractMultipleOrthogonal}
+    full = implicit_basis(basis)
+    return SA.coeffs(algebra_element(SA.coeffs(p, full), full), basis)
+end
+
+function SA.coeffs(
+    p::Polynomial{B},
+    ::FullBasis{Monomial},
+) where {B<:AbstractMultipleOrthogonal}
+    return sparse_coefficients(
+        prod(
+            univariate_orthogonal_basis(B, var, deg)[deg+1] for
+            (var, deg) in MP.powers(p.monomial)
+        ),
     )
 end
diff --git a/src/polynomial.jl b/src/polynomial.jl
index 2671314..d1f718b 100644
--- a/src/polynomial.jl
+++ b/src/polynomial.jl
@@ -1,4 +1,5 @@
 # TODO Add to MultivariatePolynomials
+MP.variables(p::SA.AlgebraElement) = MP.variables(SA.basis(p))
 Base.keytype(p::MP.AbstractPolynomialLike) = MP.monomial_type(p)
 Base.valtype(p::MP.AbstractPolynomialLike) = MP.coefficient_type(p)
 #Base.keys(p::MP.AbstractPolynomial) = MP.monomials(p)
@@ -6,6 +7,9 @@ SA.nonzero_pairs(p::MP.AbstractPolynomialLike) = MP.terms(p)
 function Base.similar(p::PT, ::Type{T}) where {PT<:MP.AbstractPolynomial,T}
     return convert(MP.similar_type(PT, T), copy(p)) # Missing the `copy` in MP
 end
+function Base.getindex(p::MP.AbstractPolynomialLike, mono::MP.AbstractMonomial)
+    return MP.coefficient(p, mono)
+end
 Base.iterate(t::MP.Term) = iterate(t, 1)
 function Base.iterate(t::MP.Term, state)
     if state == 1
@@ -16,11 +20,14 @@ function Base.iterate(t::MP.Term, state)
         return nothing
     end
 end
+function SA.unsafe_push!(p::MP.AbstractPolynomial, mono::MP.AbstractMonomial, α)
+    return MA.operate!(MA.add_mul, p, α, mono)
+end
 function MA.operate!(
     ::SA.UnsafeAddMul{typeof(*)},
     mc::MP.AbstractPolynomial,
     val,
-    c::MP.AbstractPolynomial,
+    c::MP.AbstractPolynomialLike,
 )
     return MA.operate!(MA.add_mul, mc, val, c)
 end
@@ -41,6 +48,18 @@ struct Polynomial{B<:AbstractMonomialIndexed,M<:MP.AbstractMonomial}
     end
 end
 
+function Base.hash(p::Polynomial{B}, u::UInt) where {B}
+    return hash(B, hash(p.monomial, u))
+end
+
+function Base.isequal(p::Polynomial{B}, q::Polynomial{B}) where {B}
+    return isequal(p.monomial, q.monomial)
+end
+
+# Needed for `BoundsError`
+Base.iterate(p::Polynomial) = p, nothing
+Base.iterate(::Polynomial, ::Nothing) = nothing
+
 function Polynomial{B}(v::MP.AbstractVariable) where {B}
     return Polynomial{B}(MP.monomial(v))
 end
@@ -52,31 +71,126 @@ end
 MP.variables(p::Polynomial) = MP.variables(p.monomial)
 MP.nvariables(p::Polynomial) = MP.nvariables(p.monomial)
 
-function _algebra_element(p, ::Type{B}) where {B}
-    return SA.AlgebraElement(p, _algebra(FullBasis{B,MP.monomial_type(p)}()))
+MP.monomial_type(::Type{<:SA.SparseCoefficients{K}}) where {K} = K
+MP.polynomial(p::Polynomial) = MP.polynomial(algebra_element(p))
+
+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,MP.monomial_type(typeof(p))}(),
+    )
+end
+
+function algebra_element(p::Polynomial{B,M}) where {B,M}
+    return _algebra_element(p.monomial, B)
 end
 
 function Base.:*(a::Polynomial{B}, b::Polynomial{B}) where {B}
-    return _algebra_element(Mul{B}()(a.monomial, b.monomial), B)
+    return algebra_element(
+        Mul{B}()(a.monomial, b.monomial),
+        FullBasis{B,promote_type(typeof(a.monomial), typeof(b.monomial))}(),
+    )
+end
+
+function Base.:*(a::Polynomial{B}, b::SA.AlgebraElement) where {B}
+    return _algebra_element(a) * b
 end
 
 function _show(io::IO, mime::MIME, p::Polynomial{B}) where {B}
-    print(io, B)
-    print(io, "(")
-    show(io, mime, p.monomial)
-    return print(io, ")")
+    if B != Monomial
+        print(io, B)
+        print(io, "(")
+    end
+    print(io, SA.trim_LaTeX(mime, sprint(show, mime, p.monomial)))
+    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, p::Polynomial)
-    return show(io, MIME"text/plain"(), p)
+
+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)
+
 function Base.zero(::Type{Polynomial{B,M}}) where {B,M}
     return _algebra_element(zero(MP.polynomial_type(M, Rational{Int})), B)
 end
 
 Base.zero(p::Polynomial) = zero(typeof(p))
 
+function convert_basis(basis::SA.AbstractBasis, p::MP.AbstractPolynomialLike)
+    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
+
 struct Mul{B<:AbstractMonomialIndexed} <: SA.MultiplicativeStructure end
+
+function MA.operate_to!(
+    p::MP.AbstractPolynomial,
+    op::Mul,
+    args::Vararg{MP.AbstractPolynomialLike,N},
+) where {N}
+    MA.operate!(zero, p)
+    MA.operate!(SA.UnsafeAddMul(op), p, args...)
+    MA.operate!(SA.canonical, p)
+    return p
+end
+
+function MP.polynomial(a::SA.AbstractCoefficients)
+    return MP.polynomial(collect(SA.values(a)), collect(SA.keys(a)))
+end
+
+function MP.polynomial(a::SA.AlgebraElement)
+    return MP.polynomial(
+        SA.coeffs(a, FullBasis{Monomial,MP.monomial_type(typeof(a))}()),
+    )
+end
+
+function Base.isapprox(
+    p::MP.AbstractPolynomialLike,
+    a::SA.AlgebraElement;
+    kws...,
+)
+    return isapprox(p, MP.polynomial(a); kws...)
+end
+
+function Base.isapprox(a::SA.AlgebraElement, b::SA.AlgebraElement; kws...)
+    return isapprox(MP.polynomial(a), b; kws...)
+end
+
+function Base.isapprox(
+    a::SA.AlgebraElement,
+    p::MP.AbstractPolynomialLike;
+    kws...,
+)
+    return isapprox(p, a; kws...)
+end
+
+function Base.isapprox(a::SA.AlgebraElement, α::Number; kws...)
+    return isapprox(
+        a,
+        α * constant_algebra_element(typeof(SA.basis(a)), typeof(α));
+        kws...,
+    )
+end
diff --git a/src/quotient.jl b/src/quotient.jl
index 48c970b..54295a8 100644
--- a/src/quotient.jl
+++ b/src/quotient.jl
@@ -3,8 +3,36 @@ struct QuotientBasis{T,I,B<:SA.AbstractBasis{T,I},D} <: SA.ExplicitBasis{T,I}
     divisor::D
 end
 
+implicit_basis(basis::QuotientBasis) = implicit_basis(basis.basis)
+
+Base.iterate(basis::QuotientBasis) = iterate(basis.basis)
+Base.iterate(basis::QuotientBasis, s) = iterate(basis.basis, s)
 Base.length(basis::QuotientBasis) = length(basis.basis)
 
+_object(basis::QuotientBasis) = _object(basis.basis)
+
 function MP.coefficients(p, basis::QuotientBasis)
     return MP.coefficients(rem(p, basis.divisor), basis.basis)
 end
+
+function SA.coeffs(p, ::FullBasis{Monomial}, basis::QuotientBasis)
+    return MP.coefficients(MP.polynomial(p), basis)
+end
+
+function SA.coeffs(coeffs, sub::SubBasis{Monomial}, basis::QuotientBasis)
+    return SA.coeffs(MP.polynomial(coeffs, sub.monomials), parent(sub), basis)
+end
+
+function SA.adjoint_coeffs(
+    coeffs,
+    src::SubBasis{Monomial},
+    dest::QuotientBasis{<:Polynomial{Monomial}},
+)
+    return map(src.monomials) do mono
+        return sum(
+            MP.coefficient(t) *
+            coeffs[dest.basis[Polynomial{Monomial}(MP.monomial(t))]] for
+            t in MP.terms(rem(mono, dest.divisor))
+        )
+    end
+end
diff --git a/src/scaled.jl b/src/scaled.jl
index 6c3630c..60403d9 100644
--- a/src/scaled.jl
+++ b/src/scaled.jl
@@ -1,5 +1,5 @@
 """
-    struct Scaled <: AbstractMonomial end
+    struct ScaledMonomial <: AbstractMonomial end
 
 *Scaled monomial basis* (see [Section 3.1.5, BPT12]) with the monomials of the vector `monomials`.
 Given a monomial ``x^\\alpha = x_1^{\\alpha_1} \\cdots x_n^{\\alpha_n}`` of degree ``d = \\sum_{i=1}^n \\alpha_i``,
@@ -33,16 +33,28 @@ Foundations of Computational Mathematics 7.2 (2007): 229-244.
 """
 struct ScaledMonomial <: AbstractMonomial end
 
-function MP.polynomial(p::Polynomial{ScaledMonomial})
-    return scaling(p.monomial) * p.monomial
+function (::Mul{ScaledMonomial})(a::MP.AbstractMonomial, b::MP.AbstractMonomial)
+    mono = a * b
+    α = prod(
+        MP.variables(mono);
+        init = inv(binomial(MP.degree(mono), MP.degree(a))),
+    ) do v
+        return binomial(MP.degree(mono, v), MP.degree(a, v))
+    end
+    return sparse_coefficients(MP.term(√α, mono))
 end
 
-function MP.polynomial_type(
-    ::Type{FullBasis{ScaledMonomial,M}},
-    T::Type,
-) where {M}
-    return MP.polynomial_type(M, float(T))
+function SA.coeffs(p::Polynomial{ScaledMonomial}, ::FullBasis{Monomial})
+    return scaling(p.monomial) * p.monomial
 end
+
+_float(::Type{T}) where {T<:Number} = float(T)
+# Could be for instance `MathOptInterface.ScalarAffineFunction{Float64}`
+# which does not implement `float`
+_float(::Type{F}) where {F} = F
+
+_promote_coef(::Type{T}, ::Type{ScaledMonomial}) where {T} = _float(T)
+
 function MP.polynomial(f::Function, basis::SubBasis{ScaledMonomial})
     return MP.polynomial(
         i -> scaling(basis.monomials[i]) * f(i),
@@ -57,32 +69,66 @@ function Base.promote_rule(
     return SubBasis{Monomial,M,V}
 end
 
-function change_basis(
-    Q::AbstractMatrix,
-    basis::SubBasis{ScaledMonomial,M},
-    ::Type{Monomial},
-) where {M}
-    n = length(basis)
-    scalings = map(scaling, basis.monomials)
-    scaled_Q = [Q[i, j] * scalings[i] * scalings[j] for i in 1:n, j in 1:n]
-    return scaled_Q, SubBasis{Monomial}(basis.monomials)
-end
-
-function MP.polynomial(
-    Q::AbstractMatrix,
-    basis::SubBasis{ScaledMonomial,M,V},
-    ::Type{T},
-) where {M,V<:AbstractVector{M},T}
-    return MP.polynomial(change_basis(Q, basis, Monomial)..., T)
-end
-
 function scaling(m::MP.AbstractMonomial)
     return √(factorial(MP.degree(m)) / prod(factorial, MP.exponents(m)))
 end
 unscale_coef(t::MP.AbstractTerm) = MP.coefficient(t) / scaling(MP.monomial(t))
+function SA.coeffs(
+    t::MP.AbstractTermLike,
+    ::FullBasis{ScaledMonomial},
+    ::FullBasis{Monomial},
+)
+    mono = MP.monomial(t)
+    return MP.term(mono * MP.coefficient(t), mono)
+end
 function MP.coefficients(p, ::FullBasis{ScaledMonomial})
     return unscale_coef.(MP.terms(p))
 end
 function MP.coefficients(p, basis::SubBasis{ScaledMonomial})
     return MP.coefficients(p, basis.monomials) ./ scaling.(MP.monomials(p))
 end
+function SA.coeffs(
+    p::MP.AbstractPolynomialLike,
+    ::FullBasis{Monomial},
+    basis::FullBasis{ScaledMonomial},
+)
+    return MP.polynomial(MP.coefficients(p, basis), MP.monomials(p))
+end
+
+function SA.coeffs!(
+    res,
+    cfs,
+    source::MonomialIndexedBasis{ScaledMonomial},
+    target::MonomialIndexedBasis{Monomial},
+)
+    MA.operate!(zero, res)
+    for (k, v) in SA.nonzero_pairs(cfs)
+        mono = source[k].monomial
+        SA.unsafe_push!(
+            res,
+            target[Polynomial{Monomial}(mono)],
+            v * scaling(mono),
+        )
+    end
+    MA.operate!(SA.canonical, res)
+    return res
+end
+
+function SA.coeffs!(
+    res,
+    cfs,
+    source::MonomialIndexedBasis{Monomial},
+    target::MonomialIndexedBasis{ScaledMonomial},
+)
+    MA.operate!(zero, res)
+    for (k, v) in SA.nonzero_pairs(cfs)
+        mono = source[k].monomial
+        SA.unsafe_push!(
+            res,
+            target[Polynomial{ScaledMonomial}(mono)],
+            v / scaling(mono),
+        )
+    end
+    MA.operate!(SA.canonical, res)
+    return res
+end
diff --git a/test/Project.toml b/test/Project.toml
index 5fdc69e..303677f 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -2,5 +2,6 @@
 DynamicPolynomials = "7c1d4256-1411-5781-91ec-d7bc3513ac07"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MultivariateBases = "be282fd4-ad43-11e9-1d11-8bd9d7e43378"
+MutableArithmetics = "d8a4904e-b15c-11e9-3269-09a3773c0cb0"
 StarAlgebras = "0c0c59c1-dc5f-42e9-9a8b-b5dc384a6cd1"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
diff --git a/test/chebyshev.jl b/test/chebyshev.jl
index 2d0a312..68ebdf1 100644
--- a/test/chebyshev.jl
+++ b/test/chebyshev.jl
@@ -6,9 +6,12 @@ using DynamicPolynomials
     @polyvar x
     a = MB.Polynomial{MB.Chebyshev}(x)
     b = a * a
-    @test b.coeffs == 1 // 2 + 1 // 2 * x^2
+    @test b.coeffs == MB.sparse_coefficients(1 // 2 + 1 // 2 * x^2)
     c = b * b
-    @test c.coeffs == 3 // 8 + 1 // 2 * x^2 + 1 // 8 * x^4
+    @test c.coeffs ==
+          MB.sparse_coefficients(3 // 8 + 1 // 2 * x^2 + 1 // 8 * x^4)
+    @test a * MB.Polynomial{MB.Chebyshev}(constant_monomial(typeof(x))) ==
+          a * MB.Polynomial{MB.Chebyshev}(x^0)
 end
 
 @testset "Orthogonal" begin
diff --git a/test/monomial.jl b/test/monomial.jl
index 546ba81..89bb4fb 100644
--- a/test/monomial.jl
+++ b/test/monomial.jl
@@ -8,8 +8,8 @@ using DynamicPolynomials
     @polyvar x
     a = MB.Polynomial{MB.Monomial}(x)
     b = a * a
-    @test b.coeffs == x^2
-    @test typeof(b.coeffs) == typeof(x^2)
+    @test b.coeffs == sparse_coefficients(x^2)
+    @test typeof(b.coeffs) == typeof(sparse_coefficients(x^2))
     @test MB.Polynomial{MB.Monomial}(x^2) == MB.Polynomial{MB.Monomial}(x^2)
     @test MB.Polynomial{MB.Monomial}(x^3) != MB.Polynomial{MB.Monomial}(x^2)
 end
@@ -65,6 +65,18 @@ end
     @test basis.monomials == [y^2, x * y, x^2]
     @test I1 == [1, 0, 2]
     @test I2 == [1, 2, 0]
+    for i in eachindex(basis)
+        mono = basis.monomials[i]
+        @test monomial_index(basis, mono) == i
+        for (I, b) in [(I1, basis1), (I2, basis2)]
+            idx = monomial_index(b, basis.monomials[i])
+            if iszero(I[i])
+                @test isnothing(idx)
+            else
+                @test idx == I[i]
+            end
+        end
+    end
 end
 
 @testset "API degree = $degree" for degree in 0:3
diff --git a/test/quotient.jl b/test/quotient.jl
index 417edf3..fcca64f 100644
--- a/test/quotient.jl
+++ b/test/quotient.jl
@@ -1,19 +1,27 @@
-struct PlusMinusOne # Reduce modulo `x^2 = 1`
+struct ImRem # Reduce modulo `x^2 = -1`
 end
-function Base.rem(m::AbstractMonomial, ::PlusMinusOne)
-    return prod(v^mod(e, 2) for (v, e) in powers(m))
+function Base.rem(m::AbstractMonomial, ::ImRem)
+    return prod((-1)^rem(e, 2) * v^mod(e, 2) for (v, e) in powers(m))
 end
-function Base.rem(t::AbstractTerm, d::PlusMinusOne)
+function Base.rem(t::AbstractTerm, d::ImRem)
     return coefficient(t) * rem(monomial(t), d)
 end
-function Base.rem(p::AbstractPolynomial, d::PlusMinusOne)
+function Base.rem(p::AbstractPolynomial, d::ImRem)
     return sum(rem(t, d) for t in terms(p))
 end
 
-@testset "PlusMinusOne" begin
+@testset "ImRem" begin
     @polyvar x y
-    basis =
-        MB.QuotientBasis(MB.SubBasis{MB.Monomial}([1, y, x]), PlusMinusOne())
+    basis = MB.QuotientBasis(MB.SubBasis{MB.Monomial}([1, y, x]), ImRem())
     @test length(basis) == 3
-    @test coefficients(x^3 - 2x^2 * y + 3x^2, basis) == [3, -2, 1]
+    p = x^3 - 2x^2 * y + 3x^2
+    coeffs = [3, 2, -1]
+    @test coefficients(p, basis) == coeffs
+    a = MB.algebra_element(p)
+    @test SA.coeffs(a, basis) == coeffs
+    explicit = MB.explicit_basis(a)
+    exp = MB.algebra_element(SA.coeffs(a, explicit), explicit)
+    @test SA.coeffs(exp, basis) == coeffs
+    q = MB.algebra_element(coeffs, basis)
+    @test SA.adjoint_coeffs(q, explicit) == [3, -2, 1]
 end
diff --git a/test/runtests.jl b/test/runtests.jl
index 984dce0..ec8fc2f 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,5 +1,7 @@
 using Test
 
+import MutableArithmetics as MA
+import StarAlgebras as SA
 using MultivariateBases
 const MB = MultivariateBases
 using LinearAlgebra
@@ -9,11 +11,30 @@ function api_test(B::Type{<:MB.AbstractMonomialIndexed}, degree)
     @polyvar x[1:2]
     M = typeof(prod(x))
     full_basis = FullBasis{B,M}()
+    @test sprint(show, MB.algebra(full_basis)) ==
+          "Polynomial algebra of $B basis"
+    @test typeof(MB.algebra(full_basis)) ==
+          MA.promote_operation(MB.algebra, typeof(full_basis))
     for basis in [
         maxdegree_basis(full_basis, x, degree),
-        basis_covering_monomials(full_basis, monomials(x, 0:degree)),
+        explicit_basis_covering(
+            full_basis,
+            MB.SubBasis{MB.Monomial}(monomials(x, 0:degree)),
+        ),
+        explicit_basis_covering(
+            full_basis,
+            MB.SubBasis{ScaledMonomial}(monomials(x, 0:degree)),
+        ),
     ]
+        @test typeof(MB.algebra(basis)) ==
+              MA.promote_operation(MB.algebra, typeof(basis))
         @test basis isa MB.explicit_basis_type(typeof(full_basis))
+        for i in eachindex(basis)
+            mono = basis.monomials[i]
+            poly = MB.Polynomial{B}(mono)
+            @test basis[i] == poly
+            @test basis[poly] == i
+        end
         n = binomial(2 + degree, 2)
         @test length(basis) == n
         @test firstindex(basis) == 1
@@ -25,12 +46,32 @@ function api_test(B::Type{<:MB.AbstractMonomialIndexed}, degree)
         @test typeof(empty_basis(typeof(basis))) == typeof(basis)
         @test length(empty_basis(typeof(basis))) == 0
         @test polynomial_type(basis, Float64) == polynomial_type(x[1], Float64)
-        @test polynomial(i -> 0.0, basis) isa polynomial_type(basis, Float64)
-        @test polynomial(zeros(n, n), basis, Float64) isa
-              polynomial_type(basis, Float64)
-        @test polynomial(ones(n, n), basis, Float64) isa
-              polynomial_type(basis, Float64)
+        #@test polynomial(i -> 0.0, basis) isa polynomial_type(basis, Float64)
     end
+    mono = x[1]^2 * x[2]^3
+    p = MB.Polynomial{B}(mono)
+    @test full_basis[p] == mono
+    @test full_basis[mono] == p
+    @test polynomial_type(mono, String) == polynomial_type(typeof(p), String)
+    a = MB.algebra_element(p)
+    @test typeof(polynomial(a)) == polynomial_type(typeof(a))
+    @test typeof(polynomial(a)) == polynomial_type(typeof(p), Int)
+    @test a ≈ a
+    if B == MB.Monomial
+        @test a ≈ p.monomial
+        @test p.monomial ≈ a
+    else
+        @test !(a ≈ p.monomial)
+        @test !(p.monomial ≈ a)
+    end
+    _wrap(s) = (B == MB.Monomial ? s : "$B($s)")
+    @test sprint(show, p) == _wrap(sprint(show, p.monomial))
+    @test sprint(print, p) == _wrap(sprint(print, p.monomial))
+    mime = MIME"text/latex"()
+    @test sprint(show, mime, p) ==
+          "\$\$ " *
+          _wrap(MB.SA.trim_LaTeX(mime, sprint(show, mime, p.monomial))) *
+          " \$\$"
 end
 
 function univ_orthogonal_test(
@@ -67,10 +108,10 @@ function orthogonal_test(
         end
     end
 
-    @testset "basis_covering_monomials" begin
-        basis = basis_covering_monomials(
+    @testset "explicit_basis_covering" begin
+        basis = explicit_basis_covering(
             FullBasis{B,typeof(x * y)}(),
-            monomial_vector([x^2 * y, y^2]),
+            SubBasis{MB.Monomial}(monomial_vector([x^2 * y, y^2])),
         )
         if even_odd_separated
             exps = [(0, 0), (0, 1), (0, 2), (2, 1)]
@@ -81,9 +122,9 @@ function orthogonal_test(
             @test polynomial(basis[i]) ==
                   univariate_x[exps[i][1]+1] * univariate_y[exps[i][2]+1]
         end
-        basis = basis_covering_monomials(
+        basis = explicit_basis_covering(
             FullBasis{B,typeof(x^2)}(),
-            monomial_vector([x^4, x^2, x]),
+            SubBasis{MB.Monomial}(monomial_vector([x^4, x^2, x])),
         )
         if even_odd_separated
             exps = [0, 1, 2, 4]
@@ -99,7 +140,7 @@ end
 function coefficient_test(basis::SubBasis, p, coefs; kwargs...)
     cc = coefficients(p, basis)
     @test isapprox(coefs, cc; kwargs...)
-    @test isapprox(p, polynomial(cc, basis); kwargs...)
+    @test isapprox(p, algebra_element(cc, basis); kwargs...)
 end
 
 function coefficient_test(basis::SubBasis, coefs::AbstractVector; kwargs...)
@@ -118,7 +159,10 @@ function coefficient_test(
 )
     @polyvar x y
     p = x^4 * y^2 + x^2 * y^4 - 3 * x^2 * y^2 + 1
-    basis = basis_covering_monomials(FullBasis{B,typeof(x * y)}(), monomials(p))
+    basis = explicit_basis_covering(
+        FullBasis{B,typeof(x * y)}(),
+        SubBasis{MB.Monomial}(monomials(p)),
+    )
     coefficient_test(basis, p, coefs; kwargs...)
     return
 end