From 8ca687c72d211e9e97c7910d0ccfc20d58ca418e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Tue, 30 Apr 2024 10:08:08 +0200 Subject: [PATCH 1/4] Fix monomials for LexOrder --- src/monomial_vector.jl | 70 +++++++++++++++++++++++++++++++++--------- test/comp.jl | 16 ++++++++++ 2 files changed, 71 insertions(+), 15 deletions(-) diff --git a/src/monomial_vector.jl b/src/monomial_vector.jl index 7186503..b2a73f1 100644 --- a/src/monomial_vector.jl +++ b/src/monomial_vector.jl @@ -125,7 +125,31 @@ function _error_for_negative_degree(deg) end end -function fillZfordeg!(Z, n, deg, ::Type{Commutative}, filter::Function, ::Int) +function _fill_exponents!(Z, n, degs, ::Type{Commutative}, ::Type{MP.LexOrder}, filter::Function) + _error_for_negative_degree.(degs) + maxdeg = maximum(degs, init = 0) + 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]), n:-1:2) + if isnothing(i) + break + end + j = (n:-1:2)[i] + z[j] = 0 + z[j-1] += 1 + else + z[end] += 1 + end + end +end + +function _fill_exponents!(Z, n, deg, ::Type{Commutative}, ::Type{MP.LexOrder}, filter::Function, ::Int) _error_for_negative_degree(deg) z = zeros(Int, n) z[end] = deg @@ -134,10 +158,10 @@ function fillZfordeg!(Z, n, deg, ::Type{Commutative}, filter::Function, ::Int) push!(Z, z) z = copy(z) end - if z[1] == deg + i = findfirst(i -> !iszero(z[i]), n:-1:2) + if isnothing(i) break end - i = findfirst(i -> !iszero(z[i]), n:-1:2) j = (n:-1:2)[i] p = z[j] z[j] = 0 @@ -145,7 +169,7 @@ function fillZfordeg!(Z, n, deg, ::Type{Commutative}, filter::Function, ::Int) z[j-1] += 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)) @@ -153,39 +177,55 @@ function fillZrec!(Z, z, i, n, deg, filter::Function) else for i in i:i+n-1 z[i] += 1 - fillZrec!(Z, z, i, n, deg - 1, filter) + fillZrec!(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{M}, filter::Function, maxdeg::Int, -) +) where {M} _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) + fillZrec!(Z, z, 1, n, deg, M, filter) return reverse!(view(Z, start:length(Z))) end # List exponents in decreasing Graded Lexicographic Order -function getZfordegs( +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 + @assert issorted(Z, lt = _isless) + return Z +end +function _all_exponents( + n, + degs::AbstractVector{Int}, + ::Type{V}, + ::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 _isless = let M = M (a, b) -> MP.compare(a, b, M) < 0 @@ -202,7 +242,7 @@ function MonomialVector( vars = unique!(sort(vars, rev = true)) return MonomialVector( vars, - getZfordegs( + _all_exponents( length(vars), degs, Commutative, @@ -222,7 +262,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 +288,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 diff --git a/test/comp.jl b/test/comp.jl index 2e72721..dfbcc9b 100644 --- a/test/comp.jl +++ b/test/comp.jl @@ -44,3 +44,19 @@ 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 + +@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 From 06f7e1b2689b9645efc32442acf8cd1e48d311bd Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Tue, 30 Apr 2024 14:10:02 +0200 Subject: [PATCH 2/4] Fixes --- src/monomial_vector.jl | 77 +++++++++++++++++++++++++----------------- test/comp.jl | 52 ++++++++++++++++++++++++++-- 2 files changed, 95 insertions(+), 34 deletions(-) diff --git a/src/monomial_vector.jl b/src/monomial_vector.jl index b2a73f1..2fe71eb 100644 --- a/src/monomial_vector.jl +++ b/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,9 +122,19 @@ function _error_for_negative_degree(deg) end end -function _fill_exponents!(Z, n, degs, ::Type{Commutative}, ::Type{MP.LexOrder}, filter::Function) +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) @@ -136,39 +143,41 @@ function _fill_exponents!(Z, n, degs, ::Type{Commutative}, ::Type{MP.LexOrder}, z = copy(z) end if deg == maxdeg - i = findfirst(i -> !iszero(z[i]), n:-1:2) + i = findfirst(i -> !iszero(z[i]), I) if isnothing(i) break end - j = (n:-1:2)[i] + j = I[i] z[j] = 0 - z[j-1] += 1 + z[_prev_lex_index(j, M)] += 1 else - z[end] += 1 + z[_last_lex_index(n, M)] += 1 end end end -function _fill_exponents!(Z, n, deg, ::Type{Commutative}, ::Type{MP.LexOrder}, filter::Function, ::Int) +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 - i = findfirst(i -> !iszero(z[i]), n:-1:2) + i = findfirst(i -> !iszero(z[i]), I) if isnothing(i) break end - 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 _fill_noncomm_exponents_rec!(Z, z, i, n, deg, ::Type{MP.LexOrder}, filter::Function) if deg == 0 if filter(z) @@ -182,15 +191,16 @@ function _fill_noncomm_exponents_rec!(Z, z, i, n, deg, ::Type{MP.LexOrder}, filt end end end + function _fill_exponents!( Z, n, deg, ::Type{NonCommutative}, - ::Type{M}, + ::Type{MP.LexOrder}, filter::Function, maxdeg::Int, -) where {M} +) _error_for_negative_degree(deg) _error_for_negative_degree(maxdeg) z = zeros(Int, maxdeg * n - maxdeg + 1) @@ -198,35 +208,40 @@ function _fill_exponents!( fillZrec!(Z, z, 1, n, deg, M, filter) return reverse!(view(Z, start:length(Z))) end -# List exponents in decreasing Graded Lexicographic Order -function _all_exponents( + +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}}() - _fill_exponents!(Z, n, degs, V, M, filter) - _isless = let M = M - (a, b) -> MP.compare(a, b, M) < 0 + # For non-commutative, lower degree need to create a vector of exponent as large as for the highest degree + maxdeg = maximum(degs, init = 0) + for deg in sort(degs) + _fill_exponents!(Z, n, deg, V, M, filter, maxdeg) end - @assert issorted(Z, lt = _isless) - return Z + return end + +# List exponents in decreasing Graded Lexicographic Order function _all_exponents( n, degs::AbstractVector{Int}, ::Type{V}, - ::Type{MP.Graded{M}}, + ::Type{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 = maximum(degs, init = 0) - for deg in sort(degs) - _fill_exponents!(Z, n, deg, V, M, filter, maxdeg) - end + _fill_exponents!(Z, n, degs, V, M, filter) _isless = let M = M (a, b) -> MP.compare(a, b, M) < 0 end diff --git a/test/comp.jl b/test/comp.jl index dfbcc9b..ff425c8 100644 --- a/test/comp.jl +++ b/test/comp.jl @@ -50,13 +50,59 @@ function _test_less(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] + _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 From e0f95947bc31bf77d41d5a9eb60fb8787a987eda Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Thu, 2 May 2024 15:23:36 +0200 Subject: [PATCH 3/4] Fix --- src/monomial_vector.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/monomial_vector.jl b/src/monomial_vector.jl index 2fe71eb..3979ade 100644 --- a/src/monomial_vector.jl +++ b/src/monomial_vector.jl @@ -186,7 +186,7 @@ function _fill_noncomm_exponents_rec!(Z, z, i, n, deg, ::Type{MP.LexOrder}, filt else for i in i:i+n-1 z[i] += 1 - fillZrec!(Z, z, i, n, deg - 1, LexOrder, filter) + _fill_noncomm_exponents_rec!(Z, z, i, n, deg - 1, LexOrder, filter) z[i] -= 1 end end @@ -205,7 +205,7 @@ function _fill_exponents!( _error_for_negative_degree(maxdeg) z = zeros(Int, maxdeg * n - maxdeg + 1) start = length(Z) + 1 - fillZrec!(Z, z, 1, n, deg, M, filter) + _fill_noncomm_exponents_rec!(Z, z, 1, n, deg, M, filter) return reverse!(view(Z, start:length(Z))) end From 30d73a52960d70982cd2d55a90488425ed9d87c3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Thu, 2 May 2024 17:15:12 +0200 Subject: [PATCH 4/4] Fix --- src/monomial_vector.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/monomial_vector.jl b/src/monomial_vector.jl index 3979ade..991f4b1 100644 --- a/src/monomial_vector.jl +++ b/src/monomial_vector.jl @@ -205,7 +205,7 @@ function _fill_exponents!( _error_for_negative_degree(maxdeg) z = zeros(Int, maxdeg * n - maxdeg + 1) start = length(Z) + 1 - _fill_noncomm_exponents_rec!(Z, z, 1, n, deg, M, filter) + _fill_noncomm_exponents_rec!(Z, z, 1, n, deg, MP.LexOrder, filter) return reverse!(view(Z, start:length(Z))) end