diff --git a/.github/workflows/format_check.yml b/.github/workflows/format_check.yml
index a90c551..e62b574 100644
--- a/.github/workflows/format_check.yml
+++ b/.github/workflows/format_check.yml
@@ -19,7 +19,7 @@ jobs:
         run: |
           using Pkg
           # If you update the version, also update the style guide docs.
-          Pkg.add(PackageSpec(name="JuliaFormatter", version="1"))
+          Pkg.add(PackageSpec(name="JuliaFormatter", version="2"))
           using JuliaFormatter
           format("src", verbose=true)
           format("test", verbose=true)
diff --git a/Project.toml b/Project.toml
index 8b0bf03..39423e4 100644
--- a/Project.toml
+++ b/Project.toml
@@ -5,11 +5,13 @@ version = "0.2.2"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+MergeSorted = "5f6ebb72-d4a9-4954-9a69-49758639d560"
 MultivariatePolynomials = "102ac46a-7ee4-5c85-9060-abc95bfdeaa3"
 MutableArithmetics = "d8a4904e-b15c-11e9-3269-09a3773c0cb0"
 StarAlgebras = "0c0c59c1-dc5f-42e9-9a8b-b5dc384a6cd1"
 
 [compat]
+MergeSorted = "2.0.2"
 MultivariatePolynomials = "0.5.3"
 MutableArithmetics = "1"
 julia = "1.6"
diff --git a/src/MultivariateBases.jl b/src/MultivariateBases.jl
index 1ff3ba6..2364451 100644
--- a/src/MultivariateBases.jl
+++ b/src/MultivariateBases.jl
@@ -8,39 +8,50 @@ export FullBasis, SubBasis
 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)
+struct Variables{B,V}
+    variables::V
 end
-struct Algebra{BT,B,M} <:
-       SA.AbstractStarAlgebra{Polynomial{B,M},Polynomial{B,M}}
-    basis::BT
+
+Variables{B}(vars) where {B} = Variables{B,typeof(vars)}(vars)
+
+function variable_index(v::Variables, var)
+    return findfirst(isequal(var), v.variables)
 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)
+
+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
-function MA.promote_operation(
-    ::typeof(SA.basis),
-    ::Type{<:Algebra{B}},
-) where {B}
-    return B
+
+function Base.:(==)(v::Variables{B}, w::Variables{B}) where {B}
+    # Testing `===` allows speeding up a typical situations
+    return v.variables === w.variables || v.variables == w.variables
 end
-SA.basis(a::Algebra) = a.basis
 
-#Base.:(==)(::Algebra{BT1,B1,M}, ::Algebra{BT2,B2,M}) where {BT1,B1,BT2,B2,M} = true
-function Base.:(==)(a::Algebra, b::Algebra)
-    # `===` is a shortcut for speedup
-    return a.basis === b.basis || a.basis == b.basis
+MP.monomial_type(::Type{Variables{B,V}}) where {B,V} = MP.monomial_type(V)
+
+constant_monomial_exponents(v::Variables) = map(_ -> 0, v.variables)
+
+function (v::Variables)(exponents)
+    return Polynomial(v, exponents)
 end
 
-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")
+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
+MP.monomial_type(::Type{<:SA.StarAlgebra{O}}) where {O} = MP.monomial_type(O)
+function MP.polynomial_type(::Type{A}, ::Type{T}) where {A<:SA.StarAlgebra,T}
+    return MP.polynomial_type(MA.promote_operation(SA.basis, A), T)
 end
 
+include("bases.jl")
+include("mstructures.jl")
 include("monomial.jl")
 include("scaled.jl")
 
@@ -67,18 +78,26 @@ include("lagrange.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)
+    basis::Union{QuotientBasis{<:Polynomial{B}},FullBasis{B},SubBasis{B}},
+) where {B}
+    return SA.StarAlgebra(Variables{B}(MP.variables(basis)), MStruct(basis))
 end
 
 function MA.promote_operation(
     ::typeof(algebra),
     BT::Type{
-        <:Union{QuotientBasis{Polynomial{B,M}},FullBasis{B,M},SubBasis{B,M}},
+        <:Union{
+            QuotientBasis{Polynomial{B,V,E}},
+            FullBasis{B,V,E},
+            SubBasis{B,V,E},
+        },
     },
-) where {B,M}
-    return Algebra{BT,B,M}
+) where {B,V,E}
+    return SA.StarAlgebra{
+        Variables{B,V},
+        Polynomial{B,V,E},
+        MStruct{B,V,E,SA.key_type(BT),BT},
+    }
 end
 
 include("arithmetic.jl")
diff --git a/src/arithmetic.jl b/src/arithmetic.jl
index d94e626..a8da053 100644
--- a/src/arithmetic.jl
+++ b/src/arithmetic.jl
@@ -1,6 +1,21 @@
 const _APL = MP.AbstractPolynomialLike
 # We don't define it for all `AlgebraElement` as this would be type piracy
-const _AE = SA.AlgebraElement{<:Algebra}
+const _AE = SA.AlgebraElement{<:SA.StarAlgebra{<:Variables}}
+
+_collect_if_tuple(t::Tuple) = collect(t)
+_collect_if_tuple(v::AbstractVector) = v
+
+function _polynomial(b::FullBasis{Monomial}, c::SA.SparseCoefficients)
+    return MP.polynomial(
+        collect(SA.values(c)),
+        _collect_if_tuple(MP.monomial.(getindex.(Ref(b), SA.keys(c)))),
+    )
+end
+
+function MP.polynomial(a::SA.AlgebraElement)
+    b = FullBasis{Monomial}(MP.variables(a))
+    return _polynomial(b, SA.coeffs(a, b))
+end
 
 for op in [:+, :-, :*]
     @eval begin
@@ -52,7 +67,7 @@ for op in [:+, :-]
         end
         function Base.$op(p, q::_AE)
             i = implicit(q)
-            return $op(constant_algebra_element(typeof(SA.basis(i)), p), i)
+            return $op(constant_algebra_element(SA.basis(i), p), i)
         end
         function MA.promote_operation(
             ::typeof($op),
@@ -71,15 +86,15 @@ for op in [:+, :-]
         end
         function Base.$op(p::_AE, q)
             i = implicit(p)
-            return $op(i, constant_algebra_element(typeof(SA.basis(i)), q))
+            return $op(i, constant_algebra_element(SA.basis(i), q))
         end
     end
 end
 
-function term_element(α, p::Polynomial{B,M}) where {B,M}
+function term_element(α, p::Polynomial{B}) where {B}
     return algebra_element(
-        sparse_coefficients(MP.term(α, p.monomial)),
-        FullBasis{B,M}(),
+        SA.SparseCoefficients((p.exponents,), (α,)),
+        FullBasis{B}(MP.variables(p)),
     )
 end
 
diff --git a/src/bases.jl b/src/bases.jl
new file mode 100644
index 0000000..aacafa2
--- /dev/null
+++ b/src/bases.jl
@@ -0,0 +1,283 @@
+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 MP.ordering(::Type{MonomialIndexedBasis{B,V,E}}) where {B,V,E}
+    return MP.ordering(E)
+end
+MP.ordering(b::MonomialIndexedBasis) = MP.ordering(typeof(b))
+
+function FullBasis{B}(vars) where {B}
+    O = MP.ordering(vars)
+    v = Variables{B,typeof(vars)}(vars)
+    exps = MP.ExponentsIterator{O}(constant_monomial_exponents(v))
+    return SA.MappedBasis{Polynomial{B,typeof(vars),eltype(exps)}}(
+        exps,
+        v,
+        MP.exponents,
+    )
+end
+
+function FullBasis{B}(p::MP.AbstractPolynomialLike) where {B}
+    return FullBasis{B}(MP.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)
+
+MP.monomial_type(::Type{<:FullBasis{B,V}}) where {B,V} = MP.monomial_type(V)
+function MP.polynomial_type(basis::FullBasis, ::Type{T}) where {T}
+    return MP.polynomial_type(typeof(basis), T)
+end
+
+MP.nvariables(v::Variables) = length(v.variables)
+MP.nvariables(basis::FullBasis) = MP.nvariables(basis.map)
+MP.nvariables(basis::SubBasis) = MP.nvariables(parent(basis))
+MP.variables(v::Variables) = v.variables
+MP.variables(basis::FullBasis) = MP.variables(basis.map)
+MP.variables(basis::SubBasis) = MP.variables(parent(basis))
+
+function explicit_basis_covering(
+    full::FullBasis{B},
+    target::SubBasis{B},
+) where {B}
+    return SA.SubBasis(full, target.keys)
+end
+
+function MP.monomial_type(::Type{<:SA.SubBasis{T,I,K,B}}) where {T,I,K,B}
+    return MP.monomial_type(B)
+end
+
+# The `i`th index of output is the index of occurence of `x[i]` in `y`,
+# or `0` if it does not occur.
+function multi_findsorted(x, y)
+    I = zeros(Int, length(x))
+    j = 1
+    for i in eachindex(x)
+        while j ≤ length(y) && x[i] > y[j]
+            j += 1
+        end
+        if j ≤ length(y) && x[i] == y[j]
+            I[i] = j
+        end
+    end
+    return I
+end
+
+import MergeSorted
+
+function merge_bases(basis1::MB, basis2::MB) where {MB<:SubBasis}
+    @assert basis1.parent_basis == basis2.parent_basis
+    @assert basis1.is_sorted
+    @assert basis2.is_sorted
+    keys = unique!(
+        MergeSorted.mergesorted(
+            basis1.keys,
+            basis2.keys;
+            lt = SA.comparable(basis1),
+        ),
+    )
+    I1 = multi_findsorted(keys, basis1.keys)
+    I2 = multi_findsorted(keys, basis2.keys)
+    return SA.SubBasis(basis1.parent_basis, keys), I1, I2
+end
+
+# Unsafe because we don't check that `monomials` is sorted and without duplicates
+function unsafe_basis(
+    ::Type{B},
+    monomials::AbstractVector{<:MP.AbstractMonomial},
+) where {B<:AbstractMonomialIndexed}
+    # TODO We should add a way to directly get the vector of exponents inside `DynamicPolynomials.MonomialVector`.
+    # `MP.exponents.(monomials)` should work in the meantime even if it's not the most efficient
+    return SA.SubBasis(
+        FullBasis{B}(MP.variables(monomials)),
+        MP.exponents.(monomials),
+    )
+end
+
+function Base.getindex(
+    ::FullBasis{B},
+    monomials::AbstractVector{M},
+) where {B,M<:MP.AbstractMonomial}
+    return unsafe_basis(B, MP.monomial_vector(monomials)::AbstractVector{M})
+end
+
+function SubBasis{B}(
+    monomials::AbstractVector{<:MP.AbstractTermLike},
+) where {B<:AbstractMonomialIndexed}
+    return unsafe_basis(
+        B,
+        MP.monomial_vector(monomials)::AbstractVector{<:MP.AbstractMonomial},
+    )
+end
+
+SubBasis{B}(monos::Tuple) where {B} = SubBasis{B}([monos...])
+
+function algebra_type(::Type{BT}) where {B,M,BT<:MonomialIndexedBasis{B,M}}
+    return Algebra{BT,B,M}
+end
+
+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 _algebra_element_type(
+    ::Type{C}, # type of coefficient vector
+    ::Type{B}, # basis type
+) where {C,B}
+    A = MA.promote_operation(algebra, B)
+    T = eltype(C) # Even works for `NTuple`!
+    E = SA.key_type(B)
+    return SA.AlgebraElement{
+        A,
+        T,
+        SA.SparseCoefficients{E,T,_similar_type(C, E),C,typeof(isless)},
+    }
+end
+
+function MA.promote_operation(
+    ::typeof(implicit),
+    ::Type{AE},
+) where {AG,T,AE<:SA.AlgebraElement{AG,T}}
+    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::MP.AbstractPolynomialLike)
+    return SubBasis{Monomial}(MP.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, MP.empty_monomial_vector(M))
+end
+
+function maxdegree_basis(
+    ::FullBasis{B},
+    variables,
+    maxdegree::Int,
+) where {B<:AbstractMonomialIndexed}
+    return unsafe_basis(B, MP.monomials(variables, 0:maxdegree))
+end
+
+MP.variables(c::SA.AbstractCoefficients) = MP.variables(SA.keys(c))
+
+_lazy_collect(v::AbstractVector) = collect(v)
+_lazy_collect(v::Tuple) = collect(v)
+_lazy_collect(v::Vector) = v
+
+function sparse_coefficients(p::MP.AbstractPolynomial)
+    return SA.SparseCoefficients(
+        MP.exponents.(MP.monomials(p)),
+        _lazy_collect(MP.coefficients(p)),
+    )
+end
+
+function sparse_coefficients(t::MP.AbstractTermLike)
+    return SA.SparseCoefficients((MP.exponents(t),), (MP.coefficient(t),))
+end
+
+function algebra_element(p::MP.AbstractPolynomialLike)
+    return algebra_element(sparse_coefficients(p), FullBasis{Monomial}(p))
+end
+
+function algebra_element(f::Function, basis::SubBasis)
+    return algebra_element(map(f, eachindex(basis)), basis)
+end
+
+# Another option is to use `fieldcount`, see
+# https://discourse.julialang.org/t/get-tuple-length-from-type/32483/16?u=blegat
+_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 MA.promote_operation(
+    ::typeof(algebra_element),
+    P::Type{<:MP.AbstractPolynomialLike{T}},
+) where {T}
+    V = MA.promote_operation(MP.variables, P)
+    E = _similar_type(V, Int)
+    O = MP.ordering(P)
+    P = MultivariateBases.Polynomial{Monomial,V,E}
+    B = SA.MappedBasis{
+        P,
+        E,
+        MP.ExponentsIterator{O,Nothing,E},
+        Variables{Monomial,V},
+        typeof(MP.exponents),
+    }
+    return _algebra_element_type(Vector{T}, B)
+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 algebra_element(
+        SA.SparseCoefficients(
+            (constant_monomial_exponents(b),),
+            (_one_if_type(α),),
+        ),
+        b,
+    )
+end
+
+function constant_algebra_element_type(
+    ::Type{B},
+    ::Type{T},
+) where {B<:SubBasis,T}
+    A = MA.promote_operation(algebra, B)
+    return SA.AlgebraElement{A,T,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
diff --git a/src/chebyshev.jl b/src/chebyshev.jl
index aeb35d4..7e3dff3 100644
--- a/src/chebyshev.jl
+++ b/src/chebyshev.jl
@@ -24,57 +24,53 @@ const Chebyshev = ChebyshevFirstKind
 
 # https://en.wikipedia.org/wiki/Chebyshev_polynomials#Properties
 # `T_n * T_m = T_{n + m} / 2 + T_{|n - m|} / 2`
-function univariate_mul!(::Mul{Chebyshev}, terms, var, a, b)
-    I = eachindex(terms)
-    for i in I
-        mono = MP.monomial(terms[i]) * var^(a + b)
-        terms[i] = MA.mul!!(terms[i], var^abs(a - b))
-        terms[i] = MA.operate!!(/, terms[i], 2)
-        α = MA.copy_if_mutable(MP.coefficient(terms[i]))
-        push!(terms, MP.term(α, mono))
+function univariate_mul!(::Type{Chebyshev}, exps, coefs, var, a, b)
+    for i in eachindex(exps)
+        exp = _increment(exps[i], a + b, var)
+        exps[i] = _increment!(exps[i], abs(a - b), var)
+        coefs[i] = MA.operate!!(/, coefs[i], 2)
+        push!(coefs, MA.copy_if_mutable(coefs[i]))
+        push!(exps, exp)
     end
     return
 end
 
-function (mul::Mul{B})(
-    a::MP.AbstractMonomial,
-    b::MP.AbstractMonomial,
-) where {B<:AbstractMonomialIndexed}
-    terms = [MP.term(1 // 1, MP.constant_monomial(a * b))]
-    vars_a = MP.effective_variables(a)
+function (m::MStruct{B,V,E})(
+    a::E,
+    b::E,
+    ::Type{E},
+) where {B<:AbstractMonomialIndexed,V,E}
+    exps = [zero.(a)]
+    coefs = [1 // 1]
+    vars_a = findall(!iszero, a)
     var_state_a = iterate(vars_a)
-    vars_b = MP.effective_variables(b)
+    vars_b = findall(!iszero, b)
     var_state_b = iterate(vars_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))
+            for i in eachindex(exps)
+                exps[i] = _increment!(exps[i], b[var_b], var_b)
             end
             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))
+            for i in eachindex(exps)
+                exps[i] = _increment!(exps[i], a[var_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
-            univariate_mul!(
-                mul,
-                terms,
-                var_a,
-                MP.degree(a, var_a),
-                MP.degree(b, var_b),
-            )
+            univariate_mul!(B, exps, coefs, var_a, a[var_a], b[var_b])
             var_state_a = iterate(vars_a, state_a)
             var_state_b = iterate(vars_b, state_b)
         end
     end
-    return sparse_coefficients(MP.polynomial!(terms))
+    # FIXME is `canonical` needed ?
+    return MA.operate!(SA.canonical, SA.SparseCoefficients(exps, coefs))
 end
 
 function _add_mul_scalar_vector!(res, ::SubBasis, scalar, vector)
@@ -119,10 +115,10 @@ function SA.coeffs(
 )
     sub = explicit_basis_covering(target, source)
     # Need to make A square so that it's UpperTriangular
-    extended = SubBasis{Monomial}(sub.monomials)
-    ext = SA.coeffs(algebra_element(cfs, source), extended)
+    extended = SA.SubBasis(parent(source), sub.keys)
+    ext = _collect_if_tuple(SA.coeffs(algebra_element(cfs, source), extended))
     return SA.SparseCoefficients(
-        sub.monomials,
+        sub.keys,
         #transformation_to(sub, extended) \ ext, # Julia v1.6 converts the matrix to the eltype of the `result` which is bad for JuMP
         LinearAlgebra.ldiv!(
             zeros(_promote_coef(eltype(ext), Chebyshev), length(sub)),
@@ -132,10 +128,10 @@ function SA.coeffs(
     )
 end
 
-function SA.coeffs(cfs, ::FullBasis{Monomial}, target::FullBasis{Chebyshev})
+function SA.coeffs(cfs, src::FullBasis{Monomial}, target::FullBasis{Chebyshev})
     return SA.coeffs(
         SA.values(cfs),
-        SubBasis{Monomial}(collect(SA.keys(cfs))),
+        SA.SubBasis(src, collect(SA.keys(cfs))),
         target,
     )
 end
@@ -151,7 +147,7 @@ function _scalar_product_function(::Type{Chebyshev}, i::Int)
         return 0
     else
         n = div(i, 2)
-        return (π / 2^i) * prod(n+1:i) / factorial(n)
+        return (π / 2^i) * prod((n+1):i) / factorial(n)
     end
 end
 
@@ -173,6 +169,6 @@ function _scalar_product_function(::Type{<:ChebyshevSecondKind}, i::Int)
         return 0
     else
         n = div(i, 2)
-        return π / (2^(i + 1)) * prod(n+2:i) / factorial(n)
+        return π / (2^(i + 1)) * prod((n+2):i) / factorial(n)
     end
 end
diff --git a/src/fixed.jl b/src/fixed.jl
index ddba9f3..63fea52 100644
--- a/src/fixed.jl
+++ b/src/fixed.jl
@@ -1,26 +1,3 @@
-"""
-    struct FixedBasis{B,M,T,V} <:
-        SA.ExplicitBasis{SA.AlgebraElement{Algebra{FullBasis{B,M},B,M},T,V},Int}
-        elements::Vector{SA.AlgebraElement{Algebra{FullBasis{B,M},B,M},T,V}}
-    end
-
-Fixed basis with polynomials `elements`.
-"""
-struct FixedBasis{B,M,T,V} <:
-       SA.ExplicitBasis{SA.AlgebraElement{Algebra{FullBasis{B,M},B,M},T,V},Int}
-    elements::Vector{SA.AlgebraElement{Algebra{FullBasis{B,M},B,M},T,V}}
-end
-
-Base.length(b::FixedBasis) = length(b.elements)
-Base.getindex(b::FixedBasis, i::Integer) = b.elements[i]
-
-function Base.show(io::IO, b::FixedBasis)
-    print(io, "FixedBasis(")
-    _show_vector(io, MIME"text/plain"(), b.elements)
-    print(io, ")")
-    return
-end
-
 # TODO refactor with `SA.MappedBasis`
 # https://github.com/JuliaAlgebra/StarAlgebras.jl/pull/76
 
diff --git a/src/hermite.jl b/src/hermite.jl
index b98a041..2611d3c 100644
--- a/src/hermite.jl
+++ b/src/hermite.jl
@@ -30,7 +30,7 @@ function _scalar_product_function(::Type{ProbabilistsHermite}, i::Int)
         return 0
     else
         n = div(i, 2)
-        return (√(2 * π) / (2^n)) * prod(n+1:2*n)
+        return (√(2 * π) / (2^n)) * prod((n+1):(2*n))
     end
 end
 
@@ -55,6 +55,6 @@ function _scalar_product_function(::Type{PhysicistsHermite}, i::Int)
         return 0
     else
         n = div(i, 2)
-        return (√(π) / (2^i)) * prod(n+1:i)
+        return (√(π) / (2^i)) * prod((n+1):i)
     end
 end
diff --git a/src/lagrange.jl b/src/lagrange.jl
index 7e4fcae..b756928 100644
--- a/src/lagrange.jl
+++ b/src/lagrange.jl
@@ -103,13 +103,13 @@ function eval_basis!(
         end
     end
     for i in eachindex(values)
-        l = MP.maxdegree(basis.monomials, variables[i]) + 1
+        l = maximum(Base.Fix2(getindex, i), basis.keys) + 1
         univariate_eval!(B, view(univariate_buffer, 1:l, i), values[i])
     end
     for i in eachindex(basis)
         result[i] = one(eltype(result))
         for j in eachindex(values)
-            d = MP.degree(basis.monomials[i], variables[j])
+            d = basis.keys[i][j]
             result[i] = MA.operate!!(*, result[i], univariate_buffer[d+1, j])
         end
     end
diff --git a/src/monomial.jl b/src/monomial.jl
index a6d3d23..d88deba 100644
--- a/src/monomial.jl
+++ b/src/monomial.jl
@@ -1,320 +1,24 @@
-struct FullBasis{B<:AbstractMonomialIndexed,M<:MP.AbstractMonomial} <:
-       SA.ImplicitBasis{Polynomial{B,M},M} end
-
-function Base.getindex(::FullBasis{B,M}, mono::M) where {B,M}
-    return Polynomial{B}(mono)
-end
-
-function Base.getindex(::FullBasis{B,M}, p::Polynomial{B,M}) where {B,M}
-    return p.monomial
-end
-
-SA.mstructure(::FullBasis{B}) where {B} = Mul{B}()
-
-MP.monomial_type(::Type{<:FullBasis{B,M}}) where {B,M} = M
-function MP.polynomial_type(basis::FullBasis{B,M}, ::Type{T}) where {B,M,T}
-    return MP.polynomial_type(typeof(basis), T)
-end
-
-# TODO Move it to SA as :
-# struct SubBasis{T,I,B<:SA.ImplicitBasis{T,I},V<:AbstractVector{I}} <: SA.ExplicitBasis{T,Int}
-#     implicit::B
-#     indices::V
-# end
-struct SubBasis{B<:AbstractMonomialIndexed,M,V<:AbstractVector{M}} <:
-       SA.ExplicitBasis{Polynomial{B,M},Int}
-    monomials::V
-end
-
-# Overload some of the `AbstractVector` interface for convenience
-Base.isempty(basis::SubBasis) = isempty(basis.monomials)
-Base.eachindex(basis::SubBasis) = eachindex(basis.monomials)
-_iterate(::SubBasis, ::Nothing) = nothing
-function _iterate(basis::SubBasis{B}, elem_state) where {B}
-    return parent(basis)[elem_state[1]], elem_state[2]
-end
-Base.iterate(basis::SubBasis) = _iterate(basis, iterate(basis.monomials))
-Base.iterate(basis::SubBasis, s) = _iterate(basis, iterate(basis.monomials, s))
-Base.length(basis::SubBasis) = length(basis.monomials)
-Base.firstindex(basis::SubBasis) = firstindex(basis.monomials)
-Base.lastindex(basis::SubBasis) = lastindex(basis.monomials)
-
-MP.nvariables(basis::SubBasis) = MP.nvariables(basis.monomials)
-MP.variables(basis::SubBasis) = MP.variables(basis.monomials)
-
-Base.parent(::SubBasis{B,M}) where {B,M} = FullBasis{B,M}()
-
-function Base.getindex(basis::SubBasis, index::Int)
-    return parent(basis)[basis.monomials[index]]
-end
-
-function monomial_index(basis::SubBasis, mono::MP.AbstractMonomial)
-    i = searchsortedfirst(basis.monomials, mono)
-    if i in eachindex(basis.monomials) && basis.monomials[i] == mono
-        return i
-    end
-    return
-end
-
-function Base.getindex(basis::SubBasis{B,M}, value::Polynomial{B,M}) where {B,M}
-    mono = monomial_index(basis, parent(basis)[value])
-    if isnothing(mono)
-        throw(BoundsError(basis, value))
-    end
-    return mono
-end
-
-function explicit_basis_covering(::FullBasis{B}, target::SubBasis{B}) where {B}
-    return SubBasis{B}(target.monomials)
-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`,
-# or `0` if it does not occur.
-function multi_findsorted(x, y)
-    I = zeros(Int, length(x))
-    j = 1
-    for i in eachindex(x)
-        while j ≤ length(y) && x[i] > y[j]
-            j += 1
-        end
-        if j ≤ length(y) && x[i] == y[j]
-            I[i] = j
-        end
-    end
-    return I
-end
-
-function merge_bases(basis1::MB, basis2::MB) where {MB<:SubBasis}
-    monos = MP.merge_monomial_vectors([basis1.monomials, basis2.monomials])
-    I1 = multi_findsorted(monos, basis1.monomials)
-    I2 = multi_findsorted(monos, basis2.monomials)
-    return MB(monos), I1, I2
-end
-
-# Unsafe because we don't check that `monomials` is sorted and without duplicates
-function unsafe_basis(
-    ::Type{B},
-    monomials::AbstractVector{M},
-) where {B<:AbstractMonomialIndexed,M<:MP.AbstractMonomial}
-    return SubBasis{B,M,typeof(monomials)}(monomials)
-end
-
-function Base.getindex(::FullBasis{B,M}, monomials::AbstractVector) where {B,M}
-    return unsafe_basis(B, MP.monomial_vector(monomials)::AbstractVector{M})
-end
-
-function SubBasis{B}(
-    monomials::AbstractVector,
-) where {B<:AbstractMonomialIndexed}
-    return unsafe_basis(
-        B,
-        MP.monomial_vector(monomials)::AbstractVector{<:MP.AbstractMonomial},
-    )
-end
-
-SubBasis{B}(monos::Tuple) where {B} = SubBasis{B}([monos...])
-
-function Base.copy(basis::SubBasis)
-    return typeof(basis)(copy(basis.monomials))
-end
-
-function Base.:(==)(a::SubBasis{B}, b::SubBasis{B}) where {B}
-    return a.monomials == b.monomials
-end
-
-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 implicit(a::SA.AlgebraElement)
-    basis = implicit_basis(SA.basis(a))
-    return algebra_element(SA.coeffs(a, basis), basis)
-end
-
-function MA.promote_operation(
-    ::typeof(implicit),
-    ::Type{E},
-) where {AG,T,E<:SA.AlgebraElement{AG,T}}
-    BT = MA.promote_operation(implicit_basis, MA.promote_operation(SA.basis, E))
-    A = MA.promote_operation(algebra, BT)
-    M = MP.monomial_type(BT)
-    return SA.AlgebraElement{A,T,SA.SparseCoefficients{M,T,Vector{M},Vector{T}}}
-end
-
-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(p::MP.AbstractPolynomialLike)
-    return SubBasis{Monomial}(MP.monomials(p))
-end
-
-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}
-    return SubBasis{B,M,MP.monomial_vector_type(M)}
-end
-
-function empty_basis(
-    ::Type{<:SubBasis{B,M}},
-) where {B<:AbstractMonomialIndexed,M}
-    return unsafe_basis(B, MP.empty_monomial_vector(M))
-end
-
-function maxdegree_basis(
-    ::FullBasis{B},
-    variables,
-    maxdegree::Int,
-) where {B<:AbstractMonomialIndexed}
-    return unsafe_basis(B, MP.monomials(variables, 0:maxdegree))
-end
-
-MP.variables(c::SA.AbstractCoefficients) = MP.variables(SA.keys(c))
-
-_lazy_collect(v::AbstractVector) = collect(v)
-_lazy_collect(v::Vector) = v
-
-function sparse_coefficients(p::MP.AbstractPolynomial)
-    return SA.SparseCoefficients(
-        _lazy_collect(MP.monomials(p)),
-        _lazy_collect(MP.coefficients(p)),
-    )
-end
-
-function sparse_coefficients(t::MP.AbstractTermLike)
-    return SA.SparseCoefficients((MP.monomial(t),), (MP.coefficient(t),))
-end
-
-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,Vector{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
-
-_one_if_type(α) = α
-_one_if_type(::Type{T}) where {T} = one(T)
-
-function constant_algebra_element_type(
-    ::Type{BT},
-    ::Type{T},
-) where {B,M,BT<:FullBasis{B,M},T}
-    A = MA.promote_operation(algebra, BT)
-    return SA.AlgebraElement{A,T,SA.SparseCoefficients{M,T,Vector{M},Vector{T}}}
-end
-
-function constant_algebra_element(::Type{FullBasis{B,M}}, α) where {B,M}
-    return algebra_element(
-        sparse_coefficients(
-            MP.polynomial(MP.term(_one_if_type(α), MP.constant_monomial(M))),
-        ),
-        FullBasis{B,M}(),
-    )
-end
-
-function constant_algebra_element_type(
-    ::Type{B},
-    ::Type{T},
-) where {B<:SubBasis,T}
-    A = MA.promote_operation(algebra, B)
-    return SA.AlgebraElement{A,T,Vector{T}}
-end
-
-function constant_algebra_element(::Type{<:SubBasis{B,M}}, α) where {B,M}
-    return algebra_element(
-        [_one_if_type(α)],
-        SubBasis{B}([MP.constant_monomial(M)]),
-    )
-end
-
-# TODO use Base.show_vector here, maybe by wrapping the `generator` vector
-#      into something that spits objects wrapped with the `mime` type
-function _show_vector(io::IO, mime::MIME, v)
-    print(io, '[')
-    first = true
-    for el in v
-        if !first
-            print(io, ", ")
-        end
-        first = false
-        show(io, mime, 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.monomials)
-    print(io, ')')
-    return
-end
-
-function Base.show(io::IO, mime::MIME"text/plain", basis::SubBasis)
-    return _show(io, mime, basis)
-end
-
-function Base.show(io::IO, mime::MIME"text/print", basis::SubBasis)
-    return _show(io, mime, basis)
-end
-
-function Base.print(io::IO, basis::SubBasis)
-    return show(io, MIME"text/print"(), basis)
-end
-
-function Base.show(io::IO, basis::SubBasis)
-    return show(io, MIME"text/plain"(), basis)
-end
-
 abstract type AbstractMonomial <: AbstractMonomialIndexed end
 
 function explicit_basis_covering(
-    ::FullBasis{B},
+    full::FullBasis{B},
     target::SubBasis{<:AbstractMonomial},
 ) where {B<:AbstractMonomial}
-    return SubBasis{B}(target.monomials)
+    return SA.SubBasis(full, target.keys)
 end
 
 # To break ambiguity
 function explicit_basis_covering(
-    ::FullBasis{B},
+    full::FullBasis{B},
     target::SubBasis{B},
 ) where {B<:AbstractMonomial}
-    return SubBasis{B}(target.monomials)
+    @assert full == parent(target)
+    return SA.SubBasis(parent(target), target.keys)
 end
 
 function Base.adjoint(p::Polynomial{B}) where {B<:AbstractMonomialIndexed}
-    return Polynomial{B}(adjoint(p.monomial))
+    mono = adjoint(MP.monomial(p))
+    return Polynomial(Variables{B}(MP.variables(mono)), MP.exponents(mono))
 end
 
 """
@@ -336,36 +40,44 @@ function recurrence_eval(::Type{Monomial}, previous, value, degree)
     return previous[degree] * value
 end
 
-function (::Mul{Monomial})(a::MP.AbstractMonomial, b::MP.AbstractMonomial)
-    return sparse_coefficients(a * b)
+function (m::MStruct{Monomial,V,E})(a::E, b::E, ::Type{E}) where {V,E}
+    return SA.SparseCoefficients((a .+ b,), (1,))
 end
 
 SA.coeffs(p::Polynomial{Monomial}, ::FullBasis{Monomial}) = p.monomial
 
 function MP.polynomial_type(
-    ::Union{SubBasis{B,M},Type{<:SubBasis{B,M}}},
+    ::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 MP.polynomial_type(
+    ::Type{Polynomial{B,V,E}},
     ::Type{T},
-) where {B,M,T}
-    return MP.polynomial_type(FullBasis{B,M}, T)
+) where {B,V,E,T}
+    return _polynomial_type(B, V, T)
 end
 
-function MP.polynomial_type(::Type{Polynomial{B,M}}, ::Type{T}) where {B,M,T}
-    return MP.polynomial_type(FullBasis{B,M}, T)
+function keys_as_monomials(keys, mb::FullBasis)
+    return MP.monomial.(Ref(mb.map.variables), keys)
 end
+keys_as_monomials(mb::SubBasis) = keys_as_monomials(mb.keys, parent(mb))
 
 function MP.polynomial(f::Function, mb::SubBasis{Monomial})
-    return MP.polynomial(f, mb.monomials)
+    return MP.polynomial(f, keys_as_monomials(mb))
 end
 
 function MP.polynomial(Q::AbstractMatrix, mb::SubBasis{Monomial}, T::Type)
-    return MP.polynomial(Q, mb.monomials, T)
+    return MP.polynomial(Q, keys_as_monomials(mb), T)
 end
 
 function MP.coefficients(
     p::MP.AbstractPolynomialLike,
     basis::SubBasis{Monomial},
 )
-    return MP.coefficients(p, basis.monomials)
+    return MP.coefficients(p, keys_as_monomials(basis))
 end
 
 function MP.coefficients(p::MP.AbstractPolynomialLike, ::FullBasis{Monomial})
@@ -407,35 +119,46 @@ end
 #    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
+#function MA.operate!(
+#    ::SA.UnsafeAddMul{typeof(*)},
+#    a::SA.AlgebraElement{<:SA.StarAlgebra{<: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)
-MP.extdegree(basis::SubBasis{Monomial}) = MP.extdegree(basis.monomials)
+MP.mindegree(basis::SubBasis{Monomial}) = minimum(sum, basis.keys)
+MP.maxdegree(basis::SubBasis) = maximum(sum, basis.keys)
+MP.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 MP.mindegree(basis::SubBasis{Monomial}, v)
-    return MP.mindegree(basis.monomials, v)
+    i = variable_index(basis, v)
+    return minimum(Base.Fix2(getindex, i), basis.keys)
 end
 function MP.maxdegree(basis::SubBasis, v)
-    return MP.maxdegree(basis.monomials, v)
+    i = variable_index(basis, v)
+    return maximum(Base.Fix2(getindex, i), basis.keys)
 end
 function MP.extdegree(basis::SubBasis{Monomial}, v)
-    return MP.extdegree(basis.monomials, 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 MP.polynomial_type(::Type{FullBasis{B,M}}, ::Type{T}) where {T,B,M}
-    return MP.polynomial_type(M, _promote_coef(T, B))
+function MP.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 MP.polynomial_type(MP.monomial_type(V), _promote_coef(T, B))
 end
 
 _vec(v::Vector) = v
@@ -452,7 +175,7 @@ function SA.coeffs(
     if B1 === B2 && target isa FullBasis
         # The defaults initialize to zero and then sums which promotes
         # `JuMP.VariableRef` to `JuMP.AffExpr`
-        return SA.SparseCoefficients(_vec(source.monomials), _vec(cfs))
+        return SA.SparseCoefficients(_vec(source.keys), _vec(cfs))
     else
         res = SA.zero_coeffs(
             _promote_coef(_promote_coef(SA.value_type(cfs), B1), B2),
@@ -464,3 +187,54 @@ end
 
 # FIXME this assumes that the basis is invariant under adjoint
 SA.star(::SubBasis, coeffs) = SA.star.(coeffs)
+
+# TODO use Base.show_vector here, maybe by wrapping the `generator` vector
+#      into something that spits objects wrapped with the `mime` type
+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(MP.monomial, MP.variables(basis)),
+    )
+    print(io, ')')
+    return
+end
+
+function Base.show(
+    io::IO,
+    mime::MIME"text/plain",
+    basis::Union{Variables,SubBasis},
+)
+    return _show(io, mime, basis)
+end
+
+function Base.show(
+    io::IO,
+    mime::MIME"text/print",
+    basis::Union{Variables,SubBasis},
+)
+    return _show(io, mime, basis)
+end
+
+function Base.print(io::IO, basis::Union{Variables,SubBasis})
+    return show(io, MIME"text/print"(), basis)
+end
+
+function Base.show(io::IO, basis::Union{Variables,SubBasis})
+    return show(io, MIME"text/plain"(), basis)
+end
diff --git a/src/mstructures.jl b/src/mstructures.jl
new file mode 100644
index 0000000..a1b8166
--- /dev/null
+++ b/src/mstructures.jl
@@ -0,0 +1,32 @@
+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
diff --git a/src/orthogonal.jl b/src/orthogonal.jl
index be12ca4..d5b44c4 100644
--- a/src/orthogonal.jl
+++ b/src/orthogonal.jl
@@ -107,9 +107,7 @@ function univariate_orthogonal_basis(
 )
     return univariate_eval!(
         B,
-        Vector{
-            MP.polynomial_type(Polynomial{B,MP.monomial_type(variable)}, Int),
-        }(
+        Vector{_polynomial_type(B, typeof(MP.variables(variable)), Int)}(
             undef,
             degree + 1,
         ),
@@ -117,39 +115,54 @@ function univariate_orthogonal_basis(
     )
 end
 
-function _covering(::FullBasis{B,M}, monos) where {B,M}
-    to_add = collect(monos)
-    m = Set{M}(to_add)
+function _setindex(x::AbstractVector, v, i)
+    y = copy(x)
+    y[i] = v
+    return y
+end
+_setindex(x::Tuple, v, i) = Base.setindex(x, v, i)
+
+# Same as `BangBang.setindex!!`
+_setindex!(x, v, i) = Base.setindex!(x, v, i)
+_setindex!(x::Tuple, v, i) = Base.setindex(x, v, i)
+
+_increment(x, v, i) = _setindex(x, x[i] + v, i)
+_increment!(x, Δ, i) = _setindex!(x, x[i] + Δ, i)
+
+function _covering(b::FullBasis{B}, exps) where {B}
+    to_add = collect(exps)
+    m = Set{eltype(exps)}(to_add)
     while !isempty(to_add)
-        mono = pop!(to_add)
-        for v in MP.variables(mono)
+        exp = pop!(to_add)
+        for i in eachindex(exp)
             step = even_odd_separated(B) ? 2 : 1
-            vstep = v^step
-            if MP.divides(vstep, mono)
-                new_mono = MP.map_exponents(-, mono, vstep)
-                if !(new_mono in m)
-                    push!(m, new_mono)
-                    push!(to_add, new_mono)
+            if exp[i] >= step
+                new_exp = _increment(exp, -step, i)
+                if !(new_exp in m)
+                    push!(m, new_exp)
+                    push!(to_add, new_exp)
                 end
             end
         end
     end
-    return collect(m)
+    return sort!(collect(m); lt = MP.ordering(b)())
 end
 
 function explicit_basis_covering(
     full::FullBasis{BM,M},
     monos::SubBasis{B,M},
 ) where {BM<:AbstractMonomial,B<:AbstractMultipleOrthogonal,M}
-    full = FullBasis{B,M}()
-    return SubBasis{BM}(_covering(full, monos.monomials))
+    return SA.SubBasis(
+        full,
+        _covering(FullBasis{B}(MP.variables(full)), monos.keys),
+    )
 end
 
 function explicit_basis_covering(
-    full::FullBasis{B,M},
-    monos::SubBasis{<:AbstractMonomial,M},
-) where {B<:AbstractMultipleOrthogonal,M}
-    return SubBasis{B}(_covering(full, monos.monomials))
+    full::FullBasis{B,V,E},
+    monos::SubBasis{<:AbstractMonomial,V,E},
+) where {B<:AbstractMultipleOrthogonal,V,E}
+    return SA.SubBasis(full, _covering(full, monos.keys))
 end
 
 function _scalar_product_function(::Type{<:AbstractMultipleOrthogonal}, i::Int) end
@@ -198,8 +211,9 @@ function MP.coefficients(
 ) where {B<:AbstractMultipleOrthogonal,M}
     poly_p = MP.polynomial(p)
     return map(basis) do el
-        q = SA.coeffs(el, FullBasis{Monomial,M}())
-        poly_q = MP.polynomial(q)
+        full_mono = FullBasis{Monomial}(MP.variables(basis))
+        q = SA.coeffs(el, full_mono)
+        poly_q = _polynomial(full_mono, q)
         return LinearAlgebra.dot(poly_p, poly_q, B) /
                LinearAlgebra.dot(poly_q, poly_q, B)
     end
@@ -214,14 +228,12 @@ function SA.coeffs(
 end
 
 function SA.coeffs(
-    p::Polynomial{B,M},
+    p::Polynomial{B},
     ::FullBasis{Monomial},
-) where {B<:AbstractMultipleOrthogonal,M}
+) where {B<:AbstractMultipleOrthogonal}
+    mono = MP.monomial(p)
     return sparse_coefficients(
-        prod(
-            MP.powers(p.monomial);
-            init = MP.constant_monomial(M),
-        ) do (var, deg)
+        prod(MP.powers(mono); init = MP.constant_monomial(mono)) do (var, deg)
             return univariate_orthogonal_basis(B, var, deg)[deg+1]
         end,
     )
diff --git a/src/polynomial.jl b/src/polynomial.jl
index 55d2c4e..e2cc2ea 100644
--- a/src/polynomial.jl
+++ b/src/polynomial.jl
@@ -42,46 +42,54 @@ end
 abstract type AbstractMonomialIndexed end
 
 """
-    struct Polynomial{B<:AbstractMonomialIndexed,M<:MP.AbstractMonomial}
-        monomial::M
-        function Polynomial{B}(mono::MP.AbstractMonomial) where {B}
-            return new{B,typeof(mono)}(mono)
-        end
+    struct Polynomial{B<:AbstractMonomialIndexed,V,E}
+        variables::Variables{B,V}
+        exponents::E
     end
 
-Polynomial of basis `FullBasis{B,M}()` at index `monomial`.
+Polynomial of basis `FullBasis{B,V,E}(variables)` at index `exponents`.
 """
-struct Polynomial{B<:AbstractMonomialIndexed,M<:MP.AbstractMonomial}
-    monomial::M
-    function Polynomial{B}(mono::MP.AbstractMonomial) where {B}
-        return new{B,typeof(mono)}(mono)
-    end
+struct Polynomial{B<:AbstractMonomialIndexed,V,E}
+    variables::Variables{B,V}
+    exponents::E
+end
+
+function Polynomial(v::Variables{B,V}, e) where {B,V}
+    return Polynomial{B,V,typeof(e)}(v, e)
+end
+
+function Polynomial{B}(v::MP.AbstractVariable) where {B}
+    vars = MP.variables(v)
+    # MP.exponents(v) gives (1,) for DynamicPolynomials so we use
+    # a custom `map` as workaround
+    return Polynomial(Variables{B}(vars), map(_ -> 1, vars))
 end
 
+MP.exponents(p::Polynomial) = p.exponents
+
+MP.monomial(p::Polynomial) = MP.monomial(MP.variables(p), MP.exponents(p))
+
 function Base.hash(p::Polynomial{B}, u::UInt) where {B}
-    return hash(B, hash(p.monomial, u))
+    return hash(p.variables, hash(p.monomial, u))
 end
 
 function Base.isequal(p::Polynomial{B}, q::Polynomial{B}) where {B}
-    return isequal(p.monomial, q.monomial)
+    return isequal(p.variables, q.variables) &&
+           isequal(p.exponents, q.exponents)
 end
 
-Base.isone(p::Polynomial) = isone(p.monomial)
+Base.isone(p::Polynomial) = all(iszero, p.exponents)
 
 # 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
-
 function Base.:(==)(p::Polynomial{B}, q::Polynomial{B}) where {B}
-    return p.monomial == q.monomial
+    return p.variables == q.variables && p.exponents == q.exponents
 end
 
-MP.variables(p::Polynomial) = MP.variables(p.monomial)
-MP.nvariables(p::Polynomial) = MP.nvariables(p.monomial)
+MP.variables(p::Polynomial) = MP.variables(p.variables)
+MP.nvariables(p::Polynomial) = MP.nvariables(p.variables)
 
 MP.monomial_type(::Type{<:SA.SparseCoefficients{K}}) where {K} = K
 MP.polynomial(p::Polynomial) = MP.polynomial(algebra_element(p))
@@ -93,19 +101,12 @@ end
 function _algebra_element(p, ::Type{B}) where {B<:AbstractMonomialIndexed}
     return algebra_element(
         sparse_coefficients(p),
-        FullBasis{B,MP.monomial_type(typeof(p))}(),
+        FullBasis{B}(MP.variables(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),
-        FullBasis{B,promote_type(typeof(a.monomial), typeof(b.monomial))}(),
-    )
+function algebra_element(p::Polynomial{B}) where {B}
+    return _algebra_element(MP.monomial(p), B)
 end
 
 function Base.:*(a::Polynomial{B}, b::SA.AlgebraElement) where {B}
@@ -117,7 +118,13 @@ function _show(io::IO, mime::MIME, p::Polynomial{B}) where {B}
         print(io, B)
         print(io, "(")
     end
-    print(io, SA.trim_LaTeX(mime, sprint(show, mime, p.monomial)))
+    print(
+        io,
+        SA.trim_LaTeX(
+            mime,
+            sprint(show, mime, MP.monomial(p.variables.variables, p.exponents)),
+        ),
+    )
     if B != Monomial
         print(io, ")")
     end
@@ -156,29 +163,6 @@ 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;
@@ -202,7 +186,7 @@ end
 function Base.isapprox(a::SA.AlgebraElement, α::Number; kws...)
     return isapprox(
         a,
-        α * constant_algebra_element(typeof(SA.basis(a)), typeof(α));
+        α * constant_algebra_element(SA.basis(a), typeof(α));
         kws...,
     )
 end
diff --git a/src/quotient.jl b/src/quotient.jl
index 54295a8..951bc78 100644
--- a/src/quotient.jl
+++ b/src/quotient.jl
@@ -4,6 +4,7 @@ struct QuotientBasis{T,I,B<:SA.AbstractBasis{T,I},D} <: SA.ExplicitBasis{T,I}
 end
 
 implicit_basis(basis::QuotientBasis) = implicit_basis(basis.basis)
+MP.variables(basis::QuotientBasis) = MP.variables(basis.basis)
 
 Base.iterate(basis::QuotientBasis) = iterate(basis.basis)
 Base.iterate(basis::QuotientBasis, s) = iterate(basis.basis, s)
@@ -15,12 +16,15 @@ 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)
+function SA.coeffs(coeffs, b::FullBasis{Monomial}, basis::QuotientBasis)
+    return MP.coefficients(
+        MP.polynomial(values(coeffs), keys_as_monomials(keys(coeffs), b)),
+        basis,
+    )
 end
 
 function SA.coeffs(coeffs, sub::SubBasis{Monomial}, basis::QuotientBasis)
-    return SA.coeffs(MP.polynomial(coeffs, sub.monomials), parent(sub), basis)
+    return MP.coefficients(MP.polynomial(coeffs, keys_as_monomials(sub)), basis)
 end
 
 function SA.adjoint_coeffs(
@@ -28,11 +32,10 @@ function SA.adjoint_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))
-        )
+    return map(src) do poly
+        return sum(MP.terms(rem(MP.polynomial(poly), dest.divisor))) do t
+            return MP.coefficient(t) *
+                   coeffs[SA.key_index(dest.basis, MP.exponents(t))]
+        end
     end
 end
diff --git a/src/scaled.jl b/src/scaled.jl
index fde028e..ce05766 100644
--- a/src/scaled.jl
+++ b/src/scaled.jl
@@ -33,15 +33,16 @@ Foundations of Computational Mathematics 7.2 (2007): 229-244.
 """
 struct ScaledMonomial <: AbstractMonomial end
 
-function (::Mul{ScaledMonomial})(a::MP.AbstractMonomial, b::MP.AbstractMonomial)
-    mono = a * b
+function (m::MStruct{ScaledMonomial,V,E})(a::E, b::E, ::Type{E}) where {V,E}
+    @assert a.variables == b.variables
+    exp = a.exponents .+ b.exponents
     α = 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))
+        eachindex(exp);
+        init = inv(binomial(sum(mono), sum(a.exponents))),
+    ) do i
+        return binomial(exp[i], a.exponents[i])
     end
-    return sparse_coefficients(MP.term(√α, mono))
+    return SA.SparseCoefficients((exp,), (1,))
 end
 
 function SA.coeffs(p::Polynomial{ScaledMonomial}, ::FullBasis{Monomial})
@@ -57,8 +58,8 @@ _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),
-        basis.monomials,
+        i -> scaling(basis.keys[i]) * f(i),
+        keys_as_monomials(basis),
     )
 end
 
@@ -69,10 +70,8 @@ function Base.promote_rule(
     return SubBasis{Monomial,M,V}
 end
 
-function scaling(m::MP.AbstractMonomial)
-    return √(
-        factorial(MP.degree(m)) / prod(factorial, MP.exponents(m); init = 1),
-    )
+function scaling(exp)
+    return √(factorial(sum(exp)) / prod(factorial, exp; init = 1),)
 end
 unscale_coef(t::MP.AbstractTerm) = MP.coefficient(t) / scaling(MP.monomial(t))
 function SA.coeffs(
@@ -87,7 +86,7 @@ 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))
+    return MP.coefficients(p, keys_as_monomials(basis)) ./ scaling.(basis.keys)
 end
 function SA.coeffs(
     p::MP.AbstractPolynomialLike,
@@ -105,12 +104,8 @@ function SA.coeffs!(
 )
     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),
-        )
+        exp = source[k].exponents
+        SA.unsafe_push!(res, exponents_index(target, exp), v * scaling(exp))
     end
     MA.operate!(SA.canonical, res)
     return res
@@ -124,12 +119,8 @@ function SA.coeffs!(
 )
     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),
-        )
+        exp = source[k].exponents
+        SA.unsafe_push!(res, exponents_index(target, exp), v / scaling(exp))
     end
     MA.operate!(SA.canonical, res)
     return res
diff --git a/src/trigonometric.jl b/src/trigonometric.jl
index 7b6ce57..f0bd8ea 100644
--- a/src/trigonometric.jl
+++ b/src/trigonometric.jl
@@ -21,25 +21,23 @@ _id(d) = div(d + 1, 2)
 # If a > b
 # sin(a) cos(b) = sin(a + b) + sin(a - b)
 # sin(a) cos(b) = sin(a + b) + sin(a - b)
-function univariate_mul!(::Mul{Trigonometric}, terms, var, a, b)
+function univariate_mul!(::Type{Trigonometric}, exps, coefs, var, a, b)
     @assert !iszero(a)
     @assert !iszero(b)
-    I = eachindex(terms)
     da = _id(a)
     db = _id(b)
-    for i in I
+    for i in eachindex(exps)
         if _is_cos(a) == _is_cos(b)
             # Chebyshev first kind
-            mono = MP.monomial(terms[i]) * var^(_cos_id(da + db))
-            terms[i] = MA.mul!!(terms[i], var^_cos_id(abs(da - db)))
-            terms[i] = MA.operate!!(/, terms[i], 2)
-            α = MA.copy_if_mutable(MP.coefficient(terms[i]))
-            push!(terms, MP.term(α, mono))
+            push!(exps, _increment(exps[i], _cos_id(da + db), var))
+            exps[i] = _increment!(exps[i], _cos_id(abs(da - db)), var)
+            coefs[i] = MA.operate!!(/, coefs[i], 2)
+            push!(coefs, MA.copy_if_mutable(coefs[i]))
             # cos(a + b) = cos(a) cos(b) - sin(a) sin(b)
             # cos(a - b) = cos(a) cos(b) + sin(a) sin(b)
             # cos(a - b) - cos(a + b)
             if _is_sin(a)
-                terms[end] = MA.operate!!(*, terms[end], -1)
+                coefs[end] = MA.operate!!(*, coefs[end], -1)
             end
         else
             if _is_cos(a)
@@ -48,8 +46,8 @@ function univariate_mul!(::Mul{Trigonometric}, terms, var, a, b)
             # sin(da) * cos(db)
             if da == db
                 # sin(da) * cos(da) = sin(2da) / 2
-                terms[i] = MA.mul!!(terms[i], var^_cos_id(da + db))
-                terms[i] = MA.operate!!(/, terms[i], 2)
+                exps[i] = _increment!(exps[i], _cos_id(da + db), var)
+                coefs[i] = MA.operate!!(/, coefs[i], 2)
             else
                 # Using
                 # sin(a + b) = sin(a) cos(b) + cos(a) sin(b)
@@ -58,13 +56,12 @@ function univariate_mul!(::Mul{Trigonometric}, terms, var, a, b)
                 # sin(a) cos(b) = (sin(a + b) + sin(a - b)) / 2
                 # If a < b
                 # sin(a) cos(b) = (sin(b + a) - sin(b - a)) / 2
-                mono = MP.monomial(terms[i]) * var^(_sin_id(da + db))
-                terms[i] = MA.mul!!(terms[i], var^_sin_id(abs(da - db)))
-                terms[i] = MA.operate!!(/, terms[i], 2)
-                α = MA.copy_if_mutable(MP.coefficient(terms[i]))
-                push!(terms, MP.term(α, mono))
+                push!(exps, _increment(exps[i], _sin_id(da + db), var))
+                exps[i] = _increment!(exps[i], _sin_id(abs(da - db)), var)
+                coefs[i] = MA.operate!!(/, coefs[i], 2)
+                push!(coefs, MA.copy_if_mutable(MP.coefficient(terms[i])))
                 if da < db
-                    terms[i] = MA.operate!!(*, terms[i], -1)
+                    coefs[i] = MA.operate!!(*, coefs[i], -1)
                 end
             end
         end
@@ -85,8 +82,7 @@ function recurrence_eval(
     d = _id(degree)
     if _is_cos(degree)
         # Chebyshev first order
-        return 2 * value * previous[_cos_id(d - 1)+1] -
-               previous[_cos_id(d - 2)+1]
+        return 2 * value * previous[_cos_id(d-1)+1] - previous[_cos_id(d-2)+1]
     else
         return sqrt(1 - previous[degree]^2)
     end
diff --git a/test/Project.toml b/test/Project.toml
index e253ba2..40a56d6 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -2,7 +2,11 @@
 DynamicPolynomials = "7c1d4256-1411-5781-91ec-d7bc3513ac07"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MultivariateBases = "be282fd4-ad43-11e9-1d11-8bd9d7e43378"
+MultivariatePolynomials = "102ac46a-7ee4-5c85-9060-abc95bfdeaa3"
 MutableArithmetics = "d8a4904e-b15c-11e9-3269-09a3773c0cb0"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 StarAlgebras = "0c0c59c1-dc5f-42e9-9a8b-b5dc384a6cd1"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+
+[compat]
+DynamicPolynomials = "0.6.3"
diff --git a/test/chebyshev.jl b/test/chebyshev.jl
index 68ebdf1..6126f58 100644
--- a/test/chebyshev.jl
+++ b/test/chebyshev.jl
@@ -1,17 +1,24 @@
 using Test
+import StarAlgebras as SA
 using MultivariateBases
+import MultivariateBases as MB
 using DynamicPolynomials
 
 @testset "StarAlgebras" begin
     @polyvar x
-    a = MB.Polynomial{MB.Chebyshev}(x)
+    full = MB.FullBasis{MB.Chebyshev}(x)
+    a = MB.algebra_element(MB.sparse_coefficients(1 // 1 * x), full)
+    @test a == MB.convert_basis(full, MB.algebra_element(x))
     b = a * a
-    @test b.coeffs == MB.sparse_coefficients(1 // 2 + 1 // 2 * x^2)
+    expected =
+        MB.algebra_element(MB.sparse_coefficients((1 // 2) * (1 + x^2)), full)
+    @test b == expected
     c = b * b
-    @test c.coeffs ==
+    @test SA.coeffs(c) ==
           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)
+    c1 = MB.convert_basis(full, MB.algebra_element(constant_monomial(x)))
+    c2 = MB.convert_basis(full, MB.algebra_element(x^0))
+    @test a * c1 == a * c2
 end
 
 @testset "Orthogonal" begin
diff --git a/test/fixed.jl b/test/fixed.jl
index e90849a..9a071a7 100644
--- a/test/fixed.jl
+++ b/test/fixed.jl
@@ -4,46 +4,47 @@ using MultivariateBases
 const MB = MultivariateBases
 using DynamicPolynomials
 
-@testset "FixedBasis" begin
-    @polyvar x y
-    x_term = MB.term_element(1, MB.Polynomial{MB.Monomial}(x))
-    y_term = MB.term_element(1, MB.Polynomial{MB.Monomial}(y))
-    p1 = x_term + im * y_term
-    q1 = x_term - im * y_term
-    p2 = im * x_term - 2y_term
-    q2 = -im * x_term - 2y_term
-    fixed = MB.FixedBasis([p1, p2])
-    @test length(fixed) == 2
-    @test fixed[1] ≈ p1
-    @test fixed[2] ≈ p2
-    @test sprint(show, fixed) == "FixedBasis([$p1, $p2])"
-
-    semi =
-        MB.SemisimpleBasis([MB.FixedBasis([p1, p2]), MB.FixedBasis([q1, q2])])
-    @test length(semi) == 2
-    @test sprint(show, semi) ==
-          "Semisimple basis with 2 simple sub-bases:\n  FixedBasis([$p1, $p2])\n  FixedBasis([$q1, $q2])"
-    mult = semi[1]
-    @test all(mult.elements .≈ [p1, q1])
-    smult = SA.star(mult)
-    @test all(smult.elements .≈ [q1, p1])
-
-    elem, state = iterate(semi)
-    @test elem == semi[1]
-    elem, state = iterate(semi, state)
-    @test elem == semi[2]
-    @test isnothing(iterate(semi, state))
-
-    res = zero(Complex{Int}, MB.algebra(MB.FullBasis{Monomial,typeof(x * y)}()))
-    MA.operate!(SA.UnsafeAddMul(*), res, semi[1], semi[2], true)
-    MA.operate!(SA.canonical, res)
-    @test res == p1 * p2 + q1 * q2
-    MA.operate!(SA.UnsafeAddMul(*), res, semi[1], semi[2], false)
-    @test res == p1 * p2 + q1 * q2
-    MA.operate!(SA.UnsafeAddMul(*), res, semi[2], semi[1], -1)
-    MA.operate!(SA.canonical, res)
-    @test iszero(res)
-    MA.operate!(SA.UnsafeAddMul(*), res, semi[2], semi[2], -1)
-    MA.operate!(SA.canonical, res)
-    @test res == -p2^2 - q2^2
-end
+#@testset "FixedBasis" begin
+#    @polyvar x y
+#    x_term = MB.term_element(1, MB.Polynomial{MB.Monomial}(x))
+#    y_term = MB.term_element(1, MB.Polynomial{MB.Monomial}(x))
+#    p1 = x_term + im * y_term
+#    q1 = x_term - im * y_term
+#    p2 = im * x_term - 2y_term
+#    q2 = -im * x_term - 2y_term
+#    fixed = SA.FixedBasis([p1, p2])
+#    @test length(fixed) == 2
+#    @test fixed[1] ≈ p1
+#    @test fixed[2] ≈ p2
+#    @test sprint(show, fixed) == "FixedBasis([$p1, $p2])"
+#
+#    semi =
+#        MB.SemisimpleBasis([MB.FixedBasis([p1, p2]), MB.FixedBasis([q1, q2])])
+#    @test length(semi) == 2
+#    @test sprint(show, semi) ==
+#          "Semisimple basis with 2 simple sub-bases:\n  FixedBasis([$p1, $p2])\n  FixedBasis([$q1, $q2])"
+#    mult = semi[1]
+#    @test all(mult.elements .≈ [p1, q1])
+#    smult = SA.star(mult)
+#    @test all(smult.elements .≈ [q1, p1])
+#
+#    elem, state = iterate(semi)
+#    @test elem == semi[1]
+#    elem, state = iterate(semi, state)
+#    @test elem == semi[2]
+#    @test isnothing(iterate(semi, state))
+#
+#    res = zero(Complex{Int}, MB.algebra(MB.FullBasis{Monomial,typeof(x * y)}()))
+#    MA.operate!(SA.UnsafeAddMul(*), res, semi[1], semi[2], true)
+#    MA.operate!(SA.canonical, res)
+#    @test res == p1 * p2 + q1 * q2
+#    MA.operate!(SA.UnsafeAddMul(*), res, semi[1], semi[2], false)
+#    @test res == p1 * p2 + q1 * q2
+#    MA.operate!(SA.UnsafeAddMul(*), res, semi[2], semi[1], -1)
+#    MA.operate!(SA.canonical, res)
+#    @test iszero(res)
+#    MA.operate!(SA.UnsafeAddMul(*), res, semi[2], semi[2], -1)
+#    MA.operate!(SA.canonical, res)
+#    @test res == -p2^2 - q2^2
+#end
+#
diff --git a/test/hermite.jl b/test/hermite.jl
index cd28b62..e8c951c 100644
--- a/test/hermite.jl
+++ b/test/hermite.jl
@@ -45,7 +45,4 @@ end
             1.0,
         ]),
     )
-    M = typeof(x^2)
-    mono = MB.FullBasis{MB.Monomial,M}()
-    basis = MB.FullBasis{MB.PhysicistsHermite,M}()
 end
diff --git a/test/lagrange.jl b/test/lagrange.jl
index db988e7..56857d1 100644
--- a/test/lagrange.jl
+++ b/test/lagrange.jl
@@ -25,8 +25,8 @@ function _test(B::Type)
     lag = MB.explicit_basis_covering(implicit, sub)
     @test typeof(lag) == MB.explicit_basis_type(typeof(implicit))
     a = MB.algebra_element(
-        SA.SparseCoefficients(collect(monos), coeffs),
-        MB.FullBasis{B,typeof(prod(x))}(),
+        SA.SparseCoefficients(collect(exponents.(monos)), coeffs),
+        MB.FullBasis{B}(x),
     )
     @test SA.coeffs(a, lag) ≈ SA.coeffs(coeffs, sub, lag)
     @polyvar z
diff --git a/test/monomial.jl b/test/monomial.jl
index 24e569e..a2a32ec 100644
--- a/test/monomial.jl
+++ b/test/monomial.jl
@@ -6,14 +6,10 @@ using DynamicPolynomials
 
 @testset "StarAlgebras" begin
     @polyvar x
-    a = MB.Polynomial{MB.Monomial}(x)
-    @test !isone(a)
-    b = a * a
-    @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)
-    o = MB.Polynomial{MB.Monomial}(constant_monomial(x^2))
+    vars = MB.Variables{MB.Monomial}(variables(x))
+    @test vars(exponents(x^2)) == vars(exponents(x^2))
+    @test vars(exponents(x^3)) != vars(exponents(x^2))
+    o = vars(exponents(constant_monomial(x^2)))
     @test isone(o)
 end
 
@@ -26,7 +22,7 @@ end
     @test coefficients(x + 4y, basis) == [4, 1]
     @test polynomial(basis[1]) == y
     @test polynomial(basis[2]) == x
-    @test basis.monomials == [y, x]
+    @test MB.keys_as_monomials(basis) == [y, x]
     @test polynomial.(collect(basis)) == [y, x]
     @test variables(basis) == [x, y]
     @test nvariables(basis) == 2
@@ -37,6 +33,7 @@ end
           "SubBasis{Monomial}([y, x])"
     @test sprint(print, basis) == "SubBasis{Monomial}([y, x])"
 end
+
 @testset "Affine" begin
     # It will be sorted and 1 will be moved at the end
     basis = SubBasis{MB.Monomial}([1, x, y])
@@ -46,6 +43,7 @@ end
     @test coefficients(9 + x + 4y, basis) == [9, 4, 1]
     @test MB.explicit_basis(poly) == basis
 end
+
 @testset "Quadratic" begin
     basis = SubBasis{MB.Monomial}([x^2, x * y, y^2])
     @test polynomial_type(basis, Int) == polynomial_type(x, Int)
@@ -54,6 +52,7 @@ end
     @test sprint(show, basis) == "SubBasis{Monomial}([y², xy, x²])"
     @test sprint(print, basis) == "SubBasis{Monomial}([y^2, x*y, x^2])"
 end
+
 @testset "merge_bases" begin
     basis1 = SubBasis{MB.Monomial}([x^2, y^2])
     basis2 = SubBasis{MB.Monomial}([x * y, y^2])
@@ -67,14 +66,13 @@ end
     @test extdegree(basis2, x) == (0, 1)
     @test extdegree(basis2, y) == (1, 2)
     basis, I1, I2 = MultivariateBases.merge_bases(basis1, basis2)
-    @test basis.monomials == [y^2, x * y, x^2]
+    @test MB.keys_as_monomials(basis) == [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
+        @test SA.key_index(basis, basis.keys[i]) == i
         for (I, b) in [(I1, basis1), (I2, basis2)]
-            idx = monomial_index(b, basis.monomials[i])
+            idx = SA.key_index(b, basis.keys[i])
             if iszero(I[i])
                 @test isnothing(idx)
             else
diff --git a/test/runtests.jl b/test/runtests.jl
index 9c7f370..c5ad1f3 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -16,7 +16,7 @@ end
 function _test_basis(basis)
     B = typeof(basis)
     @test typeof(MB.algebra(basis)) == MA.promote_operation(MB.algebra, B)
-    @test typeof(MB.constant_algebra_element(B, 1)) ==
+    @test typeof(MB.constant_algebra_element(basis, 1)) ==
           MB.constant_algebra_element_type(B, Int)
 end
 
@@ -33,6 +33,7 @@ Base.:*(::Float64, ::TypeA) = TypeB()
 Base.:*(::TypeA, ::Float64) = TypeB()
 Base.:*(::Float64, ::TypeB) = TypeB()
 Base.:*(::TypeB, ::Float64) = TypeB()
+Base.:+(::TypeA, ::TypeA) = TypeB()
 Base.:/(::TypeA, ::Float64) = TypeB()
 Base.:/(::TypeB, ::Float64) = TypeB()
 Base.:+(::TypeB, ::TypeB) = TypeB()
@@ -41,11 +42,13 @@ Base.convert(::Type{TypeB}, ::TypeA) = TypeB()
 
 function api_test(B::Type{<:MB.AbstractMonomialIndexed}, degree)
     @polyvar x[1:2]
-    M = typeof(prod(x))
-    full_basis = FullBasis{B,M}()
+    full_basis = FullBasis{B}(x)
     _test_basis(full_basis)
-    @test sprint(show, MB.algebra(full_basis)) ==
-          "Polynomial algebra of $B basis"
+    # FIXME `AbstractStarAlgebra` is currently calling `print` anyway
+    #@test sprint(show, MIME"text/plain"(), MB.algebra(full_basis)) ==
+    #      "*-algebra of $B polynomials in the variables [x[1], x[2]]"
+    @test sprint(print, MB.algebra(full_basis)) ==
+          "*-algebra of $B polynomials in the variables [x[1], x[2]]"
     for basis in [
         maxdegree_basis(full_basis, x, degree),
         explicit_basis_covering(
@@ -64,16 +67,13 @@ function api_test(B::Type{<:MB.AbstractMonomialIndexed}, degree)
         _test_basis(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)
+            poly =
+                MB.Polynomial(MB.Variables{B}(variables(basis)), basis.keys[i])
             @test basis[i] == poly
             @test basis[poly] == i
         end
         n = binomial(2 + degree, 2)
         @test length(basis) == n
-        @test firstindex(basis) == 1
-        @test lastindex(basis) == n
-        @test typeof(copy(basis)) == typeof(basis)
         @test nvariables(basis) == 2
         @test variables(basis) == x
         @test monomial_type(typeof(basis)) == monomial_type(x[1])
@@ -86,7 +86,7 @@ function api_test(B::Type{<:MB.AbstractMonomialIndexed}, degree)
         @test SA.star(a) == a
         if B == Chebyshev || B == ScaledMonomial
             mono = explicit_basis_covering(
-                FullBasis{Monomial,eltype(basis.monomials)}(),
+                FullBasis{Monomial}(variables(basis)),
                 basis,
             )
             a = MB.algebra_element(fill(TypeA(), length(basis)), mono)
@@ -95,9 +95,9 @@ function api_test(B::Type{<:MB.AbstractMonomialIndexed}, degree)
         end
     end
     mono = x[1]^2 * x[2]^3
-    p = MB.Polynomial{B}(mono)
-    @test full_basis[p] == mono
-    @test full_basis[mono] == p
+    p = MB.Polynomial(MB.Variables{B}(variables(mono)), exponents(mono))
+    @test full_basis[p] == exponents(mono)
+    @test full_basis[exponents(mono)] == p
     @test polynomial_type(mono, B == Monomial ? TypeA : TypeB) ==
           polynomial_type(typeof(p), TypeA)
     a = MB.algebra_element(p)
@@ -106,22 +106,25 @@ function api_test(B::Type{<:MB.AbstractMonomialIndexed}, degree)
     @test typeof(polynomial(a)) == polynomial_type(typeof(p), Int)
     @test a ≈ a
     if B == MB.Monomial
-        @test a ≈ p.monomial
-        @test p.monomial ≈ a
+        @test a ≈ monomial(p)
+        @test monomial(p) ≈ a
     else
-        @test !(a ≈ p.monomial)
-        @test !(p.monomial ≈ a)
+        @test !(a ≈ monomial(p))
+        @test !(monomial(p) ≈ 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))
+    @test sprint(show, p) == _wrap(sprint(show, monomial(p)))
+    @test sprint(print, p) == _wrap(sprint(print, monomial(p)))
     mime = MIME"text/latex"()
     @test sprint(show, mime, p) ==
           "\$\$ " *
-          _wrap(MB.SA.trim_LaTeX(mime, sprint(show, mime, p.monomial))) *
+          _wrap(MB.SA.trim_LaTeX(mime, sprint(show, mime, monomial(p)))) *
           " \$\$"
     const_mono = constant_monomial(prod(x))
-    const_poly = MB.Polynomial{B}(const_mono)
+    const_poly = MB.Polynomial(
+        MB.Variables{B}(variables(const_mono)),
+        exponents(const_mono),
+    )
     const_alg_el = MB.algebra_element(const_poly)
     for other in (const_mono, 1, const_alg_el)
         @test _test_op(+, other, const_alg_el) ≈ _test_op(*, 2, other)
@@ -157,11 +160,11 @@ function univ_orthogonal_test(
     kwargs...,
 )
     @polyvar x
-    basis = maxdegree_basis(FullBasis{B,monomial_type(x)}(), [x], 4)
+    basis = maxdegree_basis(FullBasis{B}(x), [x], 4)
     for i in eachindex(basis)
         p_i = polynomial(basis[i])
         @test isapprox(dot(p_i, p_i, B), univ(maxdegree(p_i)); kwargs...)
-        for j in 1:i-1
+        for j in 1:(i-1)
             @test isapprox(dot(p_i, polynomial(basis[j]), B), 0.0; kwargs...)
         end
     end
@@ -179,7 +182,7 @@ function orthogonal_test(
 
     @testset "Univariate $var" for (var, univ) in
                                    [(x, univariate_x), (y, univariate_y)]
-        basis = maxdegree_basis(FullBasis{B,monomial_type(var)}(), (var,), 4)
+        basis = maxdegree_basis(FullBasis{B}(var), (var,), 4)
         for i in 1:5
             @test polynomial(basis[i]) == univ[i]
         end
@@ -187,8 +190,8 @@ function orthogonal_test(
 
     @testset "explicit_basis_covering" begin
         basis = explicit_basis_covering(
-            FullBasis{B,typeof(x * y)}(),
-            SubBasis{MB.Monomial}(monomial_vector([x^2 * y, y^2])),
+            FullBasis{B}(x * y),
+            FullBasis{MB.Monomial}(x * y)[[x^2 * y, y^2]],
         )
         if even_odd_separated
             exps = [(0, 0), (0, 1), (0, 2), (2, 1)]
@@ -200,8 +203,8 @@ function orthogonal_test(
                   univariate_x[exps[i][1]+1] * univariate_y[exps[i][2]+1]
         end
         basis = explicit_basis_covering(
-            FullBasis{B,typeof(x^2)}(),
-            SubBasis{MB.Monomial}(monomial_vector([x^4, x^2, x])),
+            FullBasis{B}(x^2),
+            FullBasis{Monomial}(x^2)[[x^4, x^2, x]],
         )
         if even_odd_separated
             exps = [0, 1, 2, 4]
@@ -237,7 +240,7 @@ function coefficient_test(
     @polyvar x y
     p = x^4 * y^2 + x^2 * y^4 - 3 * x^2 * y^2 + 1
     basis = explicit_basis_covering(
-        FullBasis{B,typeof(x * y)}(),
+        FullBasis{B}(variables(x * y)),
         SubBasis{MB.Monomial}(monomials(p)),
     )
     coefficient_test(basis, p, coefs; kwargs...)
diff --git a/test/scaled.jl b/test/scaled.jl
index 4975b21..bf45b2d 100644
--- a/test/scaled.jl
+++ b/test/scaled.jl
@@ -24,9 +24,11 @@ end
     @test coefficients(p, basis) == [9, 4 / √2, 1]
     a = MB.algebra_element(p)
     @test SA.coeffs(a, basis) == [9, 4 / √2, 1]
-    full = MB.FullBasis{MB.ScaledMonomial,monomial_type(p)}()
-    @test SA.coeffs(a, full) ==
-          SA.SparseCoefficients([y^2, x * y, x^2], [9, 4 / √2, 1])
+    full = MB.FullBasis{MB.ScaledMonomial}(variables(p))
+    @test SA.coeffs(a, full) == SA.SparseCoefficients(
+        exponents.(monovec([y^2, x * y, x^2])),
+        [9, 4 / √2, 1],
+    )
     @test polynomial(basis[1]) == y^2
     @test polynomial(basis[2]) == √2 * x * y
     @test polynomial(basis[3]) == x^2
diff --git a/test/trigonometric.jl b/test/trigonometric.jl
index 1c1fda7..d209e7a 100644
--- a/test/trigonometric.jl
+++ b/test/trigonometric.jl
@@ -4,12 +4,17 @@ using DynamicPolynomials
 
 @testset "StarAlgebras" begin
     @polyvar x
-    a = MB.Polynomial{MB.Trigonometric}(x)
+    full = MB.FullBasis{MB.Trigonometric}(x)
+    a = MB.algebra_element(MB.sparse_coefficients(1//1 * x), full)
     b = a * a
     @test b.coeffs == MB.sparse_coefficients(1 // 2 + 1 // 2 * x^3)
     c = b * b
     @test c.coeffs ==
           MB.sparse_coefficients(3 // 8 + 1 // 2 * x^3 + 1 // 8 * x^7)
-    @test a * MB.Polynomial{MB.Trigonometric}(constant_monomial(typeof(x))) ==
-          a * MB.Polynomial{MB.Trigonometric}(x^0)
+    d = MB.algebra_element(
+        MB.sparse_coefficients(1//1 * constant_monomial(typeof(x))),
+        full,
+    )
+    e = MB.algebra_element(MB.sparse_coefficients(1//1 * x^0), full)
+    @test a * d == a * e
 end