Skip to content

Commit

Permalink
Fix monomials for LexOrder (#159)
Browse files Browse the repository at this point in the history 
* Fix monomials for LexOrder

* Fixes

* Fix

* Fix
  • Loading branch information
@blegat
blegat authored May 2, 2024
1 parent af9ccbd commit 91f843f
Show file tree
Hide file tree
Showing 2 changed files with 141 additions and 24 deletions.
103 changes: 79 additions & 24 deletions src/monomial_vector.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -202,7 +257,7 @@ function MonomialVector(
vars = unique!(sort(vars, rev = true))
return MonomialVector(
vars,
getZfordegs(
_all_exponents(
length(vars),
degs,
Commutative,
Expand All @@ -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,
Expand All @@ -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
Expand Down
62 changes: 62 additions & 0 deletions test/comp.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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

0 comments on commit 91f843f

Please sign in to comment.