Merge branch 'master' of https://github.com/JuliaAlgebra/Multivariate…

…Polynomials.jl

CITATION.bib

@@ -1,12 +1,9 @@
-@software{legat2021multivariatepolynomials,
-  author       = {Benoît Legat and
-                  Sascha Timme and
-                  Robin Deits},
-  title        = {JuliaAlgebra/MultivariatePolynomials.jl: v0.3.18},
+@software{legat2023multivariate,
+  author       = {Legat, Beno{\^\i}t},
+  title        = {Multivariate polynomials in Julia},
   month        = jul,
-  year         = 2021,
-  publisher    = {Zenodo},
-  version      = {v0.3.18},
-  doi          = {10.5281/zenodo.5083796},
-  url          = {https://doi.org/10.5281/zenodo.5083796}
+  year         = 2023,
+  booktitle    = {JuliaCon},
+  year         = {2022},
+  url          = {https://pretalx.com/juliacon-2022/talk/TRFSJY/},
 }

Project.toml

@@ -2,7 +2,7 @@ name = "MultivariatePolynomials"
 uuid = "102ac46a-7ee4-5c85-9060-abc95bfdeaa3"
 license = "MIT"
 repo = "https://github.com/JuliaAlgebra/MultivariatePolynomials.jl"
-version = "0.5.0"
+version = "0.5.3"
 
 [deps]
 ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"

README.md

@@ -17,6 +17,11 @@ The interface contains functions for accessing the coefficients, monomials, term
 - [**STABLE**][docs-stable-url] — **most recently tagged version of the documentation.**
 - [**LATEST**][docs-latest-url] — *in-development version of the documentation.*
 
+## Citing
+
+Please cite the [JuliaCon 2022 presentation](https://pretalx.com/juliacon-2022/talk/TRFSJY/) [[Slides](https://drive.google.com/file/d/1q9UT5rpcmm0GdRmWm7llhVpbOx7OxGnZ/view)].
+See [CITATION.bib](https://github.com/JuliaPolyhedra/MultivariatePolynomials.jl/blob/master/CITATION.bib) for the BibTeX.
+
 ## Examples
 
 Below is a simple usage example

docs/src/types.md

@@ -81,6 +81,8 @@ leading_coefficient
 leading_monomial
 remove_leading_term
 remove_monomials
+filter_terms
+OfDegree
 monic
 map_coefficients
 map_coefficients!

src/default_term.jl

@@ -108,6 +108,14 @@ function MA.operate!(::typeof(*), t1::Term, t2::AbstractTermLike)
     return t1
 end
 
+# Needed to resolve ambiguity with
+# MA.operate!(::typeof(*), ::AbstractMonomial, ::AbstractMonomialLike)
+function MA.operate!(::typeof(*), t1::Term, t2::AbstractMonomialLike)
+    MA.operate!(*, t1.coefficient, coefficient(t2))
+    MA.operate!(*, t1.monomial, monomial(t2))
+    return t1
+end
+
 function MA.operate!(::typeof(one), t::Term)
     MA.operate!(one, t.coefficient)
     MA.operate!(constant_monomial, t.monomial)

src/lazy_iterators.jl

@@ -32,7 +32,7 @@ Iterator over the elements of `data` mapped by `f`. This is similar to
 `Base.Generator(f, data)` except that the `eltype` of a `LazyMap` is given at
 construction while the `eltype` of `Base.Generator(f, data)` is `Any`.
 """
-struct LazyMap{T,VT,F} <: AbstractVector{T}
+struct LazyMap{T,VT<:AbstractVector,F} <: AbstractVector{T}
     f::F
     data::VT
 end
@@ -58,7 +58,10 @@ Base.IteratorSize(it::LazyMap) = Base.IteratorSize(it.data)
 
 Base.eltype(::LazyMap{T}) where {T} = T
 
-Base.getindex(it::LazyMap, i) = it.f(getindex(it.data, i))
+Base.getindex(it::LazyMap, i::Integer) = it.f(getindex(it.data, i))
+function Base.getindex(it::LazyMap{T}, I::AbstractVector) where {T}
+    return LazyMap{T}(it.f, getindex(it.data, I))
+end
 
 Base.eachindex(it::LazyMap) = Base.eachindex(it.data)
 Base.lastindex(it::LazyMap) = Base.lastindex(it.data)

src/monomial_vector.jl

@@ -12,7 +12,7 @@ empty_monomial_vector(p) = monomial_type(p)[]
 """
     monomial_vector(X::AbstractVector{MT}) where {MT<:AbstractMonomialLike}
 
-Returns the vector of monomials `X` in decreasing order and without any duplicates.
+Returns the vector of monomials `X` in increasing order and without any duplicates.
 
 ### Examples
 
@@ -31,7 +31,7 @@ end
 """
     monomial_vector(a, X::AbstractVector{MT}) where {MT<:AbstractMonomialLike}
 
-Returns `b, Y` where `Y` is the vector of monomials of `X` in decreasing order
+Returns `b, Y` where `Y` is the vector of monomials of `X` in increasing order
 and without any duplicates and `b` is the vector of corresponding coefficients
 in `a`, where coefficients of duplicate entries are summed together.
 
@@ -66,7 +66,7 @@ end
 Returns the return type of `monomial_vector`.
 """
 function monomial_vector_type(
-    X::Union{AbstractVector{PT},Type{<:AbstractVector{PT}}},
+    ::Union{AbstractVector{PT},Type{<:AbstractVector{PT}}},
 ) where {PT<:_APL}
     return monomial_vector_type(PT)
 end
@@ -102,10 +102,20 @@ sort_monomial_vector(X::Tuple) = sort_monomial_vector(vec(X))
 """
     merge_monomial_vectors{MT<:AbstractMonomialLike, MVT<:AbstractVector{MT}}(X::AbstractVector{MVT}}
 
-Returns the vector of monomials in the entries of `X` in decreasing order and without any duplicates, i.e. `monomial_vector(vcat(X...))`
+Returns the vector of monomials in the entries of `X` in increasing order and without any duplicates, i.e. `monomial_vector(vcat(X...))`
 
 ### Examples
 
 Calling `merge_monomial_vectors` on ``[[xy, x, xy], [x^2y, x]]`` should return ``[x^2y, xy, x]``.
 """
 merge_monomial_vectors(X) = monomial_vector(reduce(vcat, X))
+
+function error_for_negative_degree(deg)
+    if deg < 0
+        throw(
+            ArgumentError(
+                "The degree should be a nonnegative number but the provided degree `$deg` is negative.",
+            ),
+        )
+    end
+end

src/polynomial.jl

@@ -272,23 +272,29 @@ end
 
 #$(SIGNATURES)
 """
-    mindegree(p::Union{AbstractPolynomialLike, AbstractVector{<:AbstractTermLike}})
+    mindegree(p::Union{AbstractPolynomialLike, AbstractVector{<:AbstractPolynomialLike}})
 
 Returns the minimal total degree of the monomials of `p`, i.e. `minimum(degree, terms(p))`.
 
-    mindegree(p::Union{AbstractPolynomialLike, AbstractVector{<:AbstractTermLike}}, v::AbstractVariable)
+    mindegree(p::Union{AbstractPolynomialLike, AbstractVector{<:AbstractPolynomialLike}}, v::AbstractVariable)
 
 Returns the minimal degree of the monomials of `p` in the variable `v`, i.e. `minimum(degree.(terms(p), v))`.
 
 ### Examples
 Calling `mindegree` on on ``4x^2y + xy + 2x`` should return 1, `mindegree(4x^2y + xy + 2x, x)` should return 1 and  `mindegree(4x^2y + xy + 2x, y)` should return 0.
 """
-function mindegree(X::AbstractVector{<:AbstractTermLike}, args...)
-    return isempty(X) ? 0 : minimum(t -> degree(t, args...), X)
+function mindegree(
+    X::AbstractVector{<:AbstractPolynomialLike},
+    args::Vararg{Any,N},
+) where {N}
+    return isempty(X) ? 0 : minimum(t -> mindegree(t, args...), X)
 end
-function mindegree(p::AbstractPolynomialLike, args...)
+function mindegree(p::AbstractPolynomialLike, args::Vararg{Any,N}) where {N}
     return mindegree(terms(p), args...)
 end
+function mindegree(t::AbstractTermLike, args::Vararg{Any,N}) where {N}
+    return degree(t, args...)
+end
 
 #$(SIGNATURES)
 """
@@ -304,30 +310,33 @@ Returns the maximal degree of the monomials of `p` in the variable `v`, i.e. `ma
 Calling `maxdegree` on ``4x^2y + xy + 2x`` should return 3, `maxdegree(4x^2y + xy + 2x, x)` should return 2 and  `maxdegree(4x^2y + xy + 2x, y)` should return 1.
 """
 function maxdegree(
-    X::AbstractVector{<:AbstractTermLike},
+    X::AbstractVector{<:AbstractPolynomialLike},
     args::Vararg{Any,N},
 ) where {N}
-    return mapreduce(t -> degree(t, args...), max, X, init = 0)
+    return mapreduce(t -> maxdegree(t, args...), max, X, init = 0)
 end
 function maxdegree(p::AbstractPolynomialLike, args::Vararg{Any,N}) where {N}
     return maxdegree(terms(p), args...)
 end
+function maxdegree(t::AbstractTermLike, args::Vararg{Any,N}) where {N}
+    return degree(t, args...)
+end
 
 #$(SIGNATURES)
 """
-    extdegree(p::Union{AbstractPolynomialLike, AbstractVector{<:AbstractTermLike}})
+    extdegree(p::Union{AbstractPolynomialLike, AbstractVector{<:AbstractPolynomialLike}})
 
 Returns the extremal total degrees of the monomials of `p`, i.e. `(mindegree(p), maxdegree(p))`.
 
-    extdegree(p::Union{AbstractPolynomialLike, AbstractVector{<:AbstractTermLike}}, v::AbstractVariable)
+    extdegree(p::Union{AbstractPolynomialLike, AbstractVector{<:AbstractPolynomialLike}}, v::AbstractVariable)
 
 Returns the extremal degrees of the monomials of `p` in the variable `v`, i.e. `(mindegree(p, v), maxdegree(p, v))`.
 
 ### Examples
 Calling `extdegree` on ``4x^2y + xy + 2x`` should return `(1, 3)`, `extdegree(4x^2y + xy + 2x, x)` should return `(1, 2)` and  `maxdegree(4x^2y + xy + 2x, y)` should return `(0, 1)`.
 """
 function extdegree(
-    p::Union{AbstractPolynomialLike,AbstractVector{<:AbstractTermLike}},
+    p::Union{AbstractPolynomialLike,AbstractVector{<:AbstractPolynomialLike}},
     args...,
 )
     return (mindegree(p, args...), maxdegree(p, args...))
@@ -465,6 +474,52 @@ function remove_monomials(
     return q
 end
 
+"""
+    function filter_terms(f::Function, p::AbstractPolynomialLike)
+
+Filter the polynomial `p` by only keep the terms `t` such that `f(p)` is
+`true`.
+
+See also [`OfDegree`](@ref).
+
+### Examples
+
+```julia
+julia> p = 1 - 2x + x * y - 3y^2 + x^2 * y
+1 - 2x - 3y² + xy + x²y
+
+julia> filter_terms(OfDegree(2), p)
+-3y² + xy
+
+julia> filter_terms(!OfDegree(2), p)
+1 - 2x + x²y
+
+julia> filter_terms(!OfDegree(0:2), p)
+x²y
+
+julia> filter_terms(iseven ∘ coefficient, p)
+-2x
+```
+"""
+function filter_terms(f::F, p::AbstractPolynomialLike) where {F<:Function}
+    return polynomial(filter(f, terms(p)), SortedUniqState())
+end
+
+"""
+    struct OfDegree{D} <: Function
+        degree::D
+    end
+
+A function `d::OfDegree` is such that `d(t)` returns
+`degree(t) == d.degree`. Note that `!d` creates the negation.
+See also [`filter_terms`](@ref).
+"""
+struct OfDegree{D} <: Function
+    degree::D
+end
+
+(d::OfDegree)(mono::AbstractTermLike) = in(degree(mono), d.degree)
+
 """
     monic(p::AbstractPolynomialLike)
 

src/substitution.jl

@@ -41,23 +41,7 @@ is equivalent to:
     subs(polynomial, (x=>1, y=>2))
 """
 function substitute(st::_AST, p::_APL, s::AbstractMultiSubstitution)
-    return substitute(st, p, pair_zip(_monomial_vector_to_variable_tuple(s)))
-end
-
-# Evaluate the stream
-# If I do s2..., then
-# subs(x, x=>x+y, y=>2) would call substitute(Subs(), x+y, y=>2)
-# subs(x, x=>1, y=>2) would call substitute(Subs(), 1, y=>2)
-# so it would force use to also define
-# sustitute(::_AST, ..., ::AbstractSubstitution...) for Any, _APL and RationalPoly.
-#substitute(st::_AST, p::_APL, s1::AbstractSubstitution, s2::AbstractSubstitution...) = substitute(st, substitute(st, p, s1), s2...)
-function substitute(
-    st::_AST,
-    p::_APL,
-    s1::AbstractSubstitution,
-    s2::AbstractSubstitution...,
-)
-    return substitute(st, substitute(st, p, s1), s2)
+    return substitute(st, p, _flatten_subs(s))
 end
 
 ## Variables
@@ -82,35 +66,70 @@ function powersubstitute(
     return powersubstitute(st, s, p) * powersubstitute(st, s, p2...)
 end
 
-function _promote_subs(S, T, s::Substitution)
-    # `T` is a constant
-    return T
+function _promote_subs(
+    S,
+    ::Type{V},
+    s::Substitution,
+) where {V<:AbstractVariable}
+    return MA.promote_operation(substitute, S, V, typeof(s))
 end
 
-function _promote_subs(S, T::Type{<:Union{RationalPoly,_APL}}, s::Substitution)
-    return MA.promote_operation(substitute, S, T, typeof(s))
+function _flatten_subs(s::AbstractMultiSubstitution)
+    return pair_zip(_monomial_vector_to_variable_tuple(s))
 end
 
-function _promote_subs(S, T, s::AbstractMultiSubstitution)
-    return _promote_subs(
-        S,
-        T,
-        pair_zip(_monomial_vector_to_variable_tuple(s))...,
-    )
+function _flatten_subs(s::Substitution)
+    return (s,)
 end
 
-function _promote_subs(
-    S,
-    T,
-    head::AbstractSubstitution,
+# Turn a tuple of `AbstractSubstitution` into a `Tuple` if `Substitution`
+function _flatten_subs(
+    s::AbstractSubstitution,
     tail::Vararg{AbstractSubstitution,N},
 ) where {N}
-    return _promote_subs(S, _promote_subs(S, T, head), tail...)
+    return (_flatten_subs(s)..., _flatten_subs(tail...)...)
+end
+
+function power_promote(
+    S,
+    ::Type{V},
+    s::Substitutions,
+) where {V<:AbstractVariable}
+    T = MA.promote_operation(substitute, S, V, typeof.(_flatten_subs(s...))...)
+    return MA.promote_operation(*, T, T)
+end
+
+function power_promote(
+    S,
+    ::Vector{V},
+    s::Substitutions,
+) where {V<:AbstractVariable}
+    return power_promote(S, V, s)
+end
+
+function power_promote(
+    S,
+    ::Tuple{V},
+    s::Substitutions,
+) where {V<:AbstractVariable}
+    return power_promote(S, V, s)
+end
+
+function power_promote(
+    S,
+    vars::Tuple{V,Vararg{AbstractVariable,N}},
+    s::Substitutions,
+) where {V<:AbstractVariable,N}
+    return MA.promote_operation(
+        *,
+        power_promote(S, V, s),
+        power_promote(S, Base.tail(vars), s),
+    )
 end
 
 function substitute(st::_AST, m::AbstractMonomial, s::Substitutions)
     if isconstant(m)
-        return one(_promote_subs(typeof(st), typeof(m), s...))
+        return one(power_promote(typeof(st), variables(m), s))
     else
         return powersubstitute(st, s, powers(m)...)
     end

test/monomial_vector.jl

@@ -100,4 +100,11 @@
             v[1]^3,
         ]
     end
+
+    @testset "Negative degree" begin
+        Mod.@polyvar x[1:3]
+        @test_throws ArgumentError monomials(x, -1)
+        @test_throws ArgumentError monomials(x, -1:1)
+        @test_throws ArgumentError monomials(x, [1, -1])
+    end
 end

test/mutable_arithmetics.jl

@@ -1,5 +1,5 @@
-import MutableArithmetics
-const MA = MutableArithmetics
+using Test
+import MutableArithmetics as MA
 
 function all_tests(a, b, c, d, e, f, g)
     a_copy = deepcopy(a)