Skip to content

Commit

Permalink
format f4/sorting.jl a bit
Browse files Browse the repository at this point in the history
  • Loading branch information
@sumiya11
sumiya11 committed Jul 14, 2023
1 parent f0218c8 commit 4b95472
Showing 1 changed file with 113 additions and 123 deletions.
236 changes: 113 additions & 123 deletions src/f4/sorting.jl
View file
Original file line number Diff line number Diff line change
@@ -1,14 +1,12 @@
# Sorting monomials, polynomials, and other things.

# TODO: use a different sorting algorithm for julia v1.9+.
#
# Unstable algorithm should be a bit faster for large arrays
# (basically, quicksort vs. mergesort)
# (basically, quicksort/radixsort vs. mergesort)
_default_sorting_alg() = Base.Sort.DEFAULT_UNSTABLE

# Sorts vector v at indices from:to
# Sorts arr at the range of indices from..to
function sort_part!(
v,
arr,
from::Integer,
to::Integer;
lt=isless,
Expand All @@ -17,111 +15,104 @@ function sort_part!(
rev=false
)
ordr = Base.Sort.ord(lt, by, rev, Base.Sort.Forward)
sort!(v, from, to, alg, ordr)
sort!(arr, from, to, alg, ordr)
end

#------------------------------------------------------------------------------

# Sorts polynomials from the basis
# by their leading monomial in the non-decreasing way
# by the given term ordering.
# Also sorts any arrays passed in `abc` in the same order as basis.
# Sorts polynomials from the basis by their leading monomial in the
# non-decreasing way by the given monomial ordering. Also sorts any arrays
# passed in the `abc` optional argument in the same order.
#
# Returns the sorting permutation.
function sort_polys_by_lead_increasing!(
basis::Basis,
ht::MonomialHashtable,
hashtable::MonomialHashtable,
abc...;
ord::Ord=ht.ord
ord::Ord=hashtable.ord
) where {Ord <: AbstractMonomialOrdering}
gens = basis.monoms
exps = ht.monoms

inds = collect(1:(basis.nfilled))

b_monoms = basis.monoms
h_monoms = hashtable.monoms
permutation = collect(1:(basis.nfilled))
cmps =
(x, y) ->
monom_isless(@inbounds(exps[gens[x][1]]), @inbounds(exps[gens[y][1]]), ord)

sort!(inds, lt=cmps, alg=_default_sorting_alg())

(x, y) -> monom_isless(
@inbounds(h_monoms[b_monoms[x][1]]),
@inbounds(h_monoms[b_monoms[y][1]]),
ord
)
sort!(permutation, lt=cmps, alg=_default_sorting_alg())
# use array assignment insted of elemewise assignment
basis.monoms[1:(basis.nfilled)] = basis.monoms[inds]
basis.coeffs[1:(basis.nfilled)] = basis.coeffs[inds]
for a in abc
a[1:(basis.nfilled)] = a[inds]
# (seems to compile to faster code)
basis.monoms[1:(basis.nfilled)] = basis.monoms[permutation]
basis.coeffs[1:(basis.nfilled)] = basis.coeffs[permutation]
@inbounds for a in abc
@invariant length(a) >= length(permutation)
a[1:(basis.nfilled)] = a[permutation]
end

inds
permutation
end

function is_sorted_by_lead_increasing(
basis::Basis,
ht::MonomialHashtable,
ord::Ord=ht.ord
hashtable::MonomialHashtable,
ord::Ord=hashtable.ord
) where {Ord <: AbstractMonomialOrdering}
gens = basis.monoms
exps = ht.monoms

inds = collect(1:(basis.nfilled))

b_monoms = basis.monoms
h_monoms = hashtable.monoms
permutation = collect(1:(basis.nfilled))
cmps =
(x, y) ->
monom_isless(@inbounds(exps[gens[x][1]]), @inbounds(exps[gens[y][1]]), ord)

issorted(inds, lt=cmps)
(x, y) -> monom_isless(
@inbounds(h_monoms[b_monoms[x][1]]),
@inbounds(h_monoms[b_monoms[y][1]]),
ord
)
issorted(permutation, lt=cmps)
end

#------------------------------------------------------------------------------

# Sorts pairs from pairset in the range [from, from+sz]
# by lcm total degrees in increasing order
# Sorts critical pairs from the pairset in the range from..from+sz by the total
# degree of their lcms in increasing order
#
# Used in update once per one f4 iteration to sort pairs in pairset;
# Also used in normal selection strategy also once per iteration
# Used in update once per one f4 iteration to sort pairs in pairset; Also used
# in normal selection strategy also once per iteration
function sort_pairset_by_degree!(ps::Pairset, from::Int, sz::Int)
sort_part!(ps.pairs, from, from + sz, by=p -> p.deg, alg=_default_sorting_alg())
sort_part!(ps.pairs, from, from + sz, by=pair -> pair.deg)
end

#------------------------------------------------------------------------------

# Sorts the first `npairs` pairs from `pairset` by non-decreasing order of
# the exponent vector of the lcm wrt. the given monomial ordering
function sort_pairset_by_lcm!(pairset::Pairset, npairs::Int, ht::MonomialHashtable)
exps = ht.monoms

cmps = (x, y) -> monom_isless(@inbounds(exps[x.lcm]), @inbounds(exps[y.lcm]), ht.ord)

sort_part!(pairset.pairs, 1, npairs, lt=cmps, alg=_default_sorting_alg())
# Sorts the first `npairs` pairs from `pairset` in the non-decreasing order of
# their lcms by the given monomial ordering
function sort_pairset_by_lcm!(pairset::Pairset, npairs::Int, hashtable::MonomialHashtable)
monoms = hashtable.monoms
cmps =
(pair1, pair2) -> monom_isless(
@inbounds(monoms[pair1.lcm]),
@inbounds(monoms[pair2.lcm]),
hashtable.ord
)
sort_part!(pairset.pairs, 1, npairs, lt=cmps)
end

#------------------------------------------------------------------------------

# Sorts generators selected in normal strategy function
# by their ordering in the current basis (identity sort)
# Sorts generators selected in normal selection strategy by their position in
# the current basis (identity sort)
function sort_generators_by_position!(gens::Vector{Int}, load::Int)
sort_part!(gens, 1, load, alg=_default_sorting_alg())
sort_part!(gens, 1, load)
end

#------------------------------------------------------------------------------
###
# Sorting matrix rows and columns.
# See f4/matrix.jl for details.

# Comparator for matrix rows a and b.
# A row matrix is encoded by the sparse array of
# its nonzero column indices.
# Compare sparse matrix rows a and b.
# A row is an array of integers, which are the indices of nonzero elements
function matrix_row_decreasing_cmp(a::Vector{T}, b::Vector{T}) where {T <: ColumnIdx}
#= a, b - rows as arrays of exponent indices =#
#= a, b - rows as arrays of nonzero indices =#
# va and vb are the leading columns
@inbounds va = a[1]
@inbounds vb = b[1]

if va > vb
return false
end
if va < vb
return true
end

# if same column index => compare the density of rows
# If the same leading column => compare the density of rows
va = length(a)
vb = length(b)
if va > vb
Expand All @@ -130,29 +121,27 @@ function matrix_row_decreasing_cmp(a::Vector{T}, b::Vector{T}) where {T <: Colum
if va < vb
return false
end

# hmm, equal rows?
# should never happen
return false
# Hmm, equal rows?
# This should never happen!
# TODO!!!: add some kind of assertion here
# @invariant false
false
end

# Comparator for matrix rows a and b.
# A row matrix is encoded by the sparse array of
# its nonzero column indices.
# Compare sparse matrix rows a and b.
# A row is an array of integers, which are the indices of nonzero elements
function matrix_row_increasing_cmp(a::Vector{T}, b::Vector{T}) where {T <: ColumnIdx}
#= a, b - rows as arrays of exponent hashes =#
#= a, b - rows as arrays of nonzero indices =#
# va and vb are the leading columns
@inbounds va = a[1]
@inbounds vb = b[1]

if va > vb
return true
end
if va < vb
return false
end

# if same column index => compare the density of rows
# If the same leading column => compare the density of rows
va = length(a)
vb = length(b)
if va > vb
Expand All @@ -161,58 +150,63 @@ function matrix_row_increasing_cmp(a::Vector{T}, b::Vector{T}) where {T <: Colum
if va < vb
return true
end

# hmm, equal rows?
# should never happen
# Hmm, equal rows?
# This should never happen!
# @invariant false
return false
end

# Sort matrix upper rows (polynomial reducers)
# by the leading column index and density.
# Sort matrix upper rows (polynomial reducers) by the leading column index and
# density.
#
# After the sort, the first row will have
# the left-most leading column index and, then, the smallest density.
function sort_matrix_upper_rows_decreasing!(matrix)
# After the sort, the first (smallest) row will have the left-most leading
# column index and, then, the smallest density.
function sort_matrix_upper_rows_decreasing!(matrix::MacaulayMatrix)
#= smaller means pivot being more left =#
#= and density being smaller =#

inds = collect(1:(matrix.nup))
cmp = (x, y) -> matrix_row_decreasing_cmp(@inbounds(matrix.uprows[x]), @inbounds(matrix.uprows[y]))
sort!(inds, lt=cmp, alg=_default_sorting_alg())

matrix.uprows[1:(matrix.nup)] = matrix.uprows[inds]
matrix.up2coef[1:(matrix.nup)] = matrix.up2coef[inds]
# TODO
permutation = collect(1:(matrix.nup))
cmp =
(x, y) -> matrix_row_decreasing_cmp(
@inbounds(matrix.uprows[x]),
@inbounds(matrix.uprows[y])
)
sort!(permutation, lt=cmp, alg=_default_sorting_alg())
matrix.uprows[1:(matrix.nup)] = matrix.uprows[permutation]
matrix.up2coef[1:(matrix.nup)] = matrix.up2coef[permutation]
# TODO: this is a bit hacky
if !isempty(matrix.up2mult)
matrix.up2mult[1:(matrix.nup)] = matrix.up2mult[inds]
matrix.up2mult[1:(matrix.nup)] = matrix.up2mult[permutation]
end
matrix
end

# Sort matrix lower rows (polynomials to be reduced)
# by the leading column index and density.
# Sort matrix lower rows (polynomials to be reduced) by the leading column index
# and density.
#
# After the sort, the first row will have
# the right-most leading column index and, then, the largest density.
function sort_matrix_lower_rows_increasing!(matrix)
# After the sort, the first (smallest) row will have the right-most leading
# column index and, then, the largest density.
function sort_matrix_lower_rows_increasing!(matrix::MacaulayMatrix)
#= smaller means pivot being more right =#
#= and density being larger =#

inds = collect(1:(matrix.nlow))
cmp = (x, y) -> matrix_row_increasing_cmp(@inbounds(matrix.lowrows[x]), @inbounds(matrix.lowrows[y]))
sort!(inds, lt=cmp, alg=_default_sorting_alg())

matrix.lowrows[1:(matrix.nlow)] = matrix.lowrows[inds]
matrix.low2coef[1:(matrix.nlow)] = matrix.low2coef[inds]
permutation = collect(1:(matrix.nlow))
cmp =
(x, y) -> matrix_row_increasing_cmp(
@inbounds(matrix.lowrows[x]),
@inbounds(matrix.lowrows[y])
)
sort!(permutation, lt=cmp, alg=_default_sorting_alg())
matrix.lowrows[1:(matrix.nlow)] = matrix.lowrows[permutation]
matrix.low2coef[1:(matrix.nlow)] = matrix.low2coef[permutation]
# TODO: this is a bit hacky
if !isempty(matrix.low2mult)
matrix.low2mult[1:(matrix.nlow)] = matrix.low2mult[inds]
matrix.low2mult[1:(matrix.nlow)] = matrix.low2mult[permutation]
end
matrix
end

# Sorts the columns of the matrix (encoded in `col2hash` vector)
# by the role of the corresponding monomial in the matrix,
# and then by the current monomial ordering decreasingly.
# Sorts the columns of the matrix (encoded in `col2hash` vector) by the role of
# the corresponding monomial in the matrix, and then by the current monomial
# ordering decreasingly.
#
# See f4/matrix.jl for details
function sort_columns_by_hash!(col2hash::Vector{T}, symbol_ht::MonomialHashtable) where {T}
Expand All @@ -235,8 +229,6 @@ function sort_columns_by_hash!(col2hash::Vector{T}, symbol_ht::MonomialHashtable
sort!(col2hash, lt=ordcmp, alg=_default_sorting_alg())
end

#------------------------------------------------------------------------------

# Given a vector of vectors of exponent vectors and coefficients,
# sort each vector wrt. the given monomial ordering `ord`.
#
Expand All @@ -262,15 +254,13 @@ function sort_input_to_change_ordering!(
nothing
end

#------------------------------------------------------------------------------

function sort_monom_indices_decreasing!(
monoms::Vector{MonomIdx},
cnt::Integer,
ht::MonomialHashtable,
hashtable::MonomialHashtable,
ord::AbstractMonomialOrdering
)
exps = ht.monoms
exps = hashtable.monoms

cmps = (x, y) -> monom_isless(@inbounds(exps[y]), @inbounds(exps[x]), ord)

Expand All @@ -280,10 +270,10 @@ end
function sort_term_indices_decreasing!(
monoms::Vector{MonomIdx},
coeffs::Vector{C},
ht::MonomialHashtable,
hashtable::MonomialHashtable,
ord::AbstractMonomialOrdering
) where {C <: Coeff}
exps = ht.monoms
exps = hashtable.monoms

cmps =
(x, y) -> monom_isless(@inbounds(exps[monoms[y]]), @inbounds(exps[monoms[x]]), ord)
Expand Down

0 comments on commit 4b95472

Please sign in to comment.