Fixes for CondensedMatterSOS (#255)

* Fixes for CondensedMatterSOS

* promote_operation_constant

* Don't assume commutativity

Author: blegat

.github/workflows/ci.yml

@@ -46,17 +46,17 @@ jobs:
         shell: julia --project=@. {0}
         run: |
           using Pkg
-          Pkg.add(PackageSpec(name="TypedPolynomials", rev="bl/rename"))
+          Pkg.add(PackageSpec(name="TypedPolynomials", rev="master"))
       - name: DP#master
         shell: julia --project=@. {0}
         run: |
           using Pkg
-          Pkg.add(PackageSpec(name="DynamicPolynomials", rev="bl/rename"))
+          Pkg.add(PackageSpec(name="DynamicPolynomials", rev="bl/coefficient_type"))
       - uses: julia-actions/julia-buildpkg@v1
       - uses: julia-actions/julia-runtest@v1
         with:
           depwarn: error
       - uses: julia-actions/julia-processcoverage@v1
-      - uses: codecov/codecov-action@v1
+      - uses: codecov/codecov-action@v3
         with:
           file: lcov.info

docs/src/types.md

@@ -31,17 +31,19 @@ isconstant
 powers
 constant_monomial
 map_exponents
+multiplication_preserves_monomial_order
 ```
 
 ## Terms
 
 ```@docs
 AbstractTermLike
 AbstractTerm
+Term
 term
 term_type
 coefficient
-coefficienttype
+coefficient_type
 monomial
 constant_term
 zero_term
@@ -52,6 +54,7 @@ zero_term
 ```@docs
 AbstractPolynomialLike
 AbstractPolynomial
+Polynomial
 polynomial
 polynomial_type
 terms

src/default_polynomial.jl

@@ -1,3 +1,11 @@
+"""
+    struct Polynomial{CoeffType,T<:AbstractTerm{CoeffType},V<:AbstractVector{T}} <: AbstractPolynomial{CoeffType}
+        terms::V
+    end
+
+Representation of a multivariate polynomial as a vector of nonzero terms sorted
+in ascending monomial order.
+"""
 struct Polynomial{CoeffType,T<:AbstractTerm{CoeffType},V<:AbstractVector{T}} <:
        AbstractPolynomial{CoeffType}
     terms::V
@@ -7,8 +15,16 @@ struct Polynomial{CoeffType,T<:AbstractTerm{CoeffType},V<:AbstractVector{T}} <:
     end
 end
 
+# It shouldn't be needed to call this. A polynomial should always be in
+# canonical format so it is only called inside an internal call.
+function _canonicalize!(p::Polynomial)
+    sort!(p.terms)
+    uniqterms!(p.terms)
+    return
+end
+
 function coefficients(p::Polynomial)
-    return LazyMap{coefficienttype(p)}(coefficient, terms(p))
+    return LazyMap{coefficient_type(p)}(coefficient, terms(p))
 end
 function monomials(p::Polynomial)
     return LazyMap{monomial_type(p)}(monomial, terms(p))
@@ -43,7 +59,7 @@ function Base.convert(
     return convert(Polynomial{C,TT,Vector{TT}}, p)
 end
 function Base.convert(::Type{Polynomial}, p::AbstractPolynomialLike)
-    return convert(Polynomial{coefficienttype(p)}, p)
+    return convert(Polynomial{coefficient_type(p)}, p)
 end
 
 _change_eltype(::Type{<:Vector}, ::Type{T}) where {T} = Vector{T}
@@ -141,6 +157,9 @@ function MA.operate_to!(
             # TODO could use MA.mul_to!! for indices that were presents in `result` before the `resize!`.
             result.terms[i] = p.terms[i] * t
         end
+        if !multiplication_preserves_monomial_order(typeof(result))
+            _canonicalize!(result)
+        end
         return result
     end
 end
@@ -158,6 +177,9 @@ function MA.operate_to!(
             # TODO could use MA.mul_to!! for indices that were presents in `result` before the `resize!`.
             result.terms[i] = t * p.terms[i]
         end
+        if !multiplication_preserves_monomial_order(typeof(result))
+            _canonicalize!(result)
+        end
         return result
     end
 end
@@ -268,7 +290,7 @@ function __polynomial_merge!(
             if tp isa Int && j isa Int
                 tp = get1(tp)
             end
-            grlex(*(monomials(q)[j], monomial.(t)...), monomial(tp))
+            grlex(monomial(*(monomials(q)[j], monomial.(t)...)), monomial(tp))
         end
     end
     get2 = let t = t, q = q
@@ -338,12 +360,28 @@ function __polynomial_merge!(
     q::Polynomial,
     buffer = nothing,
 ) where {T,TT}
+    if !multiplication_preserves_monomial_order(typeof(p))
+        # This function assumes that the monomials of `t * q` will preserve the
+        # order of the terms of `q` so if it's not the case, we fallback to
+        # another one.
+        __polynomial_merge!(
+            MA.add_sub_op(op),
+            p,
+            get1,
+            set,
+            push,
+            resize,
+            keep,
+            t * q,
+        )
+        return
+    end
     compare_monomials = let t = t, get1 = get1, q = q
         (tp, j) -> begin
             if tp isa Int && j isa Int
                 tp = get1(tp)
             end
-            grlex(monomial(t) * monomials(q)[j], monomial(tp))
+            grlex(monomial(monomial(t) * monomials(q)[j]), monomial(tp))
         end
     end
     get2 = let t = t, q = q
@@ -550,8 +588,7 @@ function MA.operate_to!(
 )
     empty!(output.terms)
     mul_to_terms!(output.terms, p, q)
-    sort!(output.terms, lt = (<))
-    uniqterms!(output.terms)
+    _canonicalize!(output)
     return output
 end
 function MA.operate!(::typeof(*), p::Polynomial, q::Polynomial)
@@ -567,6 +604,9 @@ function MA.operate!(::typeof(*), p::Polynomial, t::AbstractTermLike)
     for i in eachindex(p.terms)
         p.terms[i] = MA.operate!!(*, p.terms[i], t)
     end
+    if !multiplication_preserves_monomial_order(typeof(p))
+        _canonicalize!(p)
+    end
     return p
 end
 function map_exponents!(f, p::Polynomial, m::AbstractMonomialLike)

src/default_term.jl

@@ -1,3 +1,21 @@
+"""
+    struct Term{CoeffType,M<:AbstractMonomial} <: AbstractTerm{CoeffType}
+        coefficient::CoeffType
+        monomial::M
+    end
+
+A representation of the multiplication between a `coefficient` and a `monomial`.
+
+!!! note
+    The `coefficient` does not need to be a `Number`. It can be for instance a
+    multivariate polynomial. When computing a multivariate `gcd`, it is
+    actually reformulated as a univariate `gcd` in one of the variable with
+    coefficients being multivariate polynomials in the other variables.
+    To create such a term, use [`term`](@ref) instead of `*`.
+    For instance, if `p` is a polynomial and `m` is a monomial, `p * m` will
+    multiply each term of `p` with `m` but `term(p, m)` will create a term
+    with `p` as coefficient and `m` as monomial.
+"""
 struct Term{CoeffType,M<:AbstractMonomial} <: AbstractTerm{CoeffType}
     coefficient::CoeffType
     monomial::M

src/polynomial.jl

@@ -1,6 +1,6 @@
 export polynomial,
     polynomial!, polynomial_type, terms, nterms, coefficients, monomials
-export coefficienttype, monomial_type
+export coefficient_type, monomial_type
 export mindegree, maxdegree, extdegree, effective_variables
 export leading_term, leading_coefficient, leading_monomial
 export remove_leading_term, remove_monomials, monic
@@ -600,5 +600,28 @@ function deg_num_leading_terms(p::AbstractPolynomialLike, var)
     return deg, num
 end
 
+"""
+    multiplication_preserves_monomial_order(P::Type{<:AbstractPolynomialLike})
+
+Returns a `Bool` indicating whether the order is preserved in the
+multiplication of monomials of type `monomial_type(P)`.
+That is, if `a < b` then `a * c < b * c` for any monomial `c`.
+This returns `true` by default.
+This is used by [`Polynomial`](@ref) so a monomial type for which the
+multiplication does not preserve the monomial order can still be used with
+[`Polynomial`](@ref) if it implements a method for this function that
+returns `false`.
+"""
+function multiplication_preserves_monomial_order(
+    ::Type{P},
+) where {P<:AbstractPolynomialLike}
+    return multiplication_preserves_monomial_order(monomial_type(P))
+end
+function multiplication_preserves_monomial_order(
+    ::Type{M},
+) where {M<:AbstractMonomial}
+    return true
+end
+
 Base.ndims(::Union{Type{<:AbstractPolynomialLike},AbstractPolynomialLike}) = 0
 Base.broadcastable(p::AbstractPolynomialLike) = Ref(p)

src/promote.jl

@@ -239,32 +239,64 @@ function MA.promote_operation(
     TT::Type{<:AbstractTermLike{S}},
     ST::Type{<:AbstractTermLike{T}},
 ) where {S,T}
+    UT = MA.promote_operation(*, monomial_type(TT), monomial_type(ST))
     U = MA.promote_operation(*, S, T)
-    return term_type(promote_type(monomial_type(TT), monomial_type(ST)), U)
+    return promote_operation_constant(*, U, UT)
 end
 function MA.promote_operation(
     ::typeof(*),
     PT::Type{<:APL{S}},
     QT::Type{<:APL{T}},
 ) where {S,T}
-    ST = MA.promote_operation(*, S, T)
-    U = MA.promote_operation(+, ST, ST)
-    return polynomial_type(
-        promote_type(monomial_type(PT), monomial_type(QT)),
-        U,
-    )
+    UP = MA.promote_operation(*, monomial_type(PT), monomial_type(QT))
+    U = MA.promote_sum_mul(S, T)
+    return polynomial_type(promote_operation_constant(*, U, UP))
 end
-function MA.promote_operation(
+
+function promote_operation_constant(
+    ::typeof(*),
+    ::Type{T},
+    ::Type{M},
+) where {T,M<:AbstractMonomialLike}
+    return term_type(M, T)
+end
+
+function promote_operation_constant(
+    ::typeof(*),
+    ::Type{M},
+    ::Type{T},
+) where {T,M<:AbstractMonomialLike}
+    return term_type(M, T)
+end
+
+function promote_operation_constant(
     ::typeof(*),
     ::Type{T},
     ::Type{P},
 ) where {T,U,P<:APL{U}}
     return similar_type(P, MA.promote_operation(*, T, U))
 end
-function MA.promote_operation(
+
+function promote_operation_constant(
     ::typeof(*),
     ::Type{P},
     ::Type{T},
 ) where {T,U,P<:APL{U}}
     return similar_type(P, MA.promote_operation(*, U, T))
 end
+
+function MA.promote_operation(
+    ::typeof(*),
+    ::Type{T},
+    ::Type{P},
+) where {T,P<:APL}
+    return promote_operation_constant(*, T, P)
+end
+
+function MA.promote_operation(
+    ::typeof(*),
+    ::Type{P},
+    ::Type{T},
+) where {T,P<:APL}
+    return promote_operation_constant(*, P, T)
+end

src/substitution.jl

@@ -98,7 +98,7 @@ function MA.promote_operation(
     args::Vararg{Type,N},
 ) where {T<:AbstractTerm,N}
     M = MA.promote_operation(substitute, Subs, monomial_type(T), args...)
-    U = coefficienttype(T)
+    U = coefficient_type(T)
     return MA.promote_operation(*, U, M)
 end
 

src/term.jl

@@ -124,26 +124,26 @@ function coefficient(f::APL, m::AbstractMonomialLike, vars)
 end
 
 """
-    coefficienttype(p::AbstractPolynomialLike)
+    coefficient_type(p::AbstractPolynomialLike)
 
 Returns the coefficient type of `p`.
 
-    coefficienttype(::Type{PT}) where PT
+    coefficient_type(::Type{PT}) where PT
 
 Returns the coefficient type of a polynomial of type `PT`.
 
 ### Examples
 
-Calling `coefficienttype` on ``(4//5)x^2y`` should return `Rational{Int}`,
-calling `coefficienttype` on ``1.0x^2y + 2.0x`` should return `Float64` and
-calling `coefficienttype` on ``xy`` should return `Int`.
+Calling `coefficient_type` on ``(4//5)x^2y`` should return `Rational{Int}`,
+calling `coefficient_type` on ``1.0x^2y + 2.0x`` should return `Float64` and
+calling `coefficient_type` on ``xy`` should return `Int`.
 """
-function coefficienttype(
+function coefficient_type(
     ::Union{PT,Type{PT},AbstractVector{PT},Type{<:AbstractVector{PT}}},
 ) where {T,PT<:APL{T}}
     return T
 end
-#coefficienttype(::{T, Type{T}}) where {T} = T
+#coefficient_type(::{T, Type{T}}) where {T} = T
 
 """
     monomial(t::AbstractTermLike)

test/division.jl

@@ -190,7 +190,7 @@ function mult_test(expected, a::Number, b, algo)
         a,
         b,
         algo,
-        promote_type(typeof(a), MP.coefficienttype(b)),
+        promote_type(typeof(a), MP.coefficient_type(b)),
     )
 end
 function mult_test(expected, a, b::Number, algo)
@@ -199,7 +199,7 @@ function mult_test(expected, a, b::Number, algo)
         a,
         b,
         algo,
-        promote_type(MP.coefficienttype(a), typeof(b)),
+        promote_type(MP.coefficient_type(a), typeof(b)),
     )
 end
 function sym_test(a, b, g, algo)

test/utils.jl

@@ -56,7 +56,7 @@ end
 Base.copy(p::CustomTerms) = CustomTerms(copy(p.p))
 
 function _typetests(x, ::Type{T}) where {T}
-    @test (@inferred coefficienttype(x)) == Int
+    @test (@inferred coefficient_type(x)) == Int
 
     @test (@inferred monomial_type(x)) <: AbstractMonomial