From 121291b7292450b6d642f0e5822fefbb3c9b57aa Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Fri, 13 Jun 2025 14:25:10 +0200 Subject: [PATCH 01/19] Refactor --- src/MultivariateBases.jl | 29 +--- src/bases.jl | 228 ++++++++++++++++++++++++++ src/iterator.jl | 57 +++++++ src/monomial.jl | 336 +++++---------------------------------- src/polynomial.jl | 7 +- src/trash.jl | 49 ++++++ test/Project.toml | 1 + 7 files changed, 381 insertions(+), 326 deletions(-) create mode 100644 src/bases.jl create mode 100644 src/iterator.jl create mode 100644 src/trash.jl diff --git a/src/MultivariateBases.jl b/src/MultivariateBases.jl index 1ff3ba6..743a80c 100644 --- a/src/MultivariateBases.jl +++ b/src/MultivariateBases.jl @@ -14,33 +14,12 @@ 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 -function MA.promote_operation( - ::typeof(SA.basis), - ::Type{<:Algebra{B}}, -) where {B} - return B -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 -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") +MP.monomial_type(::Type{<:SA.StarAlgebra{O}}) where {O} = MP.monomial_type(O) +function MP.polynomial_type(::Type{SA.StarAlgebra{O,S,M}}, ::Type{T}) where {O,S,M,T} + return MP.polynomial_type(M, T) end +include("bases.jl") include("monomial.jl") include("scaled.jl") diff --git a/src/bases.jl b/src/bases.jl new file mode 100644 index 0000000..8e9094d --- /dev/null +++ b/src/bases.jl @@ -0,0 +1,228 @@ +struct Variables{B,V} + variables::V +end + +function (v::Variables{B})(exponents) where {B} + return Polynomial{B}(MP.monomial(v.variables, exponents)) +end + +const FullBasis{B,M} = SA.MappedBasis{Polynomial{B,M}} +const SubBasis{B,M} = SA.SubBasis{Polynomial{B,M}} + +function FullBasis{B,M}(vars) where {B,M} + O = typeof(MP.ordering(vars)) + return SA.MappedBasis{Polynomial{B,M}}( + MP.ExponentsIterator{O}(vars), + Variables{B,typeof(vars)}(vars), + MP.exponents, + ) +end + +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 + +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 monomial_index(basis::SubBasis, mono::MP.AbstractMonomial) + @assert MP.variables(basis) == MP.variables(mono) + return get(basis, MP.exponents(mono), nothing) +end + +function explicit_basis_covering(::FullBasis{B}, target::SubBasis{B}) where {B} + return SubBasis{B}(target.monomials) +end + +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 \ No newline at end of file diff --git a/src/iterator.jl b/src/iterator.jl new file mode 100644 index 0000000..e6fba8d --- /dev/null +++ b/src/iterator.jl @@ -0,0 +1,57 @@ +import StarAlgebras as SA + +struct ExponentsIterator end +Base.eltype(::Type{ExponentsIterator}) = NTuple{2,Int} +Base.IteratorSize(::Type{ExponentsIterator}) = Base.IsInfinite() + +function Base.iterate(::ExponentsIterator) + z = (0, 0) + return z, (z, 0) +end + +function Base.iterate(::ExponentsIterator, state) + z, deg = state + if iszero(z[2]) + deg += 1 + z = (0, deg) + else + z = (z[1] + 1, z[2] - 1) + end + return z, (z, deg) +end + +function grlex(a::NTuple{2,Int}, b::NTuple{2,Int}) + return isless((sum(a), a), (sum(b), b)) +end + +struct Monomial + exponents::NTuple{2,Int} +end + +Base.one(::Monomial) = Monomial((0, 0)) +Base.:*(a::Monomial, b::Monomial) = Monomial(a.exponents .+ b.exponents) + +monomial(exp) = Monomial(exp) +exponents(mono::Monomial) = mono.exponents + +exps = ExponentsIterator() +SA.comparable(::ExponentsIterator) = grlex + +basis = SA.MappedBasis(exps, monomial, exponents) +mstr = SA.DiracMStructure(basis, *) +object = Monomial((0, 0)) +alg = SA.StarAlgebra(object, mstr) +one(alg) +a = SA.AlgebraElement( + SA.SparseCoefficients( + collect(Iterators.take(exps, 3)), + [2, -1, 3], + ), + alg, +) +c = a * a +import MutableArithmetics as MA +MA.operate!(SA.canonical, c) +@edit MA.operate!(SA.canonical, c) +@edit MA.operate!(SA.canonical, SA.coeffs(c)) +@edit MA.operate!(SA.canonical, SA.coeffs(c), isless) \ No newline at end of file diff --git a/src/monomial.jl b/src/monomial.jl index a6d3d23..f44dee4 100644 --- a/src/monomial.jl +++ b/src/monomial.jl @@ -1,301 +1,3 @@ -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( @@ -464,3 +166,41 @@ 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) + 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 \ No newline at end of file diff --git a/src/polynomial.jl b/src/polynomial.jl index 55d2c4e..1425aea 100644 --- a/src/polynomial.jl +++ b/src/polynomial.jl @@ -102,9 +102,10 @@ function algebra_element(p::Polynomial{B,M}) where {B,M} end function Base.:*(a::Polynomial{B}, b::Polynomial{B}) where {B} + M = promote_type(typeof(a.monomial), typeof(b.monomial)) return algebra_element( - Mul{B}()(a.monomial, b.monomial), - FullBasis{B,promote_type(typeof(a.monomial), typeof(b.monomial))}(), + Mul{B,M}()(a.monomial, b.monomial), + FullBasis{B,M}(), ) end @@ -156,7 +157,7 @@ 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 +struct Mul{B<:AbstractMonomialIndexed,M} <: SA.MultiplicativeStructure{Polynomial{B,M},M} end function MA.operate_to!( p::MP.AbstractPolynomial, diff --git a/src/trash.jl b/src/trash.jl new file mode 100644 index 0000000..2fa91ce --- /dev/null +++ b/src/trash.jl @@ -0,0 +1,49 @@ +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 + +# 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) + +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 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 + +const MonomialIndexedBasis{B,M} = Union{SubBasis{B,M},FullBasis{B,M}} diff --git a/test/Project.toml b/test/Project.toml index e253ba2..47026ac 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -2,6 +2,7 @@ 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" From 436aebdb3cf5605944abccb1b884261679309b87 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Fri, 13 Jun 2025 22:57:58 +0200 Subject: [PATCH 02/19] Fixes --- src/MultivariateBases.jl | 14 ++++++++- src/arithmetic.jl | 4 +-- src/bases.jl | 43 +++++++++------------------ src/chebyshev.jl | 50 ++++++++++++++++--------------- src/fixed.jl | 23 --------------- src/monomial.jl | 27 +++++++++-------- src/orthogonal.jl | 38 ++++++++++++++++-------- src/polynomial.jl | 64 ++++++++++++++++------------------------ src/scaled.jl | 15 +++++----- src/trash.jl | 2 +- src/trigonometric.jl | 5 ++-- test/chebyshev.jl | 8 +++-- test/runtests.jl | 4 +-- 13 files changed, 138 insertions(+), 159 deletions(-) diff --git a/src/MultivariateBases.jl b/src/MultivariateBases.jl index 743a80c..5d1b564 100644 --- a/src/MultivariateBases.jl +++ b/src/MultivariateBases.jl @@ -8,6 +8,18 @@ export FullBasis, SubBasis export maxdegree_basis, explicit_basis_covering, empty_basis, monomial_index include("interface.jl") +struct Variables{B,V} + variables::V +end + +Variables{B}(vars) where {B} = Variables{B,typeof(vars)}(vars) + +constant_monomial_exponents(v::Variables) = map(_ -> 0, v.variables) + +function (v::Variables)(exponents) + return Polynomial(v, exponents) +end + export AbstractMonomialIndexed, Monomial, ScaledMonomial include("polynomial.jl") MP.monomial_type(::Type{<:SA.AlgebraElement{A}}) where {A} = MP.monomial_type(A) @@ -48,7 +60,7 @@ 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) + return SA.StarAlgebra(Variables{B}(MP.variables(basis)), basis) end function MA.promote_operation( diff --git a/src/arithmetic.jl b/src/arithmetic.jl index d94e626..1be0e15 100644 --- a/src/arithmetic.jl +++ b/src/arithmetic.jl @@ -1,6 +1,6 @@ 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{<:Variables} for op in [:+, :-, :*] @eval begin @@ -79,7 +79,7 @@ end function term_element(α, p::Polynomial{B,M}) where {B,M} return algebra_element( sparse_coefficients(MP.term(α, p.monomial)), - FullBasis{B,M}(), + FullBasis{B,M}(MP.variables(p)), ) end diff --git a/src/bases.jl b/src/bases.jl index 8e9094d..807cc42 100644 --- a/src/bases.jl +++ b/src/bases.jl @@ -1,21 +1,12 @@ -struct Variables{B,V} - variables::V -end - -function (v::Variables{B})(exponents) where {B} - return Polynomial{B}(MP.monomial(v.variables, exponents)) -end - -const FullBasis{B,M} = SA.MappedBasis{Polynomial{B,M}} -const SubBasis{B,M} = SA.SubBasis{Polynomial{B,M}} +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}} function FullBasis{B,M}(vars) where {B,M} - O = typeof(MP.ordering(vars)) - return SA.MappedBasis{Polynomial{B,M}}( - MP.ExponentsIterator{O}(vars), - Variables{B,typeof(vars)}(vars), - MP.exponents, - ) + O = typeof(MP.ordering(first(vars))) # TODO upstream MP.ordering(vars) that handles differently tuples and vectors + 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 MP.monomial_type(::Type{<:FullBasis{B,M}}) where {B,M} = M @@ -35,8 +26,8 @@ function monomial_index(basis::SubBasis, mono::MP.AbstractMonomial) return get(basis, MP.exponents(mono), nothing) end -function explicit_basis_covering(::FullBasis{B}, target::SubBasis{B}) where {B} - return SubBasis{B}(target.monomials) +function explicit_basis_covering(full::FullBasis{B}, target::SubBasis{B}) where {B} + return SA.SubBasis(full, target.keys) end MP.monomial_type(::Type{<:SubBasis{B,M}}) where {B,M} = M @@ -69,10 +60,12 @@ function unsafe_basis( ::Type{B}, monomials::AbstractVector{M}, ) where {B<:AbstractMonomialIndexed,M<:MP.AbstractMonomial} - return SubBasis{B,M,typeof(monomials)}(monomials) + # 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,M}(MP.variables(monomials)), MP.exponents.(monomials)) end -function Base.getindex(::FullBasis{B,M}, monomials::AbstractVector) where {B,M} +function Base.getindex(::FullBasis{B}, monomials::AbstractVector{M}) where {B,M<:MP.AbstractMonomial} return unsafe_basis(B, MP.monomial_vector(monomials)::AbstractVector{M}) end @@ -87,14 +80,6 @@ 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 @@ -225,4 +210,4 @@ function constant_algebra_element(::Type{<:SubBasis{B,M}}, α) where {B,M} [_one_if_type(α)], SubBasis{B}([MP.constant_monomial(M)]), ) -end \ No newline at end of file +end diff --git a/src/chebyshev.jl b/src/chebyshev.jl index aeb35d4..ef19683 100644 --- a/src/chebyshev.jl +++ b/src/chebyshev.jl @@ -24,57 +24,61 @@ 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, +function Base.:*( + a::Polynomial{B}, + b::Polynomial{B}, ) where {B<:AbstractMonomialIndexed} - terms = [MP.term(1 // 1, MP.constant_monomial(a * b))] - vars_a = MP.effective_variables(a) + exps = [constant_monomial_exponents(a.variables)] + coefs = [1 // 1] + vars_a = findall(!iszero, a.exponents) var_state_a = iterate(vars_a) - vars_b = MP.effective_variables(b) + vars_b = findall(!iszero, b.exponents) 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.exponents[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.exponents[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, + B, + exps, + coefs, var_a, - MP.degree(a, var_a), - MP.degree(b, var_b), + a.exponents[var_a], + b.exponents[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 ? + # TODO get rid of this map_keys + @show exps + return SA.map_keys(exp -> Polynomial(a.variables, exp), MA.operate!(SA.canonical, SA.SparseCoefficients(exps, coefs))) end function _add_mul_scalar_vector!(res, ::SubBasis, scalar, vector) 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/monomial.jl b/src/monomial.jl index f44dee4..c33983d 100644 --- a/src/monomial.jl +++ b/src/monomial.jl @@ -38,8 +38,9 @@ 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 Base.:*(a::Polynomial{Monomial}, b::Polynomial{Monomial}) + @assert a.variables == b.variables + return SA.SparseCoefficients((Polynomial(a.variables, a.exponents .+ b.exponents),), (1,)) end SA.coeffs(p::Polynomial{Monomial}, ::FullBasis{Monomial}) = p.monomial @@ -109,16 +110,16 @@ 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) @@ -203,4 +204,4 @@ end function Base.show(io::IO, basis::SubBasis) return show(io, MIME"text/plain"(), basis) -end \ No newline at end of file +end diff --git a/src/orthogonal.jl b/src/orthogonal.jl index be12ca4..d4cf655 100644 --- a/src/orthogonal.jl +++ b/src/orthogonal.jl @@ -117,19 +117,32 @@ 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(::FullBasis{B,M}, exps) where {B,M} + 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 @@ -141,7 +154,6 @@ 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)) end @@ -149,7 +161,7 @@ function explicit_basis_covering( full::FullBasis{B,M}, monos::SubBasis{<:AbstractMonomial,M}, ) where {B<:AbstractMultipleOrthogonal,M} - return SubBasis{B}(_covering(full, monos.monomials)) + return SA.SubBasis(full, _covering(full, monos.keys)) end function _scalar_product_function(::Type{<:AbstractMultipleOrthogonal}, i::Int) end diff --git a/src/polynomial.jl b/src/polynomial.jl index 1425aea..e3d9140 100644 --- a/src/polynomial.jl +++ b/src/polynomial.jl @@ -51,37 +51,44 @@ abstract type AbstractMonomialIndexed end Polynomial of basis `FullBasis{B,M}()` at index `monomial`. """ -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 + 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,20 +100,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.monomial_type(typeof(p))}(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} - M = promote_type(typeof(a.monomial), typeof(b.monomial)) - return algebra_element( - Mul{B,M}()(a.monomial, b.monomial), - FullBasis{B,M}(), - ) + return _algebra_element(MP.monomial(p.variables.variables, p.monomial), B) end function Base.:*(a::Polynomial{B}, b::SA.AlgebraElement) where {B} @@ -118,7 +117,7 @@ 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 @@ -157,19 +156,6 @@ function convert_basis(basis::SA.AbstractBasis, p::SA.AlgebraElement) return SA.AlgebraElement(SA.coeffs(p, basis), algebra(basis)) end -struct Mul{B<:AbstractMonomialIndexed,M} <: SA.MultiplicativeStructure{Polynomial{B,M},M} 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 diff --git a/src/scaled.jl b/src/scaled.jl index fde028e..4a1bce4 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 Base.:*(a::Polynomial{ScaledMonomial}, b::Polynomial{ScaledMonomial}) + @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((Polynomial(a.variables, exp),), (√α,)) end function SA.coeffs(p::Polynomial{ScaledMonomial}, ::FullBasis{Monomial}) diff --git a/src/trash.jl b/src/trash.jl index 2fa91ce..b3abd02 100644 --- a/src/trash.jl +++ b/src/trash.jl @@ -32,7 +32,7 @@ Base.length(basis::SubBasis) = length(basis.monomials) Base.firstindex(basis::SubBasis) = firstindex(basis.monomials) Base.lastindex(basis::SubBasis) = lastindex(basis.monomials) -Base.parent(::SubBasis{B,M}) where {B,M} = FullBasis{B,M}() +#Base.parent(b::SubBasis{B,M}) where {B,M} = FullBasis{B,M}(MP.variables(b)) function Base.getindex(basis::SubBasis, index::Int) return parent(basis)[basis.monomials[index]] diff --git a/src/trigonometric.jl b/src/trigonometric.jl index 7b6ce57..62fd025 100644 --- a/src/trigonometric.jl +++ b/src/trigonometric.jl @@ -21,13 +21,12 @@ _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)) diff --git a/test/chebyshev.jl b/test/chebyshev.jl index 68ebdf1..8d4d846 100644 --- a/test/chebyshev.jl +++ b/test/chebyshev.jl @@ -5,9 +5,11 @@ using DynamicPolynomials @testset "StarAlgebras" begin @polyvar x a = MB.Polynomial{MB.Chebyshev}(x) - b = a * a - @test b.coeffs == MB.sparse_coefficients(1 // 2 + 1 // 2 * x^2) - c = b * b + b = SA.map_keys(exponents, a * a) + @test b == SA.SparseCoefficients([[0], [2]], [1 // 2, 1 // 2]) + full = MB.FullBasis{MB.Chebyshev,typeof(x^1)}([x]) + el = SA.AlgebraElement(b, MB.algebra(full)) + c = el * el @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))) == diff --git a/test/runtests.jl b/test/runtests.jl index 9c7f370..e6e2471 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -42,7 +42,7 @@ 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,M}(x) _test_basis(full_basis) @test sprint(show, MB.algebra(full_basis)) == "Polynomial algebra of $B basis" @@ -237,7 +237,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,typeof(x * y)}(variables(x * y)), SubBasis{MB.Monomial}(monomials(p)), ) coefficient_test(basis, p, coefs; kwargs...) From 2858f067601779bf4c8c889e4fdf1c8ac028c13b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Sat, 14 Jun 2025 09:35:50 +0200 Subject: [PATCH 03/19] Fixes --- src/MultivariateBases.jl | 3 ++- src/bases.jl | 14 +++++++--- src/chebyshev.jl | 35 +++++++++--------------- src/iterator.jl | 57 ---------------------------------------- src/monomial.jl | 11 ++++---- src/mstructures.jl | 11 ++++++++ src/orthogonal.jl | 9 ++++--- src/polynomial.jl | 4 +-- src/scaled.jl | 4 +-- test/chebyshev.jl | 23 +++++++++------- 10 files changed, 64 insertions(+), 107 deletions(-) delete mode 100644 src/iterator.jl create mode 100644 src/mstructures.jl diff --git a/src/MultivariateBases.jl b/src/MultivariateBases.jl index 5d1b564..9d1fe4b 100644 --- a/src/MultivariateBases.jl +++ b/src/MultivariateBases.jl @@ -32,6 +32,7 @@ function MP.polynomial_type(::Type{SA.StarAlgebra{O,S,M}}, ::Type{T}) where {O,S end include("bases.jl") +include("mstructures.jl") include("monomial.jl") include("scaled.jl") @@ -60,7 +61,7 @@ include("quotient.jl") function algebra( basis::Union{QuotientBasis{Polynomial{B,M}},FullBasis{B,M},SubBasis{B,M}}, ) where {B,M} - return SA.StarAlgebra(Variables{B}(MP.variables(basis)), basis) + return SA.StarAlgebra(Variables{B}(MP.variables(basis)), MStruct(basis)) end function MA.promote_operation( diff --git a/src/bases.jl b/src/bases.jl index 807cc42..b305df5 100644 --- a/src/bases.jl +++ b/src/bases.jl @@ -9,6 +9,12 @@ function FullBasis{B,M}(vars) where {B,M} return SA.MappedBasis{Polynomial{B,typeof(vars),eltype(exps)}}(exps, v, MP.exponents) end +FullBasis{B}(vars) where {B} = FullBasis{B,MP.monomial_type(vars)}(vars) + +function FullBasis{B}(p::MP.AbstractPolynomialLike) where {B} + return FullBasis{B}(MP.variables(p)) +end + 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) @@ -84,7 +90,7 @@ 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::SubBasis) = parent(basis) implicit_basis(basis::FullBasis) = basis function implicit(a::SA.AlgebraElement) @@ -148,13 +154,13 @@ _lazy_collect(v::Vector) = v function sparse_coefficients(p::MP.AbstractPolynomial) return SA.SparseCoefficients( - _lazy_collect(MP.monomials(p)), + MP.exponents.(MP.monomials(p)), _lazy_collect(MP.coefficients(p)), ) end function sparse_coefficients(t::MP.AbstractTermLike) - return SA.SparseCoefficients((MP.monomial(t),), (MP.coefficient(t),)) + return SA.SparseCoefficients((MP.exponents(t),), (MP.coefficient(t),)) end function MA.promote_operation( @@ -169,7 +175,7 @@ end function algebra_element(p::MP.AbstractPolynomialLike) return algebra_element( sparse_coefficients(p), - FullBasis{Monomial,MP.monomial_type(p)}(), + FullBasis{Monomial}(p), ) end diff --git a/src/chebyshev.jl b/src/chebyshev.jl index ef19683..66c63fb 100644 --- a/src/chebyshev.jl +++ b/src/chebyshev.jl @@ -35,29 +35,26 @@ function univariate_mul!(::Type{Chebyshev}, exps, coefs, var, a, b) return end -function Base.:*( - a::Polynomial{B}, - b::Polynomial{B}, -) where {B<:AbstractMonomialIndexed} - exps = [constant_monomial_exponents(a.variables)] +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.exponents) + vars_a = findall(!iszero, a) var_state_a = iterate(vars_a) - vars_b = findall(!iszero, b.exponents) + 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(exps) - exps[i] = _increment!(exps[i], b.exponents[var_b], var_b) + 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(exps) - exps[i] = _increment!(exps[i], a.exponents[var_a], var_a) + exps[i] = _increment!(exps[i], a[var_a], var_a) end var_state_a = iterate(vars_a, state_a) else @@ -68,17 +65,15 @@ function Base.:*( exps, coefs, var_a, - a.exponents[var_a], - b.exponents[var_b], + a[var_a], + b[var_b], ) var_state_a = iterate(vars_a, state_a) var_state_b = iterate(vars_b, state_b) end end # FIXME is `canonical` needed ? - # TODO get rid of this map_keys - @show exps - return SA.map_keys(exp -> Polynomial(a.variables, exp), MA.operate!(SA.canonical, SA.SparseCoefficients(exps, coefs))) + return MA.operate!(SA.canonical, SA.SparseCoefficients(exps, coefs)) end function _add_mul_scalar_vector!(res, ::SubBasis, scalar, vector) @@ -123,10 +118,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) + extended = SA.SubBasis(parent(source), sub.keys) ext = 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)), @@ -136,12 +131,8 @@ function SA.coeffs( ) end -function SA.coeffs(cfs, ::FullBasis{Monomial}, target::FullBasis{Chebyshev}) - return SA.coeffs( - SA.values(cfs), - SubBasis{Monomial}(collect(SA.keys(cfs))), - target, - ) +function SA.coeffs(cfs, src::FullBasis{Monomial}, target::FullBasis{Chebyshev}) + return SA.coeffs(SA.values(cfs), SA.SubBasis(src, collect(SA.keys(cfs))), target) end function degree_one_univariate_polynomial(::Type{Chebyshev}, variable) diff --git a/src/iterator.jl b/src/iterator.jl deleted file mode 100644 index e6fba8d..0000000 --- a/src/iterator.jl +++ /dev/null @@ -1,57 +0,0 @@ -import StarAlgebras as SA - -struct ExponentsIterator end -Base.eltype(::Type{ExponentsIterator}) = NTuple{2,Int} -Base.IteratorSize(::Type{ExponentsIterator}) = Base.IsInfinite() - -function Base.iterate(::ExponentsIterator) - z = (0, 0) - return z, (z, 0) -end - -function Base.iterate(::ExponentsIterator, state) - z, deg = state - if iszero(z[2]) - deg += 1 - z = (0, deg) - else - z = (z[1] + 1, z[2] - 1) - end - return z, (z, deg) -end - -function grlex(a::NTuple{2,Int}, b::NTuple{2,Int}) - return isless((sum(a), a), (sum(b), b)) -end - -struct Monomial - exponents::NTuple{2,Int} -end - -Base.one(::Monomial) = Monomial((0, 0)) -Base.:*(a::Monomial, b::Monomial) = Monomial(a.exponents .+ b.exponents) - -monomial(exp) = Monomial(exp) -exponents(mono::Monomial) = mono.exponents - -exps = ExponentsIterator() -SA.comparable(::ExponentsIterator) = grlex - -basis = SA.MappedBasis(exps, monomial, exponents) -mstr = SA.DiracMStructure(basis, *) -object = Monomial((0, 0)) -alg = SA.StarAlgebra(object, mstr) -one(alg) -a = SA.AlgebraElement( - SA.SparseCoefficients( - collect(Iterators.take(exps, 3)), - [2, -1, 3], - ), - alg, -) -c = a * a -import MutableArithmetics as MA -MA.operate!(SA.canonical, c) -@edit MA.operate!(SA.canonical, c) -@edit MA.operate!(SA.canonical, SA.coeffs(c)) -@edit MA.operate!(SA.canonical, SA.coeffs(c), isless) \ No newline at end of file diff --git a/src/monomial.jl b/src/monomial.jl index c33983d..b996b5a 100644 --- a/src/monomial.jl +++ b/src/monomial.jl @@ -38,9 +38,8 @@ function recurrence_eval(::Type{Monomial}, previous, value, degree) return previous[degree] * value end -function Base.:*(a::Polynomial{Monomial}, b::Polynomial{Monomial}) - @assert a.variables == b.variables - return SA.SparseCoefficients((Polynomial(a.variables, a.exponents .+ b.exponents),), (1,)) +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 @@ -170,7 +169,7 @@ 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) +function _show_vector(io::IO, mime::MIME, b, v) print(io, '[') first = true for el in v @@ -178,14 +177,14 @@ function _show_vector(io::IO, mime::MIME, v) print(io, ", ") end first = false - show(io, mime, el) + show(io, mime, MP.monomial(MP.variables(b), 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) + _show_vector(io, mime, basis, basis.keys) print(io, ')') return end diff --git a/src/mstructures.jl b/src/mstructures.jl new file mode 100644 index 0000000..db6e033 --- /dev/null +++ b/src/mstructures.jl @@ -0,0 +1,11 @@ +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 (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 d4cf655..f3b7920 100644 --- a/src/orthogonal.jl +++ b/src/orthogonal.jl @@ -226,13 +226,14 @@ 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(MP.variables(p), p.exponents) return sparse_coefficients( prod( - MP.powers(p.monomial); - init = MP.constant_monomial(M), + 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 e3d9140..58f00d5 100644 --- a/src/polynomial.jl +++ b/src/polynomial.jl @@ -104,8 +104,8 @@ function _algebra_element(p, ::Type{B}) where {B<:AbstractMonomialIndexed} ) end -function algebra_element(p::Polynomial{B,M}) where {B,M} - return _algebra_element(MP.monomial(p.variables.variables, p.monomial), B) +function algebra_element(p::Polynomial{B}) where {B} + return _algebra_element(MP.monomial(MP.variables(p), p.exponents), B) end function Base.:*(a::Polynomial{B}, b::SA.AlgebraElement) where {B} diff --git a/src/scaled.jl b/src/scaled.jl index 4a1bce4..6e794bb 100644 --- a/src/scaled.jl +++ b/src/scaled.jl @@ -33,7 +33,7 @@ Foundations of Computational Mathematics 7.2 (2007): 229-244. """ struct ScaledMonomial <: AbstractMonomial end -function Base.:*(a::Polynomial{ScaledMonomial}, b::Polynomial{ScaledMonomial}) +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( @@ -42,7 +42,7 @@ function Base.:*(a::Polynomial{ScaledMonomial}, b::Polynomial{ScaledMonomial}) ) do i return binomial(exp[i], a.exponents[i]) end - return SA.SparseCoefficients((Polynomial(a.variables, exp),), (√α,)) + return SA.SparseCoefficients((exp,), (1,)) end function SA.coeffs(p::Polynomial{ScaledMonomial}, ::FullBasis{Monomial}) diff --git a/test/chebyshev.jl b/test/chebyshev.jl index 8d4d846..3c2e784 100644 --- a/test/chebyshev.jl +++ b/test/chebyshev.jl @@ -1,19 +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) - b = SA.map_keys(exponents, a * a) - @test b == SA.SparseCoefficients([[0], [2]], [1 // 2, 1 // 2]) - full = MB.FullBasis{MB.Chebyshev,typeof(x^1)}([x]) - el = SA.AlgebraElement(b, MB.algebra(full)) - c = el * el - @test c.coeffs == + full = MB.FullBasis{MB.Chebyshev}(x) + sub = SA.SubBasis(full, [[1]]) + a = MB.algebra_element(MB.sparse_coefficients(1 // 1 * x), full) + @test a == MB.convert_basis(full, MB.algebra_element(x)) + b = a * a + expected = MB.algebra_element(MB.sparse_coefficients((1 // 2) * (1 + x^2)), full) + @test b == expected + c = b * b + @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 From 973d00a3e8aff2e047c72d60cb76030890e0aab5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Sat, 14 Jun 2025 12:19:26 +0200 Subject: [PATCH 04/19] Fixes --- src/MultivariateBases.jl | 6 ++++-- src/arithmetic.jl | 6 +++--- src/bases.jl | 21 ++++++++++----------- src/monomial.jl | 2 +- src/orthogonal.jl | 13 +++++++------ src/polynomial.jl | 13 +++++++------ test/lagrange.jl | 2 +- test/runtests.jl | 27 +++++++++++---------------- 8 files changed, 44 insertions(+), 46 deletions(-) diff --git a/src/MultivariateBases.jl b/src/MultivariateBases.jl index 9d1fe4b..6e3129b 100644 --- a/src/MultivariateBases.jl +++ b/src/MultivariateBases.jl @@ -14,6 +14,8 @@ end Variables{B}(vars) where {B} = Variables{B,typeof(vars)}(vars) +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) @@ -27,8 +29,8 @@ 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{SA.StarAlgebra{O,S,M}}, ::Type{T}) where {O,S,M,T} - return MP.polynomial_type(M, T) +function MP.polynomial_type(::Type{A}, ::Type{T}) where {A<:SA.StarAlgebra,T} + return MP.polynomial_type(MP.monomial_type(A), T) end include("bases.jl") diff --git a/src/arithmetic.jl b/src/arithmetic.jl index 1be0e15..e4e3faa 100644 --- a/src/arithmetic.jl +++ b/src/arithmetic.jl @@ -1,6 +1,6 @@ const _APL = MP.AbstractPolynomialLike # We don't define it for all `AlgebraElement` as this would be type piracy -const _AE = SA.AlgebraElement{<:Variables} +const _AE = SA.AlgebraElement{<:SA.StarAlgebra{<:Variables}} for op in [:+, :-, :*] @eval begin @@ -76,10 +76,10 @@ for op in [:+, :-] 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}(MP.variables(p)), + FullBasis{B}(MP.variables(p)), ) end diff --git a/src/bases.jl b/src/bases.jl index b305df5..dfb8225 100644 --- a/src/bases.jl +++ b/src/bases.jl @@ -2,21 +2,19 @@ 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}} -function FullBasis{B,M}(vars) where {B,M} +function FullBasis{B}(vars) where {B} O = typeof(MP.ordering(first(vars))) # TODO upstream MP.ordering(vars) that handles differently tuples and vectors 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 -FullBasis{B}(vars) where {B} = FullBasis{B,MP.monomial_type(vars)}(vars) - function FullBasis{B}(p::MP.AbstractPolynomialLike) where {B} return FullBasis{B}(MP.variables(p)) end -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} +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 @@ -64,11 +62,11 @@ 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} + 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,M}(MP.variables(monomials)), MP.exponents.(monomials)) + 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} @@ -76,7 +74,7 @@ function Base.getindex(::FullBasis{B}, monomials::AbstractVector{M}) where {B,M< end function SubBasis{B}( - monomials::AbstractVector, + monomials::AbstractVector{<:MP.AbstractTermLike}, ) where {B<:AbstractMonomialIndexed} return unsafe_basis( B, @@ -115,8 +113,8 @@ function MA.promote_operation( return FullBasis{B,M} end -function _explicit_basis(coeffs, ::FullBasis{B}) where {B} - return SubBasis{B}(SA.keys(coeffs)) +function _explicit_basis(coeffs, basis::FullBasis{B}) where {B} + return SA.SubBasis(basis, _lazy_collect(SA.keys(coeffs))) end _explicit_basis(_, basis::SubBasis) = basis @@ -150,6 +148,7 @@ 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) diff --git a/src/monomial.jl b/src/monomial.jl index b996b5a..b217d7a 100644 --- a/src/monomial.jl +++ b/src/monomial.jl @@ -154,7 +154,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), diff --git a/src/orthogonal.jl b/src/orthogonal.jl index f3b7920..8eef147 100644 --- a/src/orthogonal.jl +++ b/src/orthogonal.jl @@ -131,7 +131,7 @@ _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(::FullBasis{B,M}, exps) where {B,M} +function _covering(b::FullBasis{B}, exps) where {B} to_add = collect(exps) m = Set{eltype(exps)}(to_add) while !isempty(to_add) @@ -147,14 +147,14 @@ function _covering(::FullBasis{B,M}, exps) where {B,M} 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} - return SubBasis{BM}(_covering(full, monos.monomials)) + return SA.SubBasis(full, _covering(FullBasis{B}(MP.variables(full)), monos.keys)) end function explicit_basis_covering( @@ -210,8 +210,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 @@ -229,7 +230,7 @@ function SA.coeffs( p::Polynomial{B}, ::FullBasis{Monomial}, ) where {B<:AbstractMultipleOrthogonal} - mono = MP.monomial(MP.variables(p), p.exponents) + mono = MP.monomial(p) return sparse_coefficients( prod( MP.powers(mono); diff --git a/src/polynomial.jl b/src/polynomial.jl index 58f00d5..b3ebc3b 100644 --- a/src/polynomial.jl +++ b/src/polynomial.jl @@ -69,6 +69,8 @@ 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(p.variables, hash(p.monomial, u)) end @@ -105,7 +107,7 @@ function _algebra_element(p, ::Type{B}) where {B<:AbstractMonomialIndexed} end function algebra_element(p::Polynomial{B}) where {B} - return _algebra_element(MP.monomial(MP.variables(p), p.exponents), B) + return _algebra_element(MP.monomial(p), B) end function Base.:*(a::Polynomial{B}, b::SA.AlgebraElement) where {B} @@ -156,14 +158,13 @@ function convert_basis(basis::SA.AbstractBasis, p::SA.AlgebraElement) return SA.AlgebraElement(SA.coeffs(p, basis), algebra(basis)) end -function MP.polynomial(a::SA.AbstractCoefficients) - return MP.polynomial(collect(SA.values(a)), collect(SA.keys(a))) +function _polynomial(b::FullBasis{Monomial}, c::SA.SparseCoefficients) + return MP.polynomial(collect(SA.values(c)), MP.monomial.(getindex.(Ref(b), SA.keys(c)))) end function MP.polynomial(a::SA.AlgebraElement) - return MP.polynomial( - SA.coeffs(a, FullBasis{Monomial,MP.monomial_type(typeof(a))}()), - ) + b = FullBasis{Monomial}(MP.variables(a)) + return _polynomial(b, SA.coeffs(a, b)) end function Base.isapprox( diff --git a/test/lagrange.jl b/test/lagrange.jl index db988e7..9a31dd5 100644 --- a/test/lagrange.jl +++ b/test/lagrange.jl @@ -26,7 +26,7 @@ function _test(B::Type) @test typeof(lag) == MB.explicit_basis_type(typeof(implicit)) a = MB.algebra_element( SA.SparseCoefficients(collect(monos), coeffs), - MB.FullBasis{B,typeof(prod(x))}(), + MB.FullBasis{B}(x), ) @test SA.coeffs(a, lag) ≈ SA.coeffs(coeffs, sub, lag) @polyvar z diff --git a/test/runtests.jl b/test/runtests.jl index e6e2471..af4f206 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -41,8 +41,7 @@ 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}(x) + full_basis = FullBasis{B}(x) _test_basis(full_basis) @test sprint(show, MB.algebra(full_basis)) == "Polynomial algebra of $B basis" @@ -64,8 +63,7 @@ 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 @@ -86,7 +84,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,7 +93,7 @@ function api_test(B::Type{<:MB.AbstractMonomialIndexed}, degree) end end mono = x[1]^2 * x[2]^3 - p = MB.Polynomial{B}(mono) + p = MB.Polynomial(MB.Variables{B}(variables(mono)), exponents(mono)) @test full_basis[p] == mono @test full_basis[mono] == p @test polynomial_type(mono, B == Monomial ? TypeA : TypeB) == @@ -121,7 +119,7 @@ function api_test(B::Type{<:MB.AbstractMonomialIndexed}, degree) _wrap(MB.SA.trim_LaTeX(mime, sprint(show, mime, p.monomial))) * " \$\$" 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,7 +155,7 @@ 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...) @@ -179,7 +177,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 +185,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)] @@ -199,10 +197,7 @@ function orthogonal_test( @test polynomial(basis[i]) == 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])), - ) + basis = explicit_basis_covering(FullBasis{B}(x^2), FullBasis{Monomial}(x^2)[[x^4, x^2, x]]) if even_odd_separated exps = [0, 1, 2, 4] else @@ -237,7 +232,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)}(variables(x * y)), + FullBasis{B}(variables(x * y)), SubBasis{MB.Monomial}(monomials(p)), ) coefficient_test(basis, p, coefs; kwargs...) From f2f02b44c25445fa5696ac188951d5b242e1f392 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Sat, 14 Jun 2025 12:19:34 +0200 Subject: [PATCH 05/19] Remove trash --- src/trash.jl | 49 ------------------------------------------------- 1 file changed, 49 deletions(-) delete mode 100644 src/trash.jl diff --git a/src/trash.jl b/src/trash.jl deleted file mode 100644 index b3abd02..0000000 --- a/src/trash.jl +++ /dev/null @@ -1,49 +0,0 @@ -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 - -# 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) - -#Base.parent(b::SubBasis{B,M}) where {B,M} = FullBasis{B,M}(MP.variables(b)) - -function Base.getindex(basis::SubBasis, index::Int) - return parent(basis)[basis.monomials[index]] -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 - -const MonomialIndexedBasis{B,M} = Union{SubBasis{B,M},FullBasis{B,M}} From de37170f7e017e8f29b3b48d9b3f667451b5228c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Sat, 14 Jun 2025 12:26:16 +0200 Subject: [PATCH 06/19] Fixes --- src/arithmetic.jl | 9 +++++++++ src/monomial.jl | 18 +++++++++++------- src/orthogonal.jl | 8 ++++---- src/polynomial.jl | 11 +---------- 4 files changed, 25 insertions(+), 21 deletions(-) diff --git a/src/arithmetic.jl b/src/arithmetic.jl index e4e3faa..3304c33 100644 --- a/src/arithmetic.jl +++ b/src/arithmetic.jl @@ -2,6 +2,15 @@ const _APL = MP.AbstractPolynomialLike # We don't define it for all `AlgebraElement` as this would be type piracy const _AE = SA.AlgebraElement{<:SA.StarAlgebra{<:Variables}} +function _polynomial(b::FullBasis{Monomial}, c::SA.SparseCoefficients) + return MP.polynomial(collect(SA.values(c)), 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 function MA.promote_operation( diff --git a/src/monomial.jl b/src/monomial.jl index b217d7a..342511f 100644 --- a/src/monomial.jl +++ b/src/monomial.jl @@ -45,14 +45,14 @@ 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,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 MP.polynomial_type(::Type{Polynomial{B,V,E}}, ::Type{T}) where {B,V,E,T} + return _polynomial_type(B, V, T) end function MP.polynomial(f::Function, mb::SubBasis{Monomial}) @@ -136,8 +136,12 @@ 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 diff --git a/src/orthogonal.jl b/src/orthogonal.jl index 8eef147..db84058 100644 --- a/src/orthogonal.jl +++ b/src/orthogonal.jl @@ -108,7 +108,7 @@ function univariate_orthogonal_basis( return univariate_eval!( B, Vector{ - MP.polynomial_type(Polynomial{B,MP.monomial_type(variable)}, Int), + _polynomial_type(B, typeof(MP.variables(variable)), Int), }( undef, degree + 1, @@ -158,9 +158,9 @@ function explicit_basis_covering( end function explicit_basis_covering( - full::FullBasis{B,M}, - monos::SubBasis{<:AbstractMonomial,M}, -) where {B<:AbstractMultipleOrthogonal,M} + full::FullBasis{B,V,E}, + monos::SubBasis{<:AbstractMonomial,V,E}, +) where {B<:AbstractMultipleOrthogonal,V,E} return SA.SubBasis(full, _covering(full, monos.keys)) end diff --git a/src/polynomial.jl b/src/polynomial.jl index b3ebc3b..040c0bd 100644 --- a/src/polynomial.jl +++ b/src/polynomial.jl @@ -102,7 +102,7 @@ end function _algebra_element(p, ::Type{B}) where {B<:AbstractMonomialIndexed} return algebra_element( sparse_coefficients(p), - FullBasis{B,MP.monomial_type(typeof(p))}(MP.variables(p)), + FullBasis{B}(MP.variables(p)), ) end @@ -158,15 +158,6 @@ function convert_basis(basis::SA.AbstractBasis, p::SA.AlgebraElement) return SA.AlgebraElement(SA.coeffs(p, basis), algebra(basis)) end -function _polynomial(b::FullBasis{Monomial}, c::SA.SparseCoefficients) - return MP.polynomial(collect(SA.values(c)), 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 - function Base.isapprox( p::MP.AbstractPolynomialLike, a::SA.AlgebraElement; From 421e54c2018f0613e1f2eadd20c3ce4051dfe85a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Sat, 14 Jun 2025 12:43:11 +0200 Subject: [PATCH 07/19] Fixes --- src/arithmetic.jl | 2 +- src/bases.jl | 14 ++++++++------ src/monomial.jl | 3 ++- 3 files changed, 11 insertions(+), 8 deletions(-) diff --git a/src/arithmetic.jl b/src/arithmetic.jl index 3304c33..bbcf37a 100644 --- a/src/arithmetic.jl +++ b/src/arithmetic.jl @@ -80,7 +80,7 @@ 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 diff --git a/src/bases.jl b/src/bases.jl index dfb8225..a906b55 100644 --- a/src/bases.jl +++ b/src/bases.jl @@ -3,7 +3,7 @@ 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}} function FullBasis{B}(vars) where {B} - O = typeof(MP.ordering(first(vars))) # TODO upstream MP.ordering(vars) that handles differently tuples and vectors + O = typeof(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) @@ -13,6 +13,10 @@ function FullBasis{B}(p::MP.AbstractPolynomialLike) where {B} return FullBasis{B}(MP.variables(p)) end +SA.star(b::MonomialIndexedBasis{B,V,E}, exp::E) where {B,V,E} = b[SA.star(b[exp])] + +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) @@ -193,12 +197,10 @@ function constant_algebra_element_type( 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} +function constant_algebra_element(b::FullBasis, α) return algebra_element( - sparse_coefficients( - MP.polynomial(MP.term(_one_if_type(α), MP.constant_monomial(M))), - ), - FullBasis{B,M}(), + SA.SparseCoefficients((constant_monomial_exponents(b),), (_one_if_type(α),)), + b, ) end diff --git a/src/monomial.jl b/src/monomial.jl index 342511f..8e1f844 100644 --- a/src/monomial.jl +++ b/src/monomial.jl @@ -16,7 +16,8 @@ function explicit_basis_covering( 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 """ From d89db99754e9f005e91c7bfb811c44342cdb1b75 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Sat, 14 Jun 2025 22:24:20 +0200 Subject: [PATCH 08/19] Fixes --- src/MultivariateBases.jl | 25 +++++++++++++++++++------ src/arithmetic.jl | 2 +- src/bases.jl | 33 ++++++++++----------------------- src/monomial.jl | 14 +++++++------- src/mstructures.jl | 4 ++++ src/polynomial.jl | 2 +- test/runtests.jl | 28 ++++++++++++++++------------ 7 files changed, 58 insertions(+), 50 deletions(-) diff --git a/src/MultivariateBases.jl b/src/MultivariateBases.jl index 6e3129b..7e5e500 100644 --- a/src/MultivariateBases.jl +++ b/src/MultivariateBases.jl @@ -14,6 +14,19 @@ end Variables{B}(vars) where {B} = Variables{B,typeof(vars)}(vars) +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 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 + MP.monomial_type(::Type{Variables{B,V}}) where {B,V} = MP.monomial_type(V) constant_monomial_exponents(v::Variables) = map(_ -> 0, v.variables) @@ -30,7 +43,7 @@ function MP.polynomial_type(::Type{<:SA.AlgebraElement{A,T}}) where {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(MP.monomial_type(A), T) + return MP.polynomial_type(MA.promote_operation(SA.basis, A), T) end include("bases.jl") @@ -61,18 +74,18 @@ include("lagrange.jl") include("quotient.jl") function algebra( - basis::Union{QuotientBasis{Polynomial{B,M}},FullBasis{B,M},SubBasis{B,M}}, -) where {B,M} + 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 bbcf37a..4d4d61a 100644 --- a/src/arithmetic.jl +++ b/src/arithmetic.jl @@ -61,7 +61,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), diff --git a/src/bases.jl b/src/bases.jl index a906b55..58cb381 100644 --- a/src/bases.jl +++ b/src/bases.jl @@ -38,7 +38,7 @@ function explicit_basis_covering(full::FullBasis{B}, target::SubBasis{B}) where return SA.SubBasis(full, target.keys) end -MP.monomial_type(::Type{<:SubBasis{B,M}}) where {B,M} = M +MP.monomial_type(::Type{<:SA.SubBasis{T,I,K,B}}) where {T,I,K,B} = MP.monomial_type(B) # The `i`th index of output is the index of occurence of `x[i]` in `y`, # or `0` if it does not occur. @@ -107,15 +107,11 @@ function MA.promote_operation( 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}}} + return SA.AlgebraElement{A,T,SA.SparseCoefficients{M,T,Vector{M},Vector{T},typeof(isless)}} end -function MA.promote_operation( - ::typeof(implicit_basis), - ::Type{<:Union{FullBasis{B,M},SubBasis{B,M}}}, -) where {B,M} - return FullBasis{B,M} -end +MA.promote_operation(::typeof(implicit_basis), B::Type{<:SA.ImplicitBasis}) = B +MA.promote_operation(::typeof(implicit_basis), ::Type{<:SA.SubBasis{T,I,K,B}}) where {T,I,K,B} = B function _explicit_basis(coeffs, basis::FullBasis{B}) where {B} return SA.SubBasis(basis, _lazy_collect(SA.keys(coeffs))) @@ -131,8 +127,8 @@ 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)} +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( @@ -166,15 +162,6 @@ function sparse_coefficients(t::MP.AbstractTermLike) return SA.SparseCoefficients((MP.exponents(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), @@ -192,9 +179,9 @@ _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} +) where {B,V,E,BT<:FullBasis{B,V,E},T} A = MA.promote_operation(algebra, BT) - return SA.AlgebraElement{A,T,SA.SparseCoefficients{M,T,Vector{M},Vector{T}}} + return SA.AlgebraElement{A,T,SA.SparseCoefficients{E,T,Tuple{E},Tuple{T},typeof(isless)}} end function constant_algebra_element(b::FullBasis, α) @@ -212,9 +199,9 @@ function constant_algebra_element_type( return SA.AlgebraElement{A,T,Vector{T}} end -function constant_algebra_element(::Type{<:SubBasis{B,M}}, α) where {B,M} +function constant_algebra_element(basis::SubBasis, α) return algebra_element( [_one_if_type(α)], - SubBasis{B}([MP.constant_monomial(M)]), + SA.SubBasis(parent(basis), [constant_monomial_exponents(parent(basis))]), ) end diff --git a/src/monomial.jl b/src/monomial.jl index 8e1f844..a2adef5 100644 --- a/src/monomial.jl +++ b/src/monomial.jl @@ -174,7 +174,7 @@ 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, b, v) +function _show_vector(io::IO, mime::MIME, v, map = identity) print(io, '[') first = true for el in v @@ -182,30 +182,30 @@ function _show_vector(io::IO, mime::MIME, b, v) print(io, ", ") end first = false - show(io, mime, MP.monomial(MP.variables(b), el)) + 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, basis.keys) + _show_vector(io, mime, basis.keys, Base.Fix1(MP.monomial, MP.variables(b))) print(io, ')') return end -function Base.show(io::IO, mime::MIME"text/plain", basis::SubBasis) +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::SubBasis) +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::SubBasis) +function Base.print(io::IO, basis::Union{Variables,SubBasis}) return show(io, MIME"text/print"(), basis) end -function Base.show(io::IO, basis::SubBasis) +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 index db6e033..8027002 100644 --- a/src/mstructures.jl +++ b/src/mstructures.jl @@ -2,6 +2,10 @@ struct MStruct{B<:AbstractMonomialIndexed,V,E,I,BT<:SA.AbstractBasis{Polynomial{ basis::BT 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 diff --git a/src/polynomial.jl b/src/polynomial.jl index 040c0bd..45e9772 100644 --- a/src/polynomial.jl +++ b/src/polynomial.jl @@ -181,7 +181,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/test/runtests.jl b/test/runtests.jl index af4f206..a7e60b8 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 @@ -35,6 +35,7 @@ Base.:*(::Float64, ::TypeB) = TypeB() Base.:*(::TypeB, ::Float64) = TypeB() Base.:/(::TypeA, ::Float64) = TypeB() Base.:/(::TypeB, ::Float64) = TypeB() +Base.:+(::TypeA, ::TypeA) = TypeB() Base.:+(::TypeB, ::TypeB) = TypeB() Base.:-(::TypeB, ::TypeB) = TypeB() Base.convert(::Type{TypeB}, ::TypeA) = TypeB() @@ -43,8 +44,11 @@ function api_test(B::Type{<:MB.AbstractMonomialIndexed}, degree) @polyvar x[1:2] 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( @@ -94,8 +98,8 @@ function api_test(B::Type{<:MB.AbstractMonomialIndexed}, degree) end mono = x[1]^2 * x[2]^3 p = MB.Polynomial(MB.Variables{B}(variables(mono)), exponents(mono)) - @test full_basis[p] == mono - @test full_basis[mono] == p + @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) @@ -104,19 +108,19 @@ 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(MB.Variables{B}(variables(const_mono)), exponents(const_mono)) From 03620a7fdc9437cd9e6347316d24cd019dd0ef35 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Sat, 14 Jun 2025 22:32:36 +0200 Subject: [PATCH 09/19] Fixes --- src/bases.jl | 10 +++++----- src/monomial.jl | 3 ++- test/runtests.jl | 3 ++- 3 files changed, 9 insertions(+), 7 deletions(-) diff --git a/src/bases.jl b/src/bases.jl index 58cb381..85344a9 100644 --- a/src/bases.jl +++ b/src/bases.jl @@ -102,12 +102,12 @@ 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)) + ::Type{AE}, +) where {AG,T,AE<:SA.AlgebraElement{AG,T}} + BT = MA.promote_operation(implicit_basis, MA.promote_operation(SA.basis, AE)) A = MA.promote_operation(algebra, BT) - M = MP.monomial_type(BT) - return SA.AlgebraElement{A,T,SA.SparseCoefficients{M,T,Vector{M},Vector{T},typeof(isless)}} + E = SA.key_type(BT) + return SA.AlgebraElement{A,T,SA.SparseCoefficients{E,T,Vector{E},Vector{T},typeof(isless)}} end MA.promote_operation(::typeof(implicit_basis), B::Type{<:SA.ImplicitBasis}) = B diff --git a/src/monomial.jl b/src/monomial.jl index a2adef5..f3326da 100644 --- a/src/monomial.jl +++ b/src/monomial.jl @@ -161,8 +161,9 @@ function SA.coeffs( # `JuMP.VariableRef` to `JuMP.AffExpr` return SA.SparseCoefficients(_vec(source.keys), _vec(cfs)) else + T = _promote_coef(_promote_coef(SA.value_type(cfs), B1), B2) res = SA.zero_coeffs( - _promote_coef(_promote_coef(SA.value_type(cfs), B1), B2), + MA.promote_operation(+, T, T), target, ) return SA.coeffs!(res, cfs, source, target) diff --git a/test/runtests.jl b/test/runtests.jl index a7e60b8..4d08576 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -29,13 +29,14 @@ Base.zero(::Type{TypeB}) = TypeB() Base.iszero(::TypeB) = false LinearAlgebra.adjoint(::TypeB) = TypeB() +Base.:+(::TypeA, ::TypeA) = TypeB() +Base.:+(::TypeB, ::TypeA) = TypeB() Base.:*(::Float64, ::TypeA) = TypeB() Base.:*(::TypeA, ::Float64) = TypeB() Base.:*(::Float64, ::TypeB) = TypeB() Base.:*(::TypeB, ::Float64) = TypeB() Base.:/(::TypeA, ::Float64) = TypeB() Base.:/(::TypeB, ::Float64) = TypeB() -Base.:+(::TypeA, ::TypeA) = TypeB() Base.:+(::TypeB, ::TypeB) = TypeB() Base.:-(::TypeB, ::TypeB) = TypeB() Base.convert(::Type{TypeB}, ::TypeA) = TypeB() From 45c2bb5726c019183802d545d07051cff9a43c3d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Sat, 14 Jun 2025 22:42:19 +0200 Subject: [PATCH 10/19] Fixes --- src/monomial.jl | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/src/monomial.jl b/src/monomial.jl index f3326da..ad7600c 100644 --- a/src/monomial.jl +++ b/src/monomial.jl @@ -161,9 +161,8 @@ function SA.coeffs( # `JuMP.VariableRef` to `JuMP.AffExpr` return SA.SparseCoefficients(_vec(source.keys), _vec(cfs)) else - T = _promote_coef(_promote_coef(SA.value_type(cfs), B1), B2) res = SA.zero_coeffs( - MA.promote_operation(+, T, T), + _promote_coef(_promote_coef(SA.value_type(cfs), B1), B2), target, ) return SA.coeffs!(res, cfs, source, target) @@ -190,7 +189,7 @@ 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(b))) + _show_vector(io, mime, basis.keys, Base.Fix1(MP.monomial, MP.variables(basis))) print(io, ')') return end From 5f6688785cae349b026ee6a5adfef3fa7a0747c3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Sat, 14 Jun 2025 22:42:27 +0200 Subject: [PATCH 11/19] Fixes --- test/runtests.jl | 2 -- 1 file changed, 2 deletions(-) diff --git a/test/runtests.jl b/test/runtests.jl index 4d08576..0b745f1 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -29,8 +29,6 @@ Base.zero(::Type{TypeB}) = TypeB() Base.iszero(::TypeB) = false LinearAlgebra.adjoint(::TypeB) = TypeB() -Base.:+(::TypeA, ::TypeA) = TypeB() -Base.:+(::TypeB, ::TypeA) = TypeB() Base.:*(::Float64, ::TypeA) = TypeB() Base.:*(::TypeA, ::Float64) = TypeB() Base.:*(::Float64, ::TypeB) = TypeB() From b47af5539a71c5d61192caebb5a754c1ba355fed Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Sat, 14 Jun 2025 22:42:35 +0200 Subject: [PATCH 12/19] Fix format --- src/MultivariateBases.jl | 12 ++++++-- src/arithmetic.jl | 5 ++- src/bases.jl | 66 ++++++++++++++++++++++++++++++---------- src/chebyshev.jl | 25 +++++++-------- src/hermite.jl | 4 +-- src/monomial.jl | 26 +++++++++++++--- src/mstructures.jl | 19 ++++++++++-- src/orthogonal.jl | 16 +++++----- src/polynomial.jl | 11 +++++-- src/scaled.jl | 5 ++- src/trigonometric.jl | 3 +- test/chebyshev.jl | 3 +- test/runtests.jl | 15 ++++++--- 13 files changed, 148 insertions(+), 62 deletions(-) diff --git a/src/MultivariateBases.jl b/src/MultivariateBases.jl index 7e5e500..790ceaa 100644 --- a/src/MultivariateBases.jl +++ b/src/MultivariateBases.jl @@ -82,10 +82,18 @@ end function MA.promote_operation( ::typeof(algebra), BT::Type{ - <:Union{QuotientBasis{Polynomial{B,V,E}},FullBasis{B,V,E},SubBasis{B,V,E}}, + <:Union{ + QuotientBasis{Polynomial{B,V,E}}, + FullBasis{B,V,E}, + SubBasis{B,V,E}, + }, }, ) where {B,V,E} - return SA.StarAlgebra{Variables{B,V},Polynomial{B,V,E},MStruct{B,V,E,SA.key_type(BT),BT}} + 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 4d4d61a..0ef5fe4 100644 --- a/src/arithmetic.jl +++ b/src/arithmetic.jl @@ -3,7 +3,10 @@ const _APL = MP.AbstractPolynomialLike const _AE = SA.AlgebraElement{<:SA.StarAlgebra{<:Variables}} function _polynomial(b::FullBasis{Monomial}, c::SA.SparseCoefficients) - return MP.polynomial(collect(SA.values(c)), MP.monomial.(getindex.(Ref(b), SA.keys(c)))) + return MP.polynomial( + collect(SA.values(c)), + MP.monomial.(getindex.(Ref(b), SA.keys(c))), + ) end function MP.polynomial(a::SA.AlgebraElement) diff --git a/src/bases.jl b/src/bases.jl index 85344a9..6fe63d2 100644 --- a/src/bases.jl +++ b/src/bases.jl @@ -6,14 +6,20 @@ function FullBasis{B}(vars) where {B} O = typeof(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) + 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 -SA.star(b::MonomialIndexedBasis{B,V,E}, exp::E) where {B,V,E} = b[SA.star(b[exp])] +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) @@ -34,11 +40,16 @@ function monomial_index(basis::SubBasis, mono::MP.AbstractMonomial) return get(basis, MP.exponents(mono), nothing) end -function explicit_basis_covering(full::FullBasis{B}, target::SubBasis{B}) where {B} +function explicit_basis_covering( + full::FullBasis{B}, + target::SubBasis{B}, +) where {B} return SA.SubBasis(full, target.keys) end -MP.monomial_type(::Type{<:SA.SubBasis{T,I,K,B}}) where {T,I,K,B} = MP.monomial_type(B) +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. @@ -70,10 +81,16 @@ function unsafe_basis( ) 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)) + 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} +function Base.getindex( + ::FullBasis{B}, + monomials::AbstractVector{M}, +) where {B,M<:MP.AbstractMonomial} return unsafe_basis(B, MP.monomial_vector(monomials)::AbstractVector{M}) end @@ -104,14 +121,24 @@ 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)) + BT = + MA.promote_operation(implicit_basis, MA.promote_operation(SA.basis, AE)) A = MA.promote_operation(algebra, BT) E = SA.key_type(BT) - return SA.AlgebraElement{A,T,SA.SparseCoefficients{E,T,Vector{E},Vector{T},typeof(isless)}} + return SA.AlgebraElement{ + A, + T, + SA.SparseCoefficients{E,T,Vector{E},Vector{T},typeof(isless)}, + } end MA.promote_operation(::typeof(implicit_basis), B::Type{<:SA.ImplicitBasis}) = B -MA.promote_operation(::typeof(implicit_basis), ::Type{<:SA.SubBasis{T,I,K,B}}) where {T,I,K,B} = 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))) @@ -163,10 +190,7 @@ function sparse_coefficients(t::MP.AbstractTermLike) end function algebra_element(p::MP.AbstractPolynomialLike) - return algebra_element( - sparse_coefficients(p), - FullBasis{Monomial}(p), - ) + return algebra_element(sparse_coefficients(p), FullBasis{Monomial}(p)) end function algebra_element(f::Function, basis::SubBasis) @@ -181,12 +205,19 @@ function constant_algebra_element_type( ::Type{T}, ) where {B,V,E,BT<:FullBasis{B,V,E},T} A = MA.promote_operation(algebra, BT) - return SA.AlgebraElement{A,T,SA.SparseCoefficients{E,T,Tuple{E},Tuple{T},typeof(isless)}} + return SA.AlgebraElement{ + A, + T, + SA.SparseCoefficients{E,T,Tuple{E},Tuple{T},typeof(isless)}, + } end function constant_algebra_element(b::FullBasis, α) return algebra_element( - SA.SparseCoefficients((constant_monomial_exponents(b),), (_one_if_type(α),)), + SA.SparseCoefficients( + (constant_monomial_exponents(b),), + (_one_if_type(α),), + ), b, ) end @@ -202,6 +233,9 @@ end function constant_algebra_element(basis::SubBasis, α) return algebra_element( [_one_if_type(α)], - SA.SubBasis(parent(basis), [constant_monomial_exponents(parent(basis))]), + SA.SubBasis( + parent(basis), + [constant_monomial_exponents(parent(basis))], + ), ) end diff --git a/src/chebyshev.jl b/src/chebyshev.jl index 66c63fb..3e3591e 100644 --- a/src/chebyshev.jl +++ b/src/chebyshev.jl @@ -35,7 +35,11 @@ function univariate_mul!(::Type{Chebyshev}, exps, coefs, var, a, b) return end -function (m::MStruct{B,V,E})(a::E, b::E, ::Type{E}) where {B<:AbstractMonomialIndexed,V,E} +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) @@ -60,14 +64,7 @@ function (m::MStruct{B,V,E})(a::E, b::E, ::Type{E}) where {B<:AbstractMonomialIn else var_a, state_a = var_state_a var_b, state_b = var_state_b - univariate_mul!( - B, - exps, - coefs, - var_a, - a[var_a], - 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 @@ -132,7 +129,11 @@ function SA.coeffs( end function SA.coeffs(cfs, src::FullBasis{Monomial}, target::FullBasis{Chebyshev}) - return SA.coeffs(SA.values(cfs), SA.SubBasis(src, collect(SA.keys(cfs))), target) + return SA.coeffs( + SA.values(cfs), + SA.SubBasis(src, collect(SA.keys(cfs))), + target, + ) end function degree_one_univariate_polynomial(::Type{Chebyshev}, variable) @@ -146,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 @@ -168,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/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/monomial.jl b/src/monomial.jl index ad7600c..a9d5a07 100644 --- a/src/monomial.jl +++ b/src/monomial.jl @@ -52,7 +52,10 @@ function MP.polynomial_type( return _polynomial_type(B, V, T) end -function MP.polynomial_type(::Type{Polynomial{B,V,E}}, ::Type{T}) where {B,V,E,T} +function MP.polynomial_type( + ::Type{Polynomial{B,V,E}}, + ::Type{T}, +) where {B,V,E,T} return _polynomial_type(B, V, T) end @@ -141,7 +144,7 @@ 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} +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 @@ -189,16 +192,29 @@ 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))) + _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}) +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}) +function Base.show( + io::IO, + mime::MIME"text/print", + basis::Union{Variables,SubBasis}, +) return _show(io, mime, basis) end diff --git a/src/mstructures.jl b/src/mstructures.jl index 8027002..d0f3487 100644 --- a/src/mstructures.jl +++ b/src/mstructures.jl @@ -1,8 +1,17 @@ -struct MStruct{B<:AbstractMonomialIndexed,V,E,I,BT<:SA.AbstractBasis{Polynomial{B,V,E},I}} <: SA.MultiplicativeStructure{Polynomial{B,V,E},I} +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 MA.promote_operation(::typeof(SA.basis), ::Type{MStruct{B,V,E,I,BT}}) where {B,V,E,I,BT} +function MA.promote_operation( + ::typeof(SA.basis), + ::Type{MStruct{B,V,E,I,BT}}, +) where {B,V,E,I,BT} return BT end @@ -10,6 +19,10 @@ 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} +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 db84058..d5b44c4 100644 --- a/src/orthogonal.jl +++ b/src/orthogonal.jl @@ -107,9 +107,7 @@ function univariate_orthogonal_basis( ) return univariate_eval!( B, - Vector{ - _polynomial_type(B, typeof(MP.variables(variable)), Int), - }( + Vector{_polynomial_type(B, typeof(MP.variables(variable)), Int)}( undef, degree + 1, ), @@ -147,14 +145,17 @@ function _covering(b::FullBasis{B}, exps) where {B} end end end - return sort!(collect(m), lt = MP.ordering(b)()) + 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} - return SA.SubBasis(full, _covering(FullBasis{B}(MP.variables(full)), monos.keys)) + return SA.SubBasis( + full, + _covering(FullBasis{B}(MP.variables(full)), monos.keys), + ) end function explicit_basis_covering( @@ -232,10 +233,7 @@ function SA.coeffs( ) where {B<:AbstractMultipleOrthogonal} mono = MP.monomial(p) return sparse_coefficients( - prod( - MP.powers(mono); - init = MP.constant_monomial(mono), - ) 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 45e9772..d2c3a79 100644 --- a/src/polynomial.jl +++ b/src/polynomial.jl @@ -76,7 +76,8 @@ function Base.hash(p::Polynomial{B}, u::UInt) where {B} end function Base.isequal(p::Polynomial{B}, q::Polynomial{B}) where {B} - return isequal(p.variables, q.variables) && isequal(p.exponents, q.exponents) + return isequal(p.variables, q.variables) && + isequal(p.exponents, q.exponents) end Base.isone(p::Polynomial) = all(iszero, p.exponents) @@ -119,7 +120,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, MP.monomial(p.variables.variables, p.exponents)))) + print( + io, + SA.trim_LaTeX( + mime, + sprint(show, mime, MP.monomial(p.variables.variables, p.exponents)), + ), + ) if B != Monomial print(io, ")") end diff --git a/src/scaled.jl b/src/scaled.jl index 6e794bb..562a998 100644 --- a/src/scaled.jl +++ b/src/scaled.jl @@ -71,9 +71,8 @@ function Base.promote_rule( end function scaling(m::MP.AbstractMonomial) - return √( - factorial(MP.degree(m)) / prod(factorial, MP.exponents(m); init = 1), - ) + return √(factorial(MP.degree(m)) / + prod(factorial, MP.exponents(m); init = 1),) end unscale_coef(t::MP.AbstractTerm) = MP.coefficient(t) / scaling(MP.monomial(t)) function SA.coeffs( diff --git a/src/trigonometric.jl b/src/trigonometric.jl index 62fd025..10828bd 100644 --- a/src/trigonometric.jl +++ b/src/trigonometric.jl @@ -84,8 +84,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/chebyshev.jl b/test/chebyshev.jl index 3c2e784..1935995 100644 --- a/test/chebyshev.jl +++ b/test/chebyshev.jl @@ -11,7 +11,8 @@ using DynamicPolynomials a = MB.algebra_element(MB.sparse_coefficients(1 // 1 * x), full) @test a == MB.convert_basis(full, MB.algebra_element(x)) b = a * a - expected = MB.algebra_element(MB.sparse_coefficients((1 // 2) * (1 + x^2)), full) + expected = + MB.algebra_element(MB.sparse_coefficients((1 // 2) * (1 + x^2)), full) @test b == expected c = b * b @test SA.coeffs(c) == diff --git a/test/runtests.jl b/test/runtests.jl index 0b745f1..f36489a 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -66,7 +66,8 @@ 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) - poly = MB.Polynomial(MB.Variables{B}(variables(basis)), basis.keys[i]) + poly = + MB.Polynomial(MB.Variables{B}(variables(basis)), basis.keys[i]) @test basis[i] == poly @test basis[poly] == i end @@ -122,7 +123,10 @@ function api_test(B::Type{<:MB.AbstractMonomialIndexed}, degree) _wrap(MB.SA.trim_LaTeX(mime, sprint(show, mime, monomial(p)))) * " \$\$" const_mono = constant_monomial(prod(x)) - const_poly = MB.Polynomial(MB.Variables{B}(variables(const_mono)), exponents(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) @@ -162,7 +166,7 @@ function univ_orthogonal_test( 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 @@ -200,7 +204,10 @@ function orthogonal_test( @test polynomial(basis[i]) == univariate_x[exps[i][1]+1] * univariate_y[exps[i][2]+1] end - basis = explicit_basis_covering(FullBasis{B}(x^2), FullBasis{Monomial}(x^2)[[x^4, x^2, x]]) + basis = explicit_basis_covering( + FullBasis{B}(x^2), + FullBasis{Monomial}(x^2)[[x^4, x^2, x]], + ) if even_odd_separated exps = [0, 1, 2, 4] else From efe2f86d7a9ac7f30fe621a2cd120a6557196498 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Thu, 28 Aug 2025 08:36:35 +0200 Subject: [PATCH 13/19] Bump test compat --- test/Project.toml | 3 +++ 1 file changed, 3 insertions(+) diff --git a/test/Project.toml b/test/Project.toml index 47026ac..40a56d6 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -7,3 +7,6 @@ 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" From 4a9cca5f00abd682df33ca42dfbb1d6ba864062a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Tue, 2 Sep 2025 10:50:57 +0200 Subject: [PATCH 14/19] Implement missing promote_operation --- src/bases.jl | 67 +++++++++++++++++++++++++++++++++++++++------------- 1 file changed, 50 insertions(+), 17 deletions(-) diff --git a/src/bases.jl b/src/bases.jl index 6fe63d2..8185a3d 100644 --- a/src/bases.jl +++ b/src/bases.jl @@ -2,8 +2,13 @@ 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}} +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 = typeof(MP.ordering(vars)) + 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)}}( @@ -117,19 +122,27 @@ function implicit(a::SA.AlgebraElement) 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)) - A = MA.promote_operation(algebra, BT) - E = SA.key_type(BT) - return SA.AlgebraElement{ - A, - T, - SA.SparseCoefficients{E,T,Vector{E},Vector{T},typeof(isless)}, - } + return _algebra_element_type(Vector{T}, BT) end MA.promote_operation(::typeof(implicit_basis), B::Type{<:SA.ImplicitBasis}) = B @@ -197,19 +210,40 @@ 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} +_similar_type(::Type{V}, ::Type{T}) where {V<:AbstractVector,T} = SA.similar_type(V, T) + +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 {B,V,E,BT<:FullBasis{B,V,E},T} - A = MA.promote_operation(algebra, BT) - return SA.AlgebraElement{ - A, - T, - SA.SparseCoefficients{E,T,Tuple{E},Tuple{T},typeof(isless)}, - } +) where {BT<:FullBasis,T} + return _algebra_element_type(Tuple{T}, BT) end function constant_algebra_element(b::FullBasis, α) @@ -226,8 +260,7 @@ 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}} + return _algebra_element_type(Vector{T}, B) end function constant_algebra_element(basis::SubBasis, α) From 9d02df32507cfba7e3e3b30c7c5d7fba0dcb6622 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Tue, 2 Sep 2025 11:59:55 +0200 Subject: [PATCH 15/19] Fix --- src/bases.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/bases.jl b/src/bases.jl index 8185a3d..00f47fe 100644 --- a/src/bases.jl +++ b/src/bases.jl @@ -260,7 +260,8 @@ function constant_algebra_element_type( ::Type{B}, ::Type{T}, ) where {B<:SubBasis,T} - return _algebra_element_type(Vector{T}, B) + A = MA.promote_operation(algebra, B) + return SA.AlgebraElement{A,T,Vector{T}} end function constant_algebra_element(basis::SubBasis, α) From 61f14a220374587e1e6a1a8459b54fc6ae497471 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Thu, 11 Sep 2025 10:39:51 +0200 Subject: [PATCH 16/19] Remove firstindex, lastindex and copy tests --- test/runtests.jl | 3 --- 1 file changed, 3 deletions(-) diff --git a/test/runtests.jl b/test/runtests.jl index f36489a..91f8705 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -73,9 +73,6 @@ function api_test(B::Type{<:MB.AbstractMonomialIndexed}, degree) 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]) From 7ee9712b9a8666e4000464abc0901c0ca02c603c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Thu, 11 Sep 2025 15:27:59 +0200 Subject: [PATCH 17/19] Fixes --- Project.toml | 2 + src/MultivariateBases.jl | 6 ++- src/arithmetic.jl | 7 +++- src/bases.jl | 21 +++++----- src/chebyshev.jl | 2 +- src/lagrange.jl | 4 +- src/monomial.jl | 37 +++++++++++------ src/mstructures.jl | 4 ++ src/polynomial.jl | 10 ++--- src/quotient.jl | 17 ++++---- src/scaled.jl | 24 +++++------ src/trigonometric.jl | 26 ++++++------ test/chebyshev.jl | 1 - test/fixed.jl | 87 ++++++++++++++++++++-------------------- test/hermite.jl | 3 -- test/lagrange.jl | 2 +- test/monomial.jl | 24 +++++------ test/runtests.jl | 1 + test/scaled.jl | 4 +- test/trigonometric.jl | 8 ++-- 20 files changed, 155 insertions(+), 135 deletions(-) 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 790ceaa..2364451 100644 --- a/src/MultivariateBases.jl +++ b/src/MultivariateBases.jl @@ -14,6 +14,10 @@ end Variables{B}(vars) where {B} = Variables{B,typeof(vars)}(vars) +function variable_index(v::Variables, var) + return findfirst(isequal(var), v.variables) +end + 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` @@ -74,7 +78,7 @@ include("lagrange.jl") include("quotient.jl") function algebra( - basis::Union{QuotientBasis{Polynomial{B}},FullBasis{B},SubBasis{B}}, + basis::Union{QuotientBasis{<:Polynomial{B}},FullBasis{B},SubBasis{B}}, ) where {B} return SA.StarAlgebra(Variables{B}(MP.variables(basis)), MStruct(basis)) end diff --git a/src/arithmetic.jl b/src/arithmetic.jl index 0ef5fe4..a8da053 100644 --- a/src/arithmetic.jl +++ b/src/arithmetic.jl @@ -2,10 +2,13 @@ const _APL = MP.AbstractPolynomialLike # We don't define it for all `AlgebraElement` as this would be type piracy 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)), - MP.monomial.(getindex.(Ref(b), SA.keys(c))), + _collect_if_tuple(MP.monomial.(getindex.(Ref(b), SA.keys(c)))), ) end @@ -90,7 +93,7 @@ end function term_element(α, p::Polynomial{B}) where {B} return algebra_element( - sparse_coefficients(MP.term(α, p.monomial)), + SA.SparseCoefficients((p.exponents,), (α,)), FullBasis{B}(MP.variables(p)), ) end diff --git a/src/bases.jl b/src/bases.jl index 00f47fe..e2ced56 100644 --- a/src/bases.jl +++ b/src/bases.jl @@ -2,6 +2,9 @@ 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 @@ -40,11 +43,6 @@ MP.variables(v::Variables) = v.variables MP.variables(basis::FullBasis) = MP.variables(basis.map) MP.variables(basis::SubBasis) = MP.variables(parent(basis)) -function monomial_index(basis::SubBasis, mono::MP.AbstractMonomial) - @assert MP.variables(basis) == MP.variables(mono) - return get(basis, MP.exponents(mono), nothing) -end - function explicit_basis_covering( full::FullBasis{B}, target::SubBasis{B}, @@ -72,11 +70,16 @@ function multi_findsorted(x, y) return I end +import MergeSorted + 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 + @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 diff --git a/src/chebyshev.jl b/src/chebyshev.jl index 3e3591e..7e3dff3 100644 --- a/src/chebyshev.jl +++ b/src/chebyshev.jl @@ -116,7 +116,7 @@ function SA.coeffs( sub = explicit_basis_covering(target, source) # Need to make A square so that it's UpperTriangular extended = SA.SubBasis(parent(source), sub.keys) - ext = SA.coeffs(algebra_element(cfs, source), extended) + ext = _collect_if_tuple(SA.coeffs(algebra_element(cfs, source), extended)) return SA.SparseCoefficients( sub.keys, #transformation_to(sub, extended) \ ext, # Julia v1.6 converts the matrix to the eltype of the `result` which is bad for JuMP 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 a9d5a07..47057e9 100644 --- a/src/monomial.jl +++ b/src/monomial.jl @@ -1,18 +1,19 @@ 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} @@ -59,19 +60,22 @@ function MP.polynomial_type( return _polynomial_type(B, V, T) end +keys_as_monomials(keys, mb::FullBasis) = MP.monomial.(Ref(mb.map.variables), keys) +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}) @@ -125,17 +129,24 @@ end #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 diff --git a/src/mstructures.jl b/src/mstructures.jl index d0f3487..a1b8166 100644 --- a/src/mstructures.jl +++ b/src/mstructures.jl @@ -8,6 +8,10 @@ struct MStruct{ 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}}, diff --git a/src/polynomial.jl b/src/polynomial.jl index d2c3a79..e2cc2ea 100644 --- a/src/polynomial.jl +++ b/src/polynomial.jl @@ -42,14 +42,12 @@ 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,V,E} variables::Variables{B,V} diff --git a/src/quotient.jl b/src/quotient.jl index 54295a8..c2b9448 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,12 @@ 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 +29,9 @@ 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 + 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 562a998..3d51076 100644 --- a/src/scaled.jl +++ b/src/scaled.jl @@ -58,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 @@ -70,9 +70,9 @@ 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 +87,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,11 +105,11 @@ function SA.coeffs!( ) MA.operate!(zero, res) for (k, v) in SA.nonzero_pairs(cfs) - mono = source[k].monomial + exp = source[k].exponents SA.unsafe_push!( res, - target[Polynomial{Monomial}(mono)], - v * scaling(mono), + exponents_index(target, exp), + v * scaling(exp), ) end MA.operate!(SA.canonical, res) @@ -124,11 +124,11 @@ function SA.coeffs!( ) MA.operate!(zero, res) for (k, v) in SA.nonzero_pairs(cfs) - mono = source[k].monomial + exp = source[k].exponents SA.unsafe_push!( res, - target[Polynomial{ScaledMonomial}(mono)], - v / scaling(mono), + exponents_index(target, exp), + v / scaling(exp), ) end MA.operate!(SA.canonical, res) diff --git a/src/trigonometric.jl b/src/trigonometric.jl index 10828bd..f0bd8ea 100644 --- a/src/trigonometric.jl +++ b/src/trigonometric.jl @@ -29,16 +29,15 @@ function univariate_mul!(::Type{Trigonometric}, exps, coefs, var, a, b) 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) @@ -47,8 +46,8 @@ function univariate_mul!(::Type{Trigonometric}, exps, coefs, 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) @@ -57,13 +56,12 @@ function univariate_mul!(::Type{Trigonometric}, exps, coefs, 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 diff --git a/test/chebyshev.jl b/test/chebyshev.jl index 1935995..6126f58 100644 --- a/test/chebyshev.jl +++ b/test/chebyshev.jl @@ -7,7 +7,6 @@ using DynamicPolynomials @testset "StarAlgebras" begin @polyvar x full = MB.FullBasis{MB.Chebyshev}(x) - sub = SA.SubBasis(full, [[1]]) a = MB.algebra_element(MB.sparse_coefficients(1 // 1 * x), full) @test a == MB.convert_basis(full, MB.algebra_element(x)) b = a * a 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 9a31dd5..56857d1 100644 --- a/test/lagrange.jl +++ b/test/lagrange.jl @@ -25,7 +25,7 @@ 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), + SA.SparseCoefficients(collect(exponents.(monos)), coeffs), MB.FullBasis{B}(x), ) @test SA.coeffs(a, lag) ≈ SA.coeffs(coeffs, sub, lag) 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 91f8705..c5ad1f3 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -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() diff --git a/test/scaled.jl b/test/scaled.jl index 4975b21..c552797 100644 --- a/test/scaled.jl +++ b/test/scaled.jl @@ -24,9 +24,9 @@ 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)}() + full = MB.FullBasis{MB.ScaledMonomial}(variables(p)) @test SA.coeffs(a, full) == - SA.SparseCoefficients([y^2, x * y, x^2], [9, 4 / √2, 1]) + 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..376f6e1 100644 --- a/test/trigonometric.jl +++ b/test/trigonometric.jl @@ -4,12 +4,14 @@ 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 From e3e552346d354c6ab8eb51eb8ca9065e53c921a7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Tue, 14 Oct 2025 09:12:22 +0200 Subject: [PATCH 18/19] Fix format --- src/bases.jl | 19 ++++++++++++------- src/monomial.jl | 4 +++- src/quotient.jl | 8 ++++++-- src/scaled.jl | 15 +++------------ test/scaled.jl | 6 ++++-- test/trigonometric.jl | 5 ++++- 6 files changed, 32 insertions(+), 25 deletions(-) diff --git a/src/bases.jl b/src/bases.jl index e2ced56..aacafa2 100644 --- a/src/bases.jl +++ b/src/bases.jl @@ -76,7 +76,13 @@ 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))) + 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 @@ -216,11 +222,13 @@ 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} -_similar_type(::Type{V}, ::Type{T}) where {V<:AbstractVector,T} = SA.similar_type(V, 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}} + P::Type{<:MP.AbstractPolynomialLike{T}}, ) where {T} V = MA.promote_operation(MP.variables, P) E = _similar_type(V, Int) @@ -230,10 +238,7 @@ function MA.promote_operation( P, E, MP.ExponentsIterator{O,Nothing,E}, - Variables{ - Monomial, - V, - }, + Variables{Monomial,V}, typeof(MP.exponents), } return _algebra_element_type(Vector{T}, B) diff --git a/src/monomial.jl b/src/monomial.jl index 47057e9..d88deba 100644 --- a/src/monomial.jl +++ b/src/monomial.jl @@ -60,7 +60,9 @@ function MP.polynomial_type( return _polynomial_type(B, V, T) end -keys_as_monomials(keys, mb::FullBasis) = MP.monomial.(Ref(mb.map.variables), keys) +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}) diff --git a/src/quotient.jl b/src/quotient.jl index c2b9448..951bc78 100644 --- a/src/quotient.jl +++ b/src/quotient.jl @@ -17,7 +17,10 @@ function MP.coefficients(p, basis::QuotientBasis) end function SA.coeffs(coeffs, b::FullBasis{Monomial}, basis::QuotientBasis) - return MP.coefficients(MP.polynomial(values(coeffs), keys_as_monomials(keys(coeffs), b)), basis) + return MP.coefficients( + MP.polynomial(values(coeffs), keys_as_monomials(keys(coeffs), b)), + basis, + ) end function SA.coeffs(coeffs, sub::SubBasis{Monomial}, basis::QuotientBasis) @@ -31,7 +34,8 @@ function SA.adjoint_coeffs( ) return map(src) do poly return sum(MP.terms(rem(MP.polynomial(poly), dest.divisor))) do t - MP.coefficient(t) * coeffs[SA.key_index(dest.basis, MP.exponents(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 3d51076..ce05766 100644 --- a/src/scaled.jl +++ b/src/scaled.jl @@ -71,8 +71,7 @@ function Base.promote_rule( end function scaling(exp) - return √(factorial(sum(exp)) / - prod(factorial, exp; init = 1),) + 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( @@ -106,11 +105,7 @@ function SA.coeffs!( MA.operate!(zero, res) for (k, v) in SA.nonzero_pairs(cfs) exp = source[k].exponents - SA.unsafe_push!( - res, - exponents_index(target, exp), - v * scaling(exp), - ) + SA.unsafe_push!(res, exponents_index(target, exp), v * scaling(exp)) end MA.operate!(SA.canonical, res) return res @@ -125,11 +120,7 @@ function SA.coeffs!( MA.operate!(zero, res) for (k, v) in SA.nonzero_pairs(cfs) exp = source[k].exponents - SA.unsafe_push!( - res, - exponents_index(target, exp), - v / scaling(exp), - ) + SA.unsafe_push!(res, exponents_index(target, exp), v / scaling(exp)) end MA.operate!(SA.canonical, res) return res diff --git a/test/scaled.jl b/test/scaled.jl index c552797..bf45b2d 100644 --- a/test/scaled.jl +++ b/test/scaled.jl @@ -25,8 +25,10 @@ end a = MB.algebra_element(p) @test SA.coeffs(a, basis) == [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 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 376f6e1..d209e7a 100644 --- a/test/trigonometric.jl +++ b/test/trigonometric.jl @@ -11,7 +11,10 @@ using DynamicPolynomials c = b * b @test c.coeffs == MB.sparse_coefficients(3 // 8 + 1 // 2 * x^3 + 1 // 8 * x^7) - d = MB.algebra_element(MB.sparse_coefficients(1//1 * constant_monomial(typeof(x))), full) + 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 From 01a00872bbae30ff53667b59e84e966930fe1a86 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Tue, 14 Oct 2025 22:57:20 +0200 Subject: [PATCH 19/19] Update formatter to v2 --- .github/workflows/format_check.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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)