Refactor

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"

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")

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