Js/nc groebner optimization stable

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
+ ```

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

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

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

src/generic/FreeAssAlgebra.jl

@@ -370,6 +370,62 @@
     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 @@
     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})