Fix monomials for LexOrder (#159)

* Fix monomials for LexOrder

* Fixes

* Fix

* Fix

Author: blegat

src/monomial_vector.jl

@@ -12,9 +12,6 @@ struct MonomialVector{V,M} <: AbstractVector{Monomial{V,M}}
         @assert !iscomm(V) || issorted(vars, rev = true)
         @assert all(z -> length(z) == length(vars), Z)
 
-        _isless = let M = M
-            (a, b) -> MP.compare(a, b, M) < 0
-        end
         return new{V,M}(vars, Z)
     end
 end
@@ -125,68 +122,126 @@ function _error_for_negative_degree(deg)
     end
 end
 
-function fillZfordeg!(Z, n, deg, ::Type{Commutative}, filter::Function, ::Int)
+const _Lex = Union{MP.LexOrder,MP.InverseLexOrder}
+
+_last_lex_index(n, ::Type{MP.LexOrder}) = n
+_prev_lex_index(i, ::Type{MP.LexOrder}) = i - 1
+_not_first_indices(n, ::Type{MP.LexOrder}) = n:-1:2
+_last_lex_index(_, ::Type{MP.InverseLexOrder}) = 1
+_prev_lex_index(i, ::Type{MP.InverseLexOrder}) = i + 1
+_not_first_indices(n, ::Type{MP.InverseLexOrder}) = 1:(n-1)
+
+function _fill_exponents!(Z, n, degs, ::Type{Commutative}, M::Type{<:_Lex}, filter::Function)
+    _error_for_negative_degree.(degs)
+    maxdeg = maximum(degs, init = 0)
+    I = _not_first_indices(n, M)
+    z = zeros(Int, n)
+    while true
+        deg = sum(z)
+        if deg in degs && filter(z)
+            push!(Z, z)
+            z = copy(z)
+        end
+        if deg == maxdeg
+            i = findfirst(i -> !iszero(z[i]), I)
+            if isnothing(i)
+                break
+            end
+            j = I[i]
+            z[j] = 0
+            z[_prev_lex_index(j, M)] += 1
+        else
+            z[_last_lex_index(n, M)] += 1
+        end
+    end
+end
+
+function _fill_exponents!(Z, n, deg, ::Type{Commutative}, M::Type{<:_Lex}, filter::Function, ::Int)
     _error_for_negative_degree(deg)
+    I = _not_first_indices(n, M)
     z = zeros(Int, n)
-    z[end] = deg
+    z[_last_lex_index(n, M)] = deg
     while true
         if filter(z)
             push!(Z, z)
             z = copy(z)
         end
-        if z[1] == deg
+        i = findfirst(i -> !iszero(z[i]), I)
+        if isnothing(i)
             break
         end
-        i = findfirst(i -> !iszero(z[i]), n:-1:2)
-        j = (n:-1:2)[i]
+        j = I[i]
         p = z[j]
         z[j] = 0
-        z[end] = p - 1
-        z[j-1] += 1
+        z[_last_lex_index(n, M)] = p - 1
+        z[_prev_lex_index(j, M)] += 1
     end
 end
-function fillZrec!(Z, z, i, n, deg, filter::Function)
+
+function _fill_noncomm_exponents_rec!(Z, z, i, n, deg, ::Type{MP.LexOrder}, filter::Function)
     if deg == 0
         if filter(z)
             push!(Z, copy(z))
         end
     else
         for i in i:i+n-1
             z[i] += 1
-            fillZrec!(Z, z, i, n, deg - 1, filter)
+            _fill_noncomm_exponents_rec!(Z, z, i, n, deg - 1, LexOrder, filter)
             z[i] -= 1
         end
     end
 end
-function fillZfordeg!(
+
+function _fill_exponents!(
     Z,
     n,
     deg,
     ::Type{NonCommutative},
+    ::Type{MP.LexOrder},
     filter::Function,
     maxdeg::Int,
 )
     _error_for_negative_degree(deg)
     _error_for_negative_degree(maxdeg)
     z = zeros(Int, maxdeg * n - maxdeg + 1)
     start = length(Z) + 1
-    fillZrec!(Z, z, 1, n, deg, filter)
+    _fill_noncomm_exponents_rec!(Z, z, 1, n, deg, MP.LexOrder, filter)
     return reverse!(view(Z, start:length(Z)))
 end
-# List exponents in decreasing Graded Lexicographic Order
-function getZfordegs(
+
+function _fill_exponents!(Z, n, deg, ::Type{V}, ::Type{MP.Reverse{M}}, args...) where {V,M}
+    prev = lastindex(Z)
+    _fill_exponents!(Z, n, deg, V, M, args...)
+    reverse!(view(Z, (prev + 1):lastindex(Z)))
+    return
+end
+
+function _fill_exponents!(
+    Z::Vector{Vector{Int}},
     n,
     degs::AbstractVector{Int},
     ::Type{V},
-    ::Type{M},
+    ::Type{MP.Graded{M}},
     filter::Function,
 ) where {V,M}
-    Z = Vector{Vector{Int}}()
     # For non-commutative, lower degree need to create a vector of exponent as large as for the highest degree
-    maxdeg = isempty(degs) ? 0 : maximum(degs)
+    maxdeg = maximum(degs, init = 0)
     for deg in sort(degs)
-        fillZfordeg!(Z, n, deg, V, filter, maxdeg)
+        _fill_exponents!(Z, n, deg, V, M, filter, maxdeg)
     end
+    return
+end
+
+# List exponents in decreasing Graded Lexicographic Order
+function _all_exponents(
+    n,
+    degs::AbstractVector{Int},
+    ::Type{V},
+    ::Type{M},
+    filter::Function,
+) where {V,M}
+    Z = Vector{Vector{Int}}()
+    _fill_exponents!(Z, n, degs, V, M, filter)
     _isless = let M = M
         (a, b) -> MP.compare(a, b, M) < 0
     end
@@ -202,7 +257,7 @@ function MonomialVector(
     vars = unique!(sort(vars, rev = true))
     return MonomialVector(
         vars,
-        getZfordegs(
+        _all_exponents(
             length(vars),
             degs,
             Commutative,
@@ -222,7 +277,7 @@ function MonomialVector(
     filter::Function = x -> true,
 ) where {M}
     vars = unique!(sort(vars, rev = true))
-    Z = getZfordegs(
+    Z = _all_exponents(
         length(vars),
         degs,
         NonCommutative,
@@ -248,11 +303,11 @@ function MP.monomials(vars::Tuple{Vararg{Variable}}, args...)
 end
 
 #function MP.monomials(vars::TupOrVec{Variable{true}}, degs::AbstractVector{Int}, filter::Function = x->true)
-#    Z = getZfordegs(length(vars), degs, true, z -> filter(Monomial(vars, z)))
+#    Z = _all_exponents(length(vars), degs, true, z -> filter(Monomial(vars, z)))
 #    [Monomial{true}(vars, z) for z in Z]
 #end
 #function MP.monomials(vars::TupOrVec{<:Variable{<:NonCommutative}}, degs::AbstractVector{Int}, filter::Function = x->true)
-#    Z = getZfordegs(length(vars), degs, false, z -> filter(Monomial(vars, z)))
+#    Z = _all_exponents(length(vars), degs, false, z -> filter(Monomial(vars, z)))
 #    v = isempty(Z) ? vars : getvarsforlength(vars, length(first(Z)))
 #    [Monomial(v, z) for z in Z]
 #end

test/comp.jl

@@ -44,3 +44,65 @@ end
     @test ordering(x[1]) == order
     @test issorted(monomials(x[1], 0:2))
 end
+
+function _test_less(a, b)
+    @test a < b
+    @test b > a
+end
+
+function _test_monomials(vars, degs, exp)
+    # Without `collect`, `exp` is promoted to a `MonomialVector`
+    # which sorts it so it doesn't test the order
+    @test collect(monomials(vars, degs)) == exp
+end
+
+@testset "LexOrder" begin
+    @polyvar x y monomial_order = LexOrder
+    _test_less(y, y^2)
+    _test_less(x^0, y)
+    _test_less(y^2, x)
+    _test_less(x * y^2, x^2)
+    _test_monomials([x, y], 3, [y^3, y^2 * x, y * x^2, x^3])
+    _test_monomials([x, y], 1:2, [y, y^2, x, x * y, x^2])
+    _test_monomials([x, y], [0, 1, 3], [1, y, y^3, x, y^2 * x, y * x^2, x^3])
+end
+
+@testset "InverseLexOrder" begin
+    @polyvar x y monomial_order = InverseLexOrder
+    _test_less(y, y^2)
+    _test_less(x^0, y)
+    _test_less(x, y^2)
+    _test_less(x^2, x * y^2)
+    _test_monomials([x, y], 3, [x^3, x^2 * y, x * y^2, y^3])
+end
+
+@testset "Reverse{LexOrder}" begin
+    @polyvar x y monomial_order = Reverse{LexOrder}
+    _test_less(y^2, y)
+    _test_less(y, x^0)
+    _test_less(x, y^2)
+    _test_less(x^2, x * y^2)
+    _test_monomials([x, y], 3, [x^3, x^2 * y, x * y^2, y^3])
+end
+
+@testset "Reverse{InverseLexOrder}" begin
+    @polyvar x y monomial_order = Reverse{InverseLexOrder}
+    _test_less(y^2, y)
+    _test_less(y, x^0)
+    _test_less(y^2, x)
+    _test_less(x * y^2, x^2)
+    _test_monomials([x, y], 3, [y^3, y^2 * x, y * x^2, x^3])
+    _test_monomials([x, y], 1:2, [y^2, x * y, y, x^2, x])
+    _test_monomials([x, y], [0, 1, 3], [y^3, x*y^2, x^2*y, y, x^3, x, 1])
+end
+
+@testset "Graded{Reverse{InverseLexOrder}}" begin
+    @polyvar x y monomial_order = Graded{Reverse{InverseLexOrder}}
+    _test_less(y, y^2)
+    _test_less(x^0, y)
+    _test_less(x, y^2)
+    _test_less(x^2, x * y^2)
+    _test_monomials([x, y], 3, [y^3, y^2 * x, y * x^2, x^3])
+    _test_monomials([x, y], 1:2, [y, x, y^2, x * y, x^2])
+    _test_monomials([x, y], [0, 1, 3], [1, y, x, y^3, x*y^2, x^2*y, x^3])
+end