diff --git a/docs/src/free_associative_algebra.md b/docs/src/free_associative_algebra.md index 4a241436bc..cca5add4df 100644 --- a/docs/src/free_associative_algebra.md +++ b/docs/src/free_associative_algebra.md @@ -1,5 +1,5 @@ ```@meta -CurrentModule = AbstractAlgebra +CurrentModule = AbstractAlgebra.Generic DocTestSetup = quote using AbstractAlgebra end @@ -218,3 +218,35 @@ julia> collect(exponent_words(3*b*a*c - b + c + 2)) [] ``` +### Groebner bases + +The function `groebner_basis` provides the computation of a Groebner basis of an ideal, given a set of +generators of that ideal. +Since such a Groebner basis is not necessarily finite, one can additionally pass a `reduction_bound` +to the function, to only compute a partial Groebner basis. + +**Examples** + +```jldoctest; setup = :(using AbstractAlgebra) +julia> R, (x, y, u, v, t, s) = free_associative_algebra(GF(2), ["x", "y", "u", "v", "t", "s"]) +(Free associative algebra on 6 indeterminates over finite field F_2, AbstractAlgebra.Generic.FreeAssAlgElem{AbstractAlgebra.GFElem{Int64}}[x, y, u, v, t, s]) + +julia> g = Generic.groebner_basis([u*(x*y)^3 + u*(x*y)^2 + u + v, (y*x)^3*t + (y*x)^2*t + t + s]) +5-element Vector{AbstractAlgebra.Generic.FreeAssAlgElem{AbstractAlgebra.GFElem{Int64}}}: + u*x*y*x*y*x*y + u*x*y*x*y + u + v + y*x*y*x*y*x*t + y*x*y*x*t + t + s + u*x*s + v*x*t + u*x*y*x*s + v*x*y*x*t + u*x*y*x*y*x*s + v*x*y*x*y*x*t +``` + +In order to check whether a given element of the algebra is in the ideal generated by a Groebner +basis `g`, one can compute its normal form. +```jldoctest; setup = :(using AbstractAlgebra) +julia> R, (x, y, u, v, t, s) = free_associative_algebra(GF(2), ["x", "y", "u", "v", "t", "s"]); + +julia> g = Generic.groebner_basis([u*(x*y)^3 + u*(x*y)^2 + u + v, (y*x)^3*t + (y*x)^2*t + t + s]); + +julia> normal_form(u*(x*y)^3*s*t + u*(x*y)^2*s*t +u*s*t + v*s*t, g) +0 + ``` diff --git a/src/AbstractAlgebra.jl b/src/AbstractAlgebra.jl index 30c8ea81f6..4f63b1ab45 100644 --- a/src/AbstractAlgebra.jl +++ b/src/AbstractAlgebra.jl @@ -644,12 +644,14 @@ import .Generic: finish import .Generic: fit! import .Generic: gcd import .Generic: gcdx +import .Generic: groebner_basis import .Generic: has_bottom_neighbor import .Generic: has_left_neighbor import .Generic: hash import .Generic: hooklength import .Generic: image_fn import .Generic: image_map +import .Generic: interreduce! import .Generic: intersection import .Generic: inv! import .Generic: inverse_fn diff --git a/src/Generic.jl b/src/Generic.jl index 6cdbe71715..da0f998707 100644 --- a/src/Generic.jl +++ b/src/Generic.jl @@ -327,6 +327,10 @@ include("generic/MapCache.jl") include("generic/Ideal.jl") +include("generic/AhoCorasick.jl") + +include("generic/FreeAssAlgebraGroebner.jl") + ############################################################################### # # Temporary miscellaneous files being moved from Hecke.jl diff --git a/src/generic/AhoCorasick.jl b/src/generic/AhoCorasick.jl new file mode 100644 index 0000000000..93044f3395 --- /dev/null +++ b/src/generic/AhoCorasick.jl @@ -0,0 +1,255 @@ +# +# FreeAssAhoCorasick.jl : implement bulk divide check for leading terms of free associative Algebra elements +# for use e.g. in Groebner Basis computation +# +############################################################################### + +export aho_corasick_automaton + +export AhoCorasickAutomaton + +export AhoCorasickMatch + +export insert_keyword! + +export search + +export last_position + +export keyword_index + +export keyword + +#export Word + +const Word = Vector{Int} + +struct Queue{T} + data::Vector{T} +end + +function Queue{T}() where T + return Queue{T}(T[]) +end + +function enqueue!(q::Queue{T}, val::T) where T + push!(q.data, val) +end +function dequeue!(q::Queue) + return popfirst!(q.data) +end +isempty(q::Queue) = isempty(q.data) + +@doc """ + AhoCorasickAutomaton + +An Aho-Corasick automaton, which can be used to efficiently search for a fixed list of keywords (vectors of Ints) in +arbitrary lists of integers + +# Examples +```jldoctest; setup = :(using AbstractAlgebra) +julia> keywords = [[1, 2, 3, 4], [1, 5, 4], [4, 1, 2], [1, 2]]; + +julia> aut = Generic.aho_corasick_automaton(keywords); + +julia> Generic.search(aut, [10, 4, 1, 2, 3, 4]) +AbstractAlgebra.Generic.AhoCorasickMatch(6, 1, [1, 2, 3, 4]) +``` +""" +mutable struct AhoCorasickAutomaton + goto::Vector{Dict{Int,Int}} + fail::Vector{Int} + """ + Output stores for each node a tuple (i, k), where i is the index of the keyword k in + the original list of keywords. If several keywords would be the output of the node, only + the one with the smallest index is stored + """ + output::Vector{Tuple{Int,Word}} +end + +@doc """ + AhoCorasickMatch(last_position::Int, keyword_index::Int, keyword::Vector{Int}) + +The return value of searching in a given word with an AhoCorasickAutomaton. Contains the position of the last letter in +the word that matches a keyword in the automaton, an index of the keyword that was matched and the keyword itself. + +# Examples +```jldoctest; setup = :(using AbstractAlgebra) +julia> keywords = [[1, 2, 3, 4], [1, 5, 4], [4, 1, 2], [1, 2]]; + +julia> aut = Generic.aho_corasick_automaton(keywords); + +julia> result = Generic.search(aut, [10, 4, 1, 2, 3, 4]) +AbstractAlgebra.Generic.AhoCorasickMatch(6, 1, [1, 2, 3, 4]) + +julia> Generic.last_position(result) +6 + +julia> Generic.keyword_index(result) +1 + +julia> Generic.keyword(result) +4-element Vector{Int64}: + 1 + 2 + 3 + 4 +``` +""" +struct AhoCorasickMatch + last_position::Int + keyword_index::Int + keyword::Word +end + +""" +returns the last position of the match in the word that was searched +""" +function last_position(match::AhoCorasickMatch) + return match.last_position +end + +""" +returns the index of the keyword in the corresponding aho corasick automaton +""" +function keyword_index(match::AhoCorasickMatch) + return match.keyword_index +end + +""" +returns the keyword corresponding to the match +""" +function keyword(match::AhoCorasickMatch) + return match.keyword +end + +function aho_corasick_match(last_position::Int, keyword_index::Int, keyword::Word) + return AhoCorasickMatch(last_position, keyword_index, keyword) +end + +Base.hash(m::AhoCorasickMatch, h::UInt) = hash(m.last_position, hash(m.keyword_index, + hash(m.keyword, h))) +function ==(m1::AhoCorasickMatch, m2::AhoCorasickMatch) + return m1.last_position == m2.last_position && + m1.keyword_index == m2.keyword_index && + m1.keyword == m2.keyword +end + +function AhoCorasickAutomaton(keywords::Vector{Word}) + automaton = AhoCorasickAutomaton([], [], []) + construct_goto!(automaton, keywords) + construct_fail!(automaton) + return automaton +end + +function aho_corasick_automaton(keywords::Vector{Word}) + return AhoCorasickAutomaton(keywords) +end + +function lookup(automaton::AhoCorasickAutomaton, current_state::Int, next_letter::Int) + ret_value = get(automaton.goto[current_state], next_letter, nothing) + if current_state == 1 && isnothing(ret_value) + return 1 + end + return ret_value +end + + +function Base.length(automaton::AhoCorasickAutomaton) + return length(automaton.goto) +end + +function new_state!(automaton) + push!(automaton.goto, Dict{Int,Int}()) + push!(automaton.output, (typemax(Int), [])) + push!(automaton.fail, 1) + return length(automaton.goto) +end + +function enter!(automaton::AhoCorasickAutomaton, keyword::Word, current_index) + current_state = 1 + for c in keyword + current_state = get!(automaton.goto[current_state], c) do + new_state!(automaton) + end + end + if automaton.output[current_state][1] > current_index + automaton.output[current_state] = (current_index, keyword) + end +end + +function construct_goto!(automaton::AhoCorasickAutomaton, keywords::Vector{Word}) + new_state!(automaton) + for (current_index, keyword) in enumerate(keywords) + enter!(automaton, keyword, current_index) + end +end + +function construct_fail!(automaton::AhoCorasickAutomaton) + q = Queue{Int}() + for v in values(automaton.goto[1]) + enqueue!(q, v) + end + while !isempty(q) + current_state = dequeue!(q) + for k in keys(automaton.goto[current_state]) + new_state = lookup(automaton, current_state, k) + enqueue!(q, new_state) + state = automaton.fail[current_state] + while isnothing(lookup(automaton, state, k)) + state = automaton.fail[state] + end + automaton.fail[new_state] = lookup(automaton, state, k) + if automaton.output[new_state][1] > + automaton.output[automaton.fail[new_state]][1] + automaton.output[new_state] = automaton.output[automaton.fail[new_state]] + end + + end + end +end + +@doc """ + insert_keyword!(aut::AhoCorasickAutomaton, keyword::Word, index::Int) + +Insert a new keyword into a given Aho-Corasick automaton to avoid having to rebuild the entire +automaton. +""" +function insert_keyword!(aut::AhoCorasickAutomaton, keyword::Word, index::Int) + enter!(aut, keyword, index) + aut.fail = ones(Int, length(aut.goto)) + construct_fail!(aut) +end + +@doc """ + search(automaton::AhoCorasickAutomaton, word::Word) + +Search for the first occurrence of a keyword that is stored in `automaton` in the given `word`. +""" +function search(automaton::AhoCorasickAutomaton, word::Word) + current_state = 1 + result = AhoCorasickMatch(typemax(Int), typemax(Int), []) + for i in 1:length(word) + c = word[i] + while true + next_state = lookup(automaton, current_state, c) + if !isnothing(next_state) + current_state = next_state + break + else + current_state = automaton.fail[current_state] + end + end + if automaton.output[current_state][1] < result.keyword_index + result = AhoCorasickMatch( + i, + automaton.output[current_state][1], + automaton.output[current_state][2], + ) + end + end + if isempty(result.keyword) + return nothing + end + return result +end diff --git a/src/generic/FreeAssAlgebra.jl b/src/generic/FreeAssAlgebra.jl index b6addbd370..2f18087a2b 100644 --- a/src/generic/FreeAssAlgebra.jl +++ b/src/generic/FreeAssAlgebra.jl @@ -370,6 +370,62 @@ function combine_like_terms!(z::FreeAssAlgElem{T}) where T return z end + +@doc """ + isless(p::FreeAssAlgElem{T}, q::FreeAssAlgElem{T}) where T + +Implements the degree lexicographic ordering on terms, i.e. +first, the degrees of the largest monomials are compared, and if they +are the same, they are compared lexicographically and if they are still the same, +the coefficients are compared. +If everything is still the same, the next largest monomial is compared +and lastly the number of monomials is compared. +Since the coefficients are also compared, this only works when the +coefficient Ring implements isless. + +The order of letters is the reverse of the order given when initialising the algebra. + +# Examples +```jldoctest; setup = :(using AbstractAlgebra) +julia> R, (x, y) = free_associative_algebra(QQ, ["x", "y"]); + +julia> x < y^2 +true + +julia> x^2 < x^2 + y +true + +julia> y < x +true + +julia> x^2 < 2*x^2 +true +``` +""" +function isless(p::FreeAssAlgElem{T}, q::FreeAssAlgElem{T}) where T + if p == q + return false + end + l = min(length(p.exps), length(q.exps)) + sort_terms!(p) + sort_terms!(q) + for i in 1:l + c = word_cmp(q.exps[i], p.exps[i]) + if c > 0 + return true + elseif c < 0 + return false + elseif p.coeffs[i] != q.coeffs[i] + return p.coeffs[i] < q.coeffs[i] + end + end + if length(p.exps) < length(q.exps) + return true + else + return false + end +end + ############################################################################### # # Arithmetic @@ -511,7 +567,6 @@ function mul_term(c::T, w::Vector{Int}, a::FreeAssAlgElem{T}, wp::Vector{Int}) w return FreeAssAlgElem{T}(parent(a), zcoeffs, zexps, a.length) end - # return (true, l, r) with a = l*b*r and length(l) minimal # or (false, junk, junk) if a is not two-sided divisible by b function word_divides_leftmost(a::Vector{Int}, b::Vector{Int}) diff --git a/src/generic/FreeAssAlgebraGroebner.jl b/src/generic/FreeAssAlgebraGroebner.jl new file mode 100644 index 0000000000..791480e2b4 --- /dev/null +++ b/src/generic/FreeAssAlgebraGroebner.jl @@ -0,0 +1,642 @@ +############################################################################### +# +# FreeAssAlgebraGroebner.jl : free associative algebra R Groebner basis +# +############################################################################### + +export groebner_basis + +export interreduce! + +export normal_form + +export normal_form_weak + +include("PriorityQueue.jl") + +const Monomial = Vector{Int} + +abstract type Obstruction{T} end +""" +Represents the overlap of the leading term of the two polynomials +`first_poly` and `second_poly`. Here, `first_index` and `second_index` +are the indices of `first_poly` and `second_poly` respectively +in the Groebner basis. +The first and second prefix and suffix satisfy +first_prefix g_i first_suffix = second_prefix g_j second_suffix +where `g_i` is `first_poly` and `g_j` is `second_poly`. +""" +struct ObstructionTriple{T} <: Obstruction{T} + first_poly::FreeAssAlgElem{T} + second_poly::FreeAssAlgElem{T} + first_prefix::Monomial + first_suffix::Monomial + second_prefix::Monomial + second_suffix::Monomial + first_index::Int + second_index::Int +end + +function ObstructionTriple{T}(first_poly::FreeAssAlgElem{T}, + second_poly::FreeAssAlgElem{T}, + pre_and_suffixes::NTuple{4, Monomial}, + first_index::Int, + second_index::Int + ) where T + return ObstructionTriple{T}(first_poly, second_poly, pre_and_suffixes[1], + pre_and_suffixes[2], pre_and_suffixes[3], + pre_and_suffixes[4], first_index, second_index) +end + +function FreeAssAlgElem{T}(R::FreeAssAlgebra{T}, mon::Monomial) where T + return FreeAssAlgElem{T}(R, [one(base_ring(R))], [mon], 1) +end + +""" +takes an obstruction triple (p, q, o) and returns the common multiple +of the leading terms of p and q defined by o +TODO documentation +""" +function common_multiple_leading_term(ot::ObstructionTriple{T}) where T + return FreeAssAlgElem{T}(parent(ot.first_poly), ot.first_prefix) * + FreeAssAlgElem{T}(parent(ot.first_poly), _leading_word(ot.first_poly)) * + FreeAssAlgElem{T}(parent(ot.first_poly), ot.first_suffix) +end + +function s_polynomial(ot::ObstructionTriple{T}) where T + first_term = + FreeAssAlgElem{T}(parent(ot.first_poly), ot.first_prefix) * + ot.first_poly * + FreeAssAlgElem{T}(parent(ot.first_poly), ot.first_suffix) + second_term = + FreeAssAlgElem{T}(parent(ot.first_poly), ot.second_prefix) * + ot.second_poly * + FreeAssAlgElem{T}(parent(ot.first_poly), ot.second_suffix) + return inv(leading_coefficient(ot.first_poly)) * first_term - + inv(leading_coefficient(ot.second_poly)) * second_term +end + +# skip all of the extra length-checking +function _leading_word(a::FreeAssAlgElem{T}) where T + return a.exps[1] +end + +@doc """ + gb_divides_leftmost(a::Word, aut::AhoCorasickAutomaton) + +If an element of the Groebner basis that is stored in `aut` divides `a`, +return (true, a1, a2, keyword_index), where `keyword_index` is the index of the +keyword that divides `a` such that `a = a1 aut[keyword_index] a2`. +""" +function gb_divides_leftmost(a::Word, aut::AhoCorasickAutomaton) + match = search(aut, a) + if isnothing(match) + return (false, Word(), Word(), -1) + end + return ( + true, + a[1:(match.last_position - length(match.keyword))], + a[(match.last_position + 1):end], + match.keyword_index, + ) +end + +# implementation of the normal form function using aho corasick to check for all Groebner basis elements in parallel +@doc """ + normal_form(f::FreeAssAlgElem{T}, g::Vector{FreeAssAlgElem{T}}, aut::AhoCorasickAutomaton) + +Assuming `g` is a Groebner basis and `aut` an Aho-Corasick automaton for the elements of `g`, +compute the normal form of `f` with respect to `g` +""" +function normal_form( + f::FreeAssAlgElem{T}, + g::Vector{FreeAssAlgElem{T}}, + aut::AhoCorasickAutomaton, +) where T + R = parent(f) + rexps = Monomial[] + rcoeffs = T[] + while length(f) > 0 + ok, left, right, match_index = gb_divides_leftmost(f.exps[1], aut) + if ok + qi = divexact(f.coeffs[1], g[match_index].coeffs[1]) + f = _sub_rest(f, mul_term(qi, left, g[match_index], right), 1) + else + push!(rcoeffs, f.coeffs[1]) + push!(rexps, f.exps[1]) + f = FreeAssAlgElem{T}(R, f.coeffs[2:end], f.exps[2:end], length(f) - 1) + end + end + return FreeAssAlgElem{T}(R, rcoeffs, rexps, length(rcoeffs)) +end + +# normal form with leftmost word divisions +function normal_form( + f::FreeAssAlgElem{T}, + g::Vector{FreeAssAlgElem{T}}, +) where T<:FieldElement + R = parent(f) + s = length(g) + rcoeffs = T[] + rexps = Monomial[] + while length(f) > 0 + i = 1 + @label again + ok, ml, mr = word_divides_leftmost(f.exps[1], g[i].exps[1]) + if !ok && i < s + i += 1 + @goto again + end + if ok + qi = divexact(f.coeffs[1], g[i].coeffs[1]) + f = _sub_rest(f, mul_term(qi, ml, g[i], mr), 1) # enforce lt cancelation + else + push!(rcoeffs, f.coeffs[1]) + push!(rexps, f.exps[1]) + f = FreeAssAlgElem{T}(R, f.coeffs[2:end], f.exps[2:end], length(f) - 1) + end + end + r = FreeAssAlgElem{T}(R, rcoeffs, rexps, length(rcoeffs)) + return r +end + +# weak normal form with leftmost word divisions +function normal_form_weak( + f::FreeAssAlgElem{T}, + g::Vector{FreeAssAlgElem{T}}, +) where T<:FieldElement + R = parent(f) + s = length(g) + while length(f) > 0 + i = 1 + @label again + ok, ml, mr = word_divides_leftmost(f.exps[1], g[i].exps[1]) + if !ok && i < s + i += 1 + @goto again + end + if ok + qi = divexact(f.coeffs[1], g[i].coeffs[1]) + f = _sub_rest(f, mul_term(qi, ml, g[i], mr), 1) # enforce lt cancelation + else + break + end + end + return f +end + +@doc raw""" + interreduce!(g::Vector{FreeAssAlgElem{T}}) where T + +Interreduce a given Groebner basis with itself, i.e. compute the normal form of each +element of `g` with respect to the rest of the elements and discard elements with +normal form $0$ and duplicates. +""" +function interreduce!(g::Vector{FreeAssAlgElem{T}}) where T + i = 1 + while length(g) > 1 && length(g) >= i + aut = AhoCorasickAutomaton([g_j.exps[1] for g_j in g[1:end .!= i]]) + r = normal_form(g[i], g[1:end .!= i], aut) + if iszero(r) + deleteat!(g, i) + elseif g[i] != r + g[i] = r + i = 1 + else + i += 1 + end + end + return g +end + +## checks whether there is an overlap between a and b at position i of b +# such that b[i:length(b)] = a[1:length(b)-i] +function check_left_overlap(a::Monomial, b::Monomial, i::Int) + if length(b) - i >= length(a) + return false # this is a not a left overlap but might be a center overlap + end + for j in 0:(length(b) - i) + if b[i + j] != a[j + 1] + return false + end + end + return true +end + +### +# find all non-trivial left-obstructions of a and b +# i.e. all words w_1 and w_2^' s.t. w_1 a = b w_2^' +# where length(w_1) < length(b) and length(w_2^') < length(a) +# the return vector is of the form [(w_1, w_2^'), ...] +# if w_1 or w_2^' is empty, the corresponding obstruction is not returned +function left_obstructions(a::Monomial, b::Monomial) + v = Tuple{Monomial,Monomial}[] + for i in 2:length(b) + if check_left_overlap(a, b, i) + if length(b) - i + 2 <= length(a) # w_2^' should not be empty! + push!(v, (b[1:(i - 1)], a[(length(b) - i + 2):length(a)])) + end + end + end + return v +end + + +### +# find all non-trivial right-obstructions of a and b +# i.e. all words w_2 and w_1^' s.t. a w_1^' = w_2 b +# where length(w_1^') < length(b) and length(w_2) < length(a) +# the return vector is of the form [(w_1^', w_2), ...] +# if w_1^' or w_2 is empty, the corresponding obstruction is not returned +function right_obstructions(a::Monomial, b::Monomial) + left_obstr = left_obstructions(b, a) + right_obstr = [(w1, w2) for (w2, w1) in left_obstr] + return right_obstr +end + +### +# check whether a is a subword of b starting at index i +# a == b is also allowed +function check_center_overlap(a::Vector{Int}, b::Vector{Int}, i::Int) + i + length(a) - 1 <= length(b) || return false + return all(j -> a[j] == b[i + j - 1], 1:length(a)) +end + + +function center_obstructions_first_in_second(a::Monomial, b::Monomial) + v = Tuple{Monomial,Monomial}[] + for i in 1:length(b)-length(a) + 1 + if check_center_overlap(a, b, i) + push!(v, (b[1:(i - 1)], b[(i + length(a)):length(b)])) + end + end + return v +end + +## +# return all center obstructions of a and b, i.e. all (w_i, w_i^') +# such that either +# w_i a w_i^' = b +# or +# w_i b w_i^' = a +# either or both of w_i and w_i^' can be empty +function center_obstructions(a::Monomial, b::Monomial) + if length(a) > length(b) + return center_obstructions_first_in_second(b, a) + else + return center_obstructions_first_in_second(a, b) + end +end + +# all non-trivial ones +function obstructions(a::Monomial, b::Monomial) + one = Int[] # the empty word + res = Tuple{Monomial,Monomial,Monomial,Monomial}[] + for x in center_obstructions_first_in_second(b, a) + push!(res, (one, one, x[1], x[2])) + end + for x in center_obstructions_first_in_second(a, b) + push!(res, (x[1], x[2], one, one)) + end + for x in left_obstructions(a, b) + push!(res, (x[1], one, one, x[2])) + end + for x in left_obstructions(b, a) + push!(res, (one, x[2], x[1], one)) + end + for x in res + @assert vcat(x[1], a, x[2]) == vcat(x[3], b, x[4]) + end + return res +end + +# all non-trivial self obstructions +function obstructions(a::Monomial) + one = Int[] # the empty word + res = Tuple{Monomial,Monomial,Monomial,Monomial}[] + for x in left_obstructions(a, a) + push!(res, (one, x[2], x[1], one)) + end + for x in res + @assert vcat(x[1], a, x[2]) == vcat(x[3], a, x[4]) + end + return res +end + + +# check whether w_2 = v w_1 for some word v +function is_subword_right(w_1::Monomial, w_2::Monomial) + if length(w_1) > length(w_2) + return false + end + for i in 0:(length(w_1) - 1) + if w_1[length(w_1) - i] != w_2[length(w_2) - i] + return false + end + end + return true +end + + +# check whether w_2 = w_1 v for some word v +function is_subword_left(w_1::Monomial, w_2::Monomial) + if length(w_1) > length(w_2) + return false + end + for i in 1:length(w_1) + if w_1[i] != w_2[i] + return false + end + end + return true +end + + +### +# check if obs2 is a subobstructon of obs1, i.e. if +# the second pre-and suffix of obs2 are right- respectively left subwords of the second pre-and suffix of obs1. +# In other words, check if for obs1 = (w_i, w_i^'; u_j, u_j^') and obs2 = (w_k, w_k^'; v_l, v_l^') +# it holds that u_j == w v_l and u_j^' = v_l^' w^' for some w, w^' +# both w and w^' might be empty +function is_subobstruction(obs1_second_prefix::Monomial, obs1_second_suffix::Monomial, + obs2_second_prefix::Monomial, obs2_second_suffix) + return is_subword_right(obs2_second_prefix, obs1_second_prefix) && is_subword_left(obs2_second_suffix, obs1_second_suffix) +end + +function is_subobstruction(obs1::ObstructionTriple{T}, obs2::ObstructionTriple{T}) where T + return is_subobstruction(obs1.second_prefix, obs1.second_suffix, obs2.second_prefix, obs2.second_suffix) +end + +""" +if obs2 is a subobstruction of obs1, i.e. obs1[3] = w obs2[3] and obs1[4] = obs2[4]w', +returns length(ww') +thus, if it returns 0 and obs2 is a subobstruction of obs1, they are equal +if obs2 is not a subobstruction of obs1 the return value is useless +""" +function get_diff_length_for_subobstruction( + obs1::ObstructionTriple{T}, + obs2::ObstructionTriple{T}, +) where T + return length(obs1.second_prefix) - length(obs2.second_prefix) + + length(obs1.second_suffix) - length(obs2.second_suffix) +end + +# check whether there exists a (possibly empty) w^'' such that +# w1 LM(g1) w2 = w1 LM(g1) w^'' LM(g2) u2 +# only returns the correct value, if all the arguments come from an obstruction +# and if that is the case, i.e. there is no overlap, returns false +# assumes that (w1, w2; u1, u2) are an obstruction of g1 and g2 +# i.e. w1 LM(g1) w2 = u1 LM(g2) u2 +function has_overlap(g1, g2, w1, w2, u1, u2) + @assert vcat(w1, _leading_word(g1), w2) == vcat(u1, _leading_word(g2), u2) + lw2 = _leading_word(g2) + return length(w2) < length(lw2) + length(u2) +end + +function has_overlap(g2, w2, u2) + lw2 = _leading_word(g2) + return length(w2) < length(lw2) + length(u2) +end + +function has_overlap(obs::ObstructionTriple{T}) where T + return has_overlap(obs.second_poly, obs.first_suffix, obs.second_suffix) +end + +function is_redundant( + obs::ObstructionTriple{T}, + new_obstructions::PriorityQueue{Obstruction{T},FreeAssAlgElem{T}}, +) where T + # case 4b from Thm. 4.2.22 in Non-Commutative Groebner Bases and Applications by Xingqiang Xiu + for obstruction_pair in new_obstructions + o = obstruction_pair[1] + if o.second_index == obs.second_index + if is_subobstruction(obs, o) + if get_diff_length_for_subobstruction(obs, o) > 0 + return true + elseif obs.first_index > o.first_index + return true + elseif obs.first_index == o.first_index && + get_diff_length_for_subobstruction(obs, o) == 0 && + word_gt(obs.first_prefix, o.first_prefix) + return true + end + end + end + end + return false +end + +""" +check, whether obs1 is a proper multiple of obs2, i.e. they belong to the same polynomials and are of the form +obs1 = [w w_i, w_i' w'; w w_j, w_j' w'] and obs2 = [w_i, w_i'; w_j, w_j'] +""" +function is_proper_multiple( + obs1::ObstructionTriple{T}, + obs2::ObstructionTriple{T}, +) where T + if obs1.first_poly != obs2.first_poly || obs1.second_poly != obs2.second_poly #TODO compare indices instead? + return false + end + if is_subword_right(obs2.first_prefix, obs1.first_prefix) && + is_subword_left(obs2.first_suffix, obs1.first_suffix) + w = copy(obs1.first_prefix) + w2 = copy(obs1.first_suffix) + for _ in 1:length(obs2.first_prefix) + pop!(w) + end + for _ in 1:length(obs2.first_suffix) + popfirst!(w2) + end + if length(w) + length(w2) == 0 + return false + end + @assert obs1.first_prefix == vcat(w, obs2.first_prefix) + @assert obs1.first_suffix == vcat(obs2.first_suffix, w2) + return obs1.second_prefix == vcat(w, obs2.second_prefix) && + obs1.second_suffix == vcat(obs2.second_suffix, w2) + else + return false + end +end + +""" +check, whether obs is a proper multiple of any of the obstructions in the priority queue +""" +function is_proper_multiple( + obs::ObstructionTriple{T}, + obstructions::PriorityQueue{Obstruction{T},FreeAssAlgElem{T}}, +) where T + for obspair in obstructions + obs2 = obspair[1] + if is_proper_multiple(obs, obs2) + return true + end + end + return false +end + +function is_redundant( + obs::ObstructionTriple{T}, + new_obstructions::PriorityQueue{Obstruction{T},FreeAssAlgElem{T}}, + newest_element::FreeAssAlgElem{T}, + newest_index::Int, +) where T + w1 = [] + w2 = [] + for i in 1:length(obs.second_poly.exps[1]) + word_to_check = + vcat(obs.second_prefix, obs.second_poly.exps[1], obs.second_suffix) + if check_center_overlap(newest_element.exps[1], word_to_check, i) + w1 = word_to_check[1:(i - 1)] + w2 = word_to_check[(i + length(newest_element.exps[1])):end] + break + end + end + if length(w1) + length(w2) == 0 + return false + end + obs1 = ObstructionTriple{T}( + obs.first_poly, + newest_element, + obs.first_prefix, obs.first_suffix, w1, w2, + obs.first_index, + newest_index, + ) + obs2 = ObstructionTriple{T}( + obs.second_poly, + newest_element, + obs.second_prefix, obs.second_suffix, w1, w2, + obs.second_index, + newest_index, + ) + o1_bool = !has_overlap(obs1) || is_proper_multiple(obs1, new_obstructions) # TODO maybe only call is_proper_multiple if both obs have no overlap for performance? + o2_bool = !has_overlap(obs2) || is_proper_multiple(obs2, new_obstructions) + return o1_bool && o2_bool +end + + +function remove_redundancies!( + all_obstructions::PriorityQueue{Obstruction{T},FreeAssAlgElem{T}}, + newest_index::Int, + newest_element::FreeAssAlgElem{T}, +) where T + del_counter = 0 + new_obstructions = PriorityQueue{Obstruction{T},FreeAssAlgElem{T}}() + old_obstructions = PriorityQueue{Obstruction{T},FreeAssAlgElem{T}}() + for obstr_pair in all_obstructions + if obstr_pair[1].second_index == newest_index + new_obstructions[obstr_pair[1]] = obstr_pair[2] + else + old_obstructions[obstr_pair[1]] = obstr_pair[2] + end + end + for obstr_pair in new_obstructions + if is_redundant(obstr_pair[1], new_obstructions) + del_counter += 1 + delete!(new_obstructions, obstr_pair[1]) + delete!(all_obstructions, obstr_pair[1]) + end + end + if length(new_obstructions) == 0 + return nothing + end + + for obstr_pair in old_obstructions + if is_redundant(obstr_pair[1], new_obstructions, newest_element, newest_index) + del_counter += 1 + delete!(all_obstructions, obstr_pair[1]) + end + end + # TODO case 4e from Thm 4.1 in Kreuzer Xiu +end + +function get_obstructions(g::Vector{FreeAssAlgElem{T}}) where T + s = length(g) + result = PriorityQueue{Obstruction{T},FreeAssAlgElem{T}}() + for i in 1:s, j in 1:i + if i == j + obs = obstructions(_leading_word(g[i])) + else + obs = obstructions(_leading_word(g[i]), _leading_word(g[j])) + end + for o in obs + triple = ObstructionTriple{T}(g[i], g[j], o, i, j) + push!(result, triple => common_multiple_leading_term(triple)) + end + end + # TODO maybe here some redundancies can be removed too, check Kreuzer Xiu + return result +end + + +function add_obstructions!( + obstruction_queue::PriorityQueue{Obstruction{T},FreeAssAlgElem{T}}, + g::Vector{FreeAssAlgElem{T}}, +) where T + s = length(g) + for i in 1:s + if i == s + obs = obstructions(_leading_word(g[i])) + else + obs = obstructions(_leading_word(g[i]), _leading_word(g[s])) + end + for o in obs + triple = ObstructionTriple{T}(g[i], g[s], o, i, s) + push!(obstruction_queue, triple => common_multiple_leading_term(triple)) + end + end + #remove_redundancies!(obstruction_queue, s, g[s]) #TODO too slow in practice +end + + +function groebner_basis_buchberger( + g::Vector{FreeAssAlgElem{T}}, + reduction_bound::Int = typemax(Int), + remove_redundancies::Bool = false +) where T<:FieldElement + g = copy(g) + # interreduce!(g) # on some small examples, this increases running time, so it might not be optimal to use this here + nonzero_reductions = 0 + # compute the aho corasick automaton + # to make normal form computation more efficient + aut = AhoCorasickAutomaton([g_i.exps[1] for g_i in g]) + # step 1 from Thm. 5.2.12 Noncommutative Groebner Bases and Applications, Xingqiang Xiu + obstruction_queue = get_obstructions(g) + while !isempty(obstruction_queue) # step 2 + obstruction = popfirst!(obstruction_queue)[1] + # step3 + S = s_polynomial(obstruction) + Sp = normal_form(S, g, aut) # or normal_form_weak + if iszero(Sp) + continue + end + nonzero_reductions += 1 + # step 4 + push!(g, Sp) + insert_keyword!(aut, Sp.exps[1], length(g)) + if nonzero_reductions >= reduction_bound + return g + end + add_obstructions!(obstruction_queue, g) + if remove_redundancies + remove_redundancies!(obstruction_queue, length(g), g[length(g)]) + end + end + return g +end + +@doc """ + groebner_basis(g::Vector{FreeAssAlgElem{T}}, reduction_bound::Int = typemax(Int), remove_redundancies::Bool = false) + +Compute a Groebner basis for the ideal generated by `g`. Stop when `reduction_bound` many +non-zero entries have been added to the Groebner basis. If the computation stops due to the bound being exceeded, +the result is in general not an actual Groebner basis, just a subset of one. However, whenever the normal form with +respect to this incomplete Groebner basis is `0`, it will also be `0` with respect to the full Groebner basis. +""" +function groebner_basis( + g::Vector{FreeAssAlgElem{T}}, + reduction_bound::Int = typemax(Int), + remove_redundancies::Bool = false +) where T<:FieldElement + return groebner_basis_buchberger(g, reduction_bound, remove_redundancies) +end diff --git a/src/generic/PriorityQueue.jl b/src/generic/PriorityQueue.jl new file mode 100644 index 0000000000..131294352e --- /dev/null +++ b/src/generic/PriorityQueue.jl @@ -0,0 +1,439 @@ +# The contents of this file were taken from the Package DataStructures.jl: https://juliacollections.github.io/DataStructures.jl/latest/ + +# This file contains code that was formerly a part of Julia. License is MIT: http://julialang.org/license + +# PriorityQueue +# ------------- + +export PriorityQueue + +using Base: Ordering, ForwardOrdering, Forward, ReverseOrdering, Reverse, lt + +# Binary heap indexing +heapleft(i::Integer) = 2i +heapright(i::Integer) = 2i + 1 +heapparent(i::Integer) = div(i, 2) + +""" + PriorityQueue{K, V}([ord]) + +Construct a new `PriorityQueue`, with keys of type `K` and values/priorities +of type `V`. If an order is not given, the priority queue is min-ordered using +the default comparison for `V`. + +A `PriorityQueue` acts like a `Dict`, mapping values to their +priorities. New elements are added using `push!` and retrieved +using `popfirst!` or `popat!` based on their priority. + +Parameters +--------- + +`K::Type` Data type for the keys + +`V::Type` Data type for the values/priorities + +`ord::Base.Ordering` Priority queue ordering + +# Examples +```jldoctest; setup = :(using AbstractAlgebra) +julia> Generic.PriorityQueue(Base.Order.Forward, "a" => 2, "b" => 3, "c" => 1) +AbstractAlgebra.Generic.PriorityQueue{String, Int64, Base.Order.ForwardOrdering} with 3 entries: + "c" => 1 + "a" => 2 + "b" => 3 +``` +""" +struct PriorityQueue{K,V,O<:Ordering} <: AbstractDict{K,V} + # Binary heap of (element, priority) pairs. + xs::Vector{Pair{K,V}} + o::O + + # Map elements to their index in xs + index::Dict{K, Int} + + function PriorityQueue{K,V,O}(o::O) where {K,V,O<:Ordering} + new{K,V,O}(Vector{Pair{K,V}}(), o, Dict{K, Int}()) + end + + PriorityQueue{K, V, O}(xs::Vector{Pair{K,V}}, o::O, index::Dict{K, Int}) where {K,V,O<:Ordering} = new(xs, o, index) + + function PriorityQueue{K,V,O}(o::O, itr) where {K,V,O<:Ordering} + xs = Vector{Pair{K,V}}(undef, length(itr)) + index = Dict{K, Int}() + for (i, (k, v)) in enumerate(itr) + xs[i] = Pair{K,V}(k, v) + if haskey(index, k) + throw(ArgumentError("PriorityQueue keys must be unique")) + end + index[k] = i + end + pq = new{K,V,O}(xs, o, index) + + # heapify + for i in heapparent(length(pq.xs)):-1:1 + percolate_down!(pq, i) + end + + return pq + end +end + +# A copy constructor +PriorityQueue(xs::Vector{Pair{K,V}}, o::O, index::Dict{K, Int}) where {K,V,O<:Ordering} = + PriorityQueue{K,V,O}(xs, o, index) + +# Any-Any constructors +PriorityQueue(o::Ordering=Forward) = PriorityQueue{Any,Any,typeof(o)}(o) + +# Construction from Pairs +PriorityQueue(ps::Pair...) = PriorityQueue(Forward, ps) +PriorityQueue(o::Ordering, ps::Pair...) = PriorityQueue(o, ps) +PriorityQueue{K,V}(ps::Pair...) where {K,V} = PriorityQueue{K,V,ForwardOrdering}(Forward, ps) +PriorityQueue{K,V}(o::Ord, ps::Pair...) where {K,V,Ord<:Ordering} = PriorityQueue{K,V,Ord}(o, ps) + +# Construction specifying Key/Value types +# e.g., PriorityQueue{Int,Float64}([1=>1, 2=>2.0]) +PriorityQueue{K,V}(kv) where {K,V} = PriorityQueue{K,V}(Forward, kv) +function PriorityQueue{K,V}(o::Ord, kv) where {K,V,Ord<:Ordering} + try + PriorityQueue{K,V,Ord}(o, kv) + catch e + if not_iterator_of_pairs(kv) + throw(ArgumentError("PriorityQueue(kv): kv needs to be an iterator of tuples or pairs")) + else + rethrow(e) + end + end +end + +# Construction inferring Key/Value types from input +# e.g. PriorityQueue{} + +PriorityQueue(o1::Ordering, o2::Ordering) = throw(ArgumentError("PriorityQueue with two parameters must be called with an Ordering and an iterable of pairs")) +PriorityQueue(kv, o::Ordering=Forward) = PriorityQueue(o, kv) +function PriorityQueue(o::Ordering, kv) + try + _priority_queue_with_eltype(o, kv, eltype(kv)) + catch e + if not_iterator_of_pairs(kv) + throw(ArgumentError("PriorityQueue(kv): kv needs to be an iterator of tuples or pairs")) + else + rethrow(e) + end + end +end + +_priority_queue_with_eltype(o::Ord, ps, ::Type{Pair{K,V}} ) where {K,V,Ord} = PriorityQueue{ K, V,Ord}(o, ps) +_priority_queue_with_eltype(o::Ord, kv, ::Type{Tuple{K,V}}) where {K,V,Ord} = PriorityQueue{ K, V,Ord}(o, kv) +_priority_queue_with_eltype(o::Ord, ps, ::Type{Pair{K}} ) where {K, Ord} = PriorityQueue{ K,Any,Ord}(o, ps) +_priority_queue_with_eltype(o::Ord, kv, ::Type ) where { Ord} = PriorityQueue{Any,Any,Ord}(o, kv) + +## TODO: It seems impossible (or at least very challenging) to create the eltype below. +## If deemed possible, please create a test and uncomment this definition. +# _priority_queue_with_eltype{ D,Ord}(o::Ord, ps, ::Type{Pair{K,V} where K}) = PriorityQueue{Any, D,Ord}(o, ps) + + +""" + length(pq::PriorityQueue) + +Return the number of pairs (`k`, `v`) in the priority queue `pq`. +""" +Base.length(pq::PriorityQueue) = length(pq.xs) + +""" + isempty(pq::PriorityQueue) + +Verify if priority queue `pq` is empty. +""" +Base.isempty(pq::PriorityQueue) = isempty(pq.xs) + +""" + haskey(pq::PriorityQueue, key) + +Verify if priority queue `pq` has `key` in its keys. + +# Example + +```jldoctest; setup = :(using AbstractAlgebra) +julia> pq = Generic.PriorityQueue("a" => 1, "b" => 2, "c" => 3) +AbstractAlgebra.Generic.PriorityQueue{String, Int64, Base.Order.ForwardOrdering} with 3 entries: + "a" => 1 + "b" => 2 + "c" => 3 + +julia> haskey(pq, "a") +true + +julia> haskey(pq, "e") +false +``` +""" +Base.haskey(pq::PriorityQueue, key) = haskey(pq.index, key) + +""" + first(pq::PriorityQueue) + +Return the lowest priority pair (`k`, `v`) from `pq` without removing it from the +priority queue. +""" +Base.first(pq::PriorityQueue) = first(pq.xs) + +function percolate_down!(pq::PriorityQueue, i::Integer) + x = pq.xs[i] + @inbounds while (l = heapleft(i)) <= length(pq) + r = heapright(i) + j = r > length(pq) || lt(pq.o, pq.xs[l].second, pq.xs[r].second) ? l : r + xj = pq.xs[j] + if lt(pq.o, xj.second, x.second) + pq.index[xj.first] = i + pq.xs[i] = xj + i = j + else + break + end + end + pq.index[x.first] = i + pq.xs[i] = x +end + + +function percolate_up!(pq::PriorityQueue, i::Integer) + x = pq.xs[i] + @inbounds while i > 1 + j = heapparent(i) + xj = pq.xs[j] + if lt(pq.o, x.second, xj.second) + pq.index[xj.first] = i + pq.xs[i] = xj + i = j + else + break + end + end + pq.index[x.first] = i + pq.xs[i] = x +end + +# Equivalent to percolate_up! with an element having lower priority than any other +function force_up!(pq::PriorityQueue, i::Integer) + x = pq.xs[i] + @inbounds while i > 1 + j = heapparent(i) + pq.index[pq.xs[j].first] = i + pq.xs[i] = pq.xs[j] + i = j + end + pq.index[x.first] = i + pq.xs[i] = x +end + +Base.getindex(pq::PriorityQueue, key) = pq.xs[pq.index[key]].second + +function Base.get(pq::PriorityQueue, key, default) + i = get(pq.index, key, 0) + i == 0 ? default : pq.xs[i].second +end + +function Base.get!(pq::PriorityQueue, key, default) + i = get(pq.index, key, 0) + if i == 0 + push!(pq, key=>default) + return default + else + return pq.xs[i].second + end +end + +# Change the priority of an existing element, or enqueue it if it isn't present. +function Base.setindex!(pq::PriorityQueue{K, V}, value, key) where {K,V} + i = get(pq.index, key, 0) + if i != 0 + @inbounds oldvalue = pq.xs[i].second + pq.xs[i] = Pair{K,V}(key, value) + if lt(pq.o, oldvalue, value) + percolate_down!(pq, i) + else + percolate_up!(pq, i) + end + else + push!(pq, key=>value) + end + return value +end + +""" + push!(pq::PriorityQueue{K,V}, pair::Pair{K,V}) where {K,V} + +Insert the a key `k` into a priority queue `pq` with priority `v`. + +# Examples + +```jldoctest; setup = :(using AbstractAlgebra) +julia> a = Generic.PriorityQueue("a" => 1, "b" => 2, "c" => 3, "e" => 5) +AbstractAlgebra.Generic.PriorityQueue{String, Int64, Base.Order.ForwardOrdering} with 4 entries: + "a" => 1 + "b" => 2 + "c" => 3 + "e" => 5 + +julia> push!(a, "d" => 4) +AbstractAlgebra.Generic.PriorityQueue{String, Int64, Base.Order.ForwardOrdering} with 5 entries: + "a" => 1 + "b" => 2 + "c" => 3 + "d" => 4 + "e" => 5 +``` +""" +function Base.push!(pq::PriorityQueue{K,V}, pair::Pair{K,V}) where {K,V} + key = pair.first + if haskey(pq, key) + throw(ArgumentError("PriorityQueue keys must be unique")) + end + push!(pq.xs, pair) + pq.index[key] = length(pq) + percolate_up!(pq, length(pq)) + + return pq +end + +Base.push!(pq::PriorityQueue{K,V}, kv::Pair) where {K,V} = push!(pq, Pair{K,V}(kv.first, kv.second)) + +""" + popfirst!(pq::PriorityQueue) + +Remove and return the lowest priority key and value from a priority queue `pq` as a pair. + +# Examples + +```jldoctest; setup = :(using AbstractAlgebra) +julia> a = Generic.PriorityQueue(Base.Order.Forward, "a" => 2, "b" => 3, "c" => 1) +AbstractAlgebra.Generic.PriorityQueue{String, Int64, Base.Order.ForwardOrdering} with 3 entries: + "c" => 1 + "a" => 2 + "b" => 3 + +julia> popfirst!(a) +"c" => 1 + +julia> a +AbstractAlgebra.Generic.PriorityQueue{String, Int64, Base.Order.ForwardOrdering} with 2 entries: + "a" => 2 + "b" => 3 +``` +""" +function Base.popfirst!(pq::PriorityQueue) + x = pq.xs[1] + y = pop!(pq.xs) + if !isempty(pq) + @inbounds pq.xs[1] = y + pq.index[y.first] = 1 + percolate_down!(pq, 1) + end + delete!(pq.index, x.first) + return x +end + +if isdefined(Base, :popat!) # We will overload if it is defined, else we define on our own + import Base: popat! +end + +function popat!(pq::PriorityQueue, key) + idx = pq.index[key] + force_up!(pq, idx) + popfirst!(pq) +end + +""" + delete!(pq::PriorityQueue, key) + +Delete the mapping for the given `key` in a priority queue `pq` and return the priority queue. + +# Examples + +```jldoctest; setup = :(using AbstractAlgebra) +julia> q = Generic.PriorityQueue(Base.Order.Forward, "a" => 2, "b" => 3, "c" => 1) +AbstractAlgebra.Generic.PriorityQueue{String, Int64, Base.Order.ForwardOrdering} with 3 entries: + "c" => 1 + "a" => 2 + "b" => 3 +julia> delete!(q, "b") +AbstractAlgebra.Generic.PriorityQueue{String, Int64, Base.Order.ForwardOrdering} with 2 entries: + "c" => 1 + "a" => 2 +``` +""" +function Base.delete!(pq::PriorityQueue, key) + popat!(pq, key) + return pq +end + +""" + empty!(pq::PriorityQueue) + +Reset priority queue `pq`. +""" +function Base.empty!(pq::PriorityQueue) + empty!(pq.xs) + empty!(pq.index) + return pq +end + +Base.empty(pq::PriorityQueue) = PriorityQueue(empty(pq.xs), pq.o, empty(pq.index)) + +#Base.merge!(d::SortedDict, other::PriorityQueue) = invoke(merge!, Tuple{AbstractDict, PriorityQueue}, d, other) + +function Base.merge!(d::AbstractDict, other::PriorityQueue) + next = iterate(other, false) + while next !== nothing + (k, v), state = next + d[k] = v + next = iterate(other, state) + end + return d +end + +function Base.merge!(combine::Function, d::AbstractDict, other::PriorityQueue) + next = iterate(other, false) + while next !== nothing + (k, v), state = next + d[k] = haskey(d, k) ? combine(d[k], v) : v + next = iterate(other, state) + end + return d +end + +# Opaque not to be exported. +mutable struct _PQIteratorState{K, V, O <: Ordering} + pq::PriorityQueue{K, V, O} + _PQIteratorState{K, V, O}(pq::PriorityQueue{K, V, O}) where {K, V, O <: Ordering} = new(pq) +end + +_PQIteratorState(pq::PriorityQueue{K, V, O}) where {K, V, O <: Ordering} = _PQIteratorState{K, V, O}(pq) + +# Unordered iteration through key value pairs in a PriorityQueue +# O(n) iteration. +function _iterate(pq::PriorityQueue, state) + (k, idx), i = state + return (pq.xs[idx], i) +end +_iterate(pq::PriorityQueue, ::Nothing) = nothing + +Base.iterate(pq::PriorityQueue, ::Nothing) = nothing + +function Base.iterate(pq::PriorityQueue, ordered::Bool=true) + if ordered + isempty(pq) && return nothing + state = _PQIteratorState(PriorityQueue(copy(pq.xs), pq.o, copy(pq.index))) + return popfirst!(state.pq), state + else + _iterate(pq, iterate(pq.index)) + end +end + +function Base.iterate(pq::PriorityQueue, state::_PQIteratorState) + isempty(state.pq) && return nothing + return popfirst!(state.pq), state +end + +Base.iterate(pq::PriorityQueue, i) = _iterate(pq, iterate(pq.index, i)) diff --git a/test/NCRings-test.jl b/test/NCRings-test.jl index 9107297539..a4e4b4666c 100644 --- a/test/NCRings-test.jl +++ b/test/NCRings-test.jl @@ -1,6 +1,7 @@ include("generic/MatrixAlgebra-test.jl") include("generic/NCPoly-test.jl") include("generic/FreeAssAlgebra-test.jl") +include("generic/FreeAssAlgebraGroebner-test.jl") @testset "NCRings.oftype" begin F = GF(3) diff --git a/test/generic/AhoCorasick-test.jl b/test/generic/AhoCorasick-test.jl new file mode 100644 index 0000000000..70f45e0b7c --- /dev/null +++ b/test/generic/AhoCorasick-test.jl @@ -0,0 +1,16 @@ +using AbstractAlgebra.Generic: AhoCorasickAutomaton, search, AhoCorasickMatch, aho_corasick_automaton +@testset "Generic.AhoCorasick" begin + keywords = [[1, 2, 3, 4], [1, 5, 4], [4, 1, 2], [1, 2]] + aut = aho_corasick_automaton(keywords) + @test search(aut, [10, 4, 1, 2, 3, 4]) == AhoCorasickMatch(6, 1, [1, 2, 3, 4]) + @test hash(search(aut, [10, 4, 1, 2, 3, 4])) == hash(AhoCorasickMatch(6, 1, [1, 2, 3, 4])) + @test isnothing(search(aut, Int[])) + @test search(aut, [1, 5, 4, 1, 1, 1, 4, 4]) == AhoCorasickMatch(3, 2, [1, 5, 4]) + @test search(aut, [1, 2, 3, 1, 4, 1, 2, 1, 4, 1, 2]) == AhoCorasickMatch(7, 3, [4, 1, 2]) + @test search(aut, [2, 1, 2, 3, 1]) == AhoCorasickMatch(3, 4, [1, 2]) + @test isnothing(search(aut, [1, 3, 1, 5, 1, 4, 8])) + @test isnothing(search(aut, [8, 8, 7, 10, 456])) + @test search(aut, [4, 1, 5, 4]) == AhoCorasickMatch(4, 2, [1, 5, 4]) + @test isnothing(search(aut, [4, 1, 5, 10])) + @test !isnothing(AhoCorasickAutomaton(Vector{Int}[])) +end diff --git a/test/generic/FreeAssAlgebraGroebner-test.jl b/test/generic/FreeAssAlgebraGroebner-test.jl new file mode 100644 index 0000000000..3e89daf20f --- /dev/null +++ b/test/generic/FreeAssAlgebraGroebner-test.jl @@ -0,0 +1,87 @@ +include("AhoCorasick-test.jl") +using AbstractAlgebra.Generic: AhoCorasickAutomaton +import AbstractAlgebra.Generic: normal_form_weak +@testset "Generic.FreeAssAlgebra.groebner" begin + + R, (x, y, u, v, t, s) = FreeAssociativeAlgebra(GF(2), ["x", "y", "u", "v", "t", "s"]) + g = AbstractAlgebra.groebner_basis([u*(x*y)^3 + u*(x*y)^2 + u + v, + (y*x)^3*t + (y*x)^2*t + t + s]) + @test length(g) >= 5 + + # Example 6.1 Kreuzer & Xiu + R, (a, b) = FreeAssociativeAlgebra(QQ, ["a", "b"]) + + g = AbstractAlgebra.groebner_basis([a^2 - 1, b^3 - 1, (a*b*a*b^2)^2 - 1]) + AbstractAlgebra.interreduce!(g) + @test length(g) == 5 + + g = AbstractAlgebra.groebner_basis([a^2 - 1, b^3 - 1, (a*b*a*b*a*b^2)^2 - 1]) + AbstractAlgebra.interreduce!(g) + @test length(g) == 15 + + g2 = AbstractAlgebra.groebner_basis([a^2 - 1, b^3 - 1, (a*b*a*b*a*b^2)^2 - 1], typemax(Int), true) + @test all(u ->iszero(normal_form(u, g2)), g) # make sure removing redundant obstructions in the computation does not change the groebner basis + @test all(u ->iszero(normal_form(u, g)), g2) + +end + +@testset "Generic.FreeAssociativeAlgebra.groebner.normal_form" begin + R, (x, y, u, v, t, s) = FreeAssociativeAlgebra(QQ, ["x", "y", "u", "v", "t", "s"]) # x > y > ... > s + ideal_generators = [x*y, u*y*t, s*t - t*s, x*y*y + x*x*y - one(R)] + aut = AhoCorasickAutomaton([g_i.exps[1] for g_i in ideal_generators]) + @test normal_form(x*y, ideal_generators, aut) == zero(R) + @test normal_form(u*y*t, ideal_generators, aut) == zero(R) + @test normal_form(s*t - t*s, ideal_generators, aut) == zero(R) + @test x*y*y + x*x*y - one(R) in ideal_generators + @test normal_form(x*y*t + y, ideal_generators, aut) == y + @test normal_form(one(R), ideal_generators, aut) == one(R) + @test normal_form(v*y, ideal_generators, aut) == v*y + @test normal_form(x*y, ideal_generators) == zero(R) + @test normal_form(u*y*t, ideal_generators) == zero(R) + @test normal_form(s*t - t*s, ideal_generators) == zero(R) + @test normal_form(x*y*t + y, ideal_generators) == y + @test normal_form(one(R), ideal_generators) == one(R) + @test normal_form(v*y, ideal_generators) == v*y + @test normal_form(x*y*v*v + t*s*x*y*v + y*s*t - y*t*s + v*x*y*y*s + v*x*x*y*s, ideal_generators, aut) == zero(R) + @test normal_form_weak(x*y*v*v + t*s*x*y*v + y*s*t - y*t*s + v*x*y*y*s + v*x*x*y*s, ideal_generators) == zero(R) + @test normal_form_weak(x*y + u*y, ideal_generators) <= x*y + u*y + +# @test gb_divides_leftmost((x*y*u*v*t).exps[1], aut) == (true, [], [3, 4, 5], 1) +# @test gb_divides_leftmost((x*u*y*t*(s*t - t*s)).exps[1], aut) == (true, [1], [5, 6], 2) +# @test gb_divides_leftmost((x*s*t).exps[1], aut) == (false, [], [], -1) +end + +@testset "Generic.FreeAssociativeAlgebra.groebner.overlaps_and_obstructions" begin + w1 = [1, 1, 2, 1, 3] + w2 = [2, 1, 3, 4, 3, 4] + w3 = [1, 3, 4] + w4 = [1, 1, 2, 1] + w5 = [5, 1, 3, 4, 4, 2, 1] + @test AbstractAlgebra.Generic.check_left_overlap(w2, w1, 3) + @test !AbstractAlgebra.Generic.check_left_overlap(w1, w2, 3) + @test AbstractAlgebra.Generic.check_center_overlap(w3, w2, 2) + @test AbstractAlgebra.Generic.check_center_overlap(w4, w1, 1) + @test !AbstractAlgebra.Generic.check_left_overlap(w4, w1, 1) + R, (x, y, z) = FreeAssociativeAlgebra(QQ, ["x", "y", "z"]) + poly1 = x*y*x*x*z + poly2 = y*x*x*z*y*y + poly3 = x*y*x*y + poly4 = y*x*y*z + poly5 = x*y + poly6 = x*y*x*y*z*x*y + lw1 = AbstractAlgebra.Generic._leading_word(poly1) + lw2 = AbstractAlgebra.Generic._leading_word(poly2) + lw3 = AbstractAlgebra.Generic._leading_word(poly3) + lw4 = AbstractAlgebra.Generic._leading_word(poly4) + lw5 = AbstractAlgebra.Generic._leading_word(poly5) + lw6 = AbstractAlgebra.Generic._leading_word(poly6) + ot = AbstractAlgebra.Generic.ObstructionTriple{Rational{BigInt}}( poly1, poly2, [], [2, 2], [1], [], 1, 2) + + @test AbstractAlgebra.Generic.has_overlap(ot) + @test length(AbstractAlgebra.Generic.left_obstructions(lw2, lw1)) == 1 + @test isempty(AbstractAlgebra.Generic.left_obstructions(lw1, lw2)) + @test length(AbstractAlgebra.Generic.right_obstructions(lw3, lw4)) == 2 + @test isempty(AbstractAlgebra.Generic.left_obstructions(lw3, lw4)) + @test length(AbstractAlgebra.Generic.center_obstructions(lw5, lw6)) == 3 + @test length(AbstractAlgebra.Generic.get_obstructions([poly1, poly2])) == 2 +end