From fa8eff49ff8039d613ae0cf93374bf48f332c375 Mon Sep 17 00:00:00 2001 From: Sasha Demin Date: Thu, 19 Sep 2024 01:51:38 +0200 Subject: [PATCH 1/6] clean up a bit --- src/Groebner.jl | 10 +- src/arithmetic/Zp.jl | 88 ----------- src/f4/basis.jl | 314 +++++++++++-------------------------- src/f4/f4.jl | 81 +++------- src/f4/hashtable.jl | 54 +------ src/f4/learn_apply.jl | 17 +- src/f4/linalg/backend.jl | 59 ------- src/f4/matrix.jl | 1 - src/f4/trace.jl | 5 - src/groebner/groebner.jl | 12 -- src/groebner/parameters.jl | 41 +---- src/interface.jl | 21 --- src/utils/keywords.jl | 34 +--- test/groebner/groebner.jl | 21 +-- test/runtests.jl | 1 - 15 files changed, 126 insertions(+), 633 deletions(-) diff --git a/src/Groebner.jl b/src/Groebner.jl index aa448309..b1717774 100644 --- a/src/Groebner.jl +++ b/src/Groebner.jl @@ -27,14 +27,6 @@ It is useful to enable this when debugging the Groebner package. """ invariants_enabled() = false -""" - logging_enabled() -> Bool - -Specifies if logging is enabled. If `false`, then all logging in Groebner is -disabled, and entails **(almost)** no runtime overhead. -""" -logging_enabled() = true - ### # Imports @@ -159,7 +151,7 @@ include("groebner/homogenization.jl") include("interface.jl") using PrecompileTools -include("precompile.jl") +# include("precompile.jl") ### # Exports diff --git a/src/arithmetic/Zp.jl b/src/arithmetic/Zp.jl index 6ac4f61b..2501f14f 100644 --- a/src/arithmetic/Zp.jl +++ b/src/arithmetic/Zp.jl @@ -309,84 +309,6 @@ function inv_mod_p(a::T, arithm::SignedCompositeArithmeticZp{T}) where {T} CompositeNumber(invmod.(a.data, arithm.ps.data)) end -### -# FloatingPointArithmeticZp - -struct FloatingPointArithmeticZp{AccumType, CoeffType} <: AbstractArithmeticZp{AccumType, CoeffType} - multiplier::AccumType - divisor::AccumType - - function FloatingPointArithmeticZp( - ::Type{AccumType}, - ::Type{CoeffType}, - p::CoeffType - ) where {AccumType <: CoeffZp, CoeffType <: CoeffZp} - @invariant AccumType === Float64 - @invariant 0 < p < 2^25 # < 52 / 2 - 1 - @invariant Primes.isprime(Int(p)) - multiplier = 1 / p - new{AccumType, CoeffType}(multiplier, p) - end -end - -divisor(arithm::FloatingPointArithmeticZp) = arithm.divisor - -@inline function mod_p(a::T, mod::FloatingPointArithmeticZp{T, C}) where {T, C} - b = a * mod.multiplier - c = floor(b) - a - mod.divisor * c # may be fused -end - -@inline function fma_mod_p(a1::T, a2::C, a3::T, mod::FloatingPointArithmeticZp{T, C}) where {T, C} - b = muladd(a1, a2, a3) - mod_p(b, mod) -end - -inv_mod_p(a::T, arithm::FloatingPointArithmeticZp{T}) where {T} = - T(invmod(Int(a), Int(divisor(arithm)))) - -### -# FloatingPointCompositeArithmeticZp - -struct FloatingPointCompositeArithmeticZp{AccumType, CoeffType, T, N} <: - AbstractArithmeticZp{AccumType, CoeffType} - multiplier::CompositeNumber{N, T} - divisor::CompositeNumber{N, T} - - function FloatingPointCompositeArithmeticZp( - ::Type{CompositeNumber{N, AT}}, - ::Type{CompositeNumber{N, CT}}, - ps::CompositeNumber{N, CT} - ) where {N, AT <: CoeffZp, CT <: CoeffZp} - @invariant AT === Float64 - @invariant all(0 .< ps.data .< 2^25) # < 52 / 2 - 1 - @invariant all(Primes.isprime.(Int.(ps.data))) - multiplier = inv(ps) - new{CompositeNumber{N, AT}, CompositeNumber{N, CT}, AT, N}(multiplier, ps) - end -end - -divisor(arithm::FloatingPointCompositeArithmeticZp) = arithm.divisor - -@inline function mod_p(a::T, mod::FloatingPointCompositeArithmeticZp{T, C}) where {T, C} - b = a * mod.multiplier - c = T(floor.(b.data)) - a - mod.divisor * c # may be fused -end - -@inline function fma_mod_p( - a1::T, - a2::C, - a3::T, - mod::FloatingPointCompositeArithmeticZp{T, C} -) where {T, C} - b = CompositeNumber(muladd.(a1.data, a2.data, a3.data)) - mod_p(b, mod) -end - -inv_mod_p(a::T, arithm::FloatingPointCompositeArithmeticZp{T}) where {T} = - T(invmod.(Int.(a.data), Int.(divisor(arithm).data))) - ### # Selection of arithmetic @@ -418,21 +340,11 @@ function select_arithmetic( if CoeffType <: CompositeCoeffZp if hint === :signed || CoeffType <: CompositeNumber{N, T} where {N, T <: Signed} return SignedCompositeArithmeticZp(AccumType, CoeffType, CoeffType(characteristic)) - elseif hint === :floating - return FloatingPointCompositeArithmeticZp( - AccumType, - CoeffType, - CoeffType(characteristic) - ) else return CompositeArithmeticZp(AccumType, CoeffType, CoeffType(characteristic)) end end - if hint === :floating - return FloatingPointArithmeticZp(AccumType, CoeffType, CoeffType(characteristic)) - end - if hint === :signed return SignedArithmeticZp(AccumType, CoeffType, CoeffType(characteristic)) end diff --git a/src/f4/basis.jl b/src/f4/basis.jl index 75fdcaa4..e73d263d 100644 --- a/src/f4/basis.jl +++ b/src/f4/basis.jl @@ -15,6 +15,11 @@ struct CriticalPair poly2::Int32 # Index of lcm(lm(poly1), lm(poly2)) in the hashtable lcm::MonomId + + function CriticalPair(poly1, poly2, lcm) + @invariant poly1 < poly2 + new(poly1, poly2, lcm) + end end mutable struct Pairset{ExponentType <: Integer} @@ -73,6 +78,9 @@ function pairset_find_smallest_degree_pair(ps::Pairset) pair_idx, pair_min_deg end +# Returns N, the number of critical pairs of the smallest degree. Sorts the +# critical pairs so that the first N pairs in the pairset are the smallest with +# respect to degree. function pairset_partition_by_degree!(ps::Pairset) @invariant ps.load > 0 _, pair_min_deg = pairset_find_smallest_degree_pair(ps) @@ -97,15 +105,6 @@ function pairset_partition_by_degree!(ps::Pairset) i - 1 end -# Returns N, the number of critical pairs of the smallest degree. Sorts the -# critical pairs so that the first N pairs in the pairset are the smallest with -# respect to degree. -function pairset_lowest_degree_pairs!(pairset::Pairset) - n_lowest_degree_pairs = pairset_partition_by_degree!(pairset) - @invariant n_lowest_degree_pairs > 0 - n_lowest_degree_pairs -end - ### # Basis @@ -190,7 +189,6 @@ function basis_changematrix_initialize!( basis.changematrix[i] = Dict{Int, Dict{MonomId, C}}() basis.changematrix[i][i] = Dict{MonomId, C}(id_of_1 => one(C)) end - nothing end function basis_changematrix_mul!( @@ -209,7 +207,6 @@ function basis_changematrix_mul!( new_row[pos] = new_poly end basis.changematrix[idx] = new_row - nothing end function basis_changematrix_addmul!( @@ -241,7 +238,6 @@ function basis_changematrix_addmul!( poly[new_monom_id] = mod_p(poly[new_monom_id] + ref_cf * AccumType(cf), arithmetic) end end - nothing end function basis_changematrix_deep_copy_with_new_type( @@ -402,28 +398,6 @@ function basis_make_monic!( basis end -function basis_make_monic!( - basis::Basis{C}, - arithmetic::Union{FloatingPointCompositeArithmeticZp{A, C}, FloatingPointArithmeticZp{A, C}}, - changematrix::Bool -) where {A <: Union{CoeffZp, CompositeCoeffZp}, C <: Union{CoeffZp, CompositeCoeffZp}} - cfs = basis.coeffs - @inbounds for i in 1:(basis.n_filled) - !isassigned(cfs, i) && continue - isone(cfs[i][1]) && continue - mul = inv_mod_p(A(cfs[i][1]), arithmetic) - cfs[i][1] = one(C) - for j in 2:length(cfs[i]) - cfs[i][j] = mod_p(A(cfs[i][j]) * A(mul), arithmetic) - end - @invariant isone(cfs[i][1]) - if changematrix - basis_changematrix_mul!(basis, i, A(mul), arithmetic) - end - end - basis -end - function basis_make_monic!( basis::Basis{C}, arithmetic::AbstractArithmeticQQ, @@ -442,9 +416,8 @@ function basis_make_monic!( basis end -# Generate new S-pairs from pairs of polynomials -# (basis[idx], basis[i]) -# for every i < idx +# Follows paragraph 4.4. from +# Gebauer, Möller, On an Installation of Buchberger's Algorithm function pairset_update!( pairset::Pairset{D}, basis::Basis{C}, @@ -452,100 +425,111 @@ function pairset_update!( update_ht::MonomialHashtable{M}, idx::Int ) where {D, C <: Coeff, M <: Monom} - pl, bl = pairset.load, idx - ps = pairset.pairs + pairset_resize_lcms_if_needed!(pairset, basis.n_filled) + pairset_resize_if_needed!(pairset, basis.n_filled) + + pairs = pairset.pairs lcms = pairset.lcms degs = pairset.degrees - @inbounds new_lead = basis.monoms[idx][1] + @inbounds lead_idx = basis.monoms[idx][1] # Generate new pairs. - hashtable_resize_if_needed!(update_ht, bl) - @inbounds for i in 1:(bl - 1) - newidx = pl + i + hashtable_resize_if_needed!(update_ht, idx) + @inbounds for i in 1:(idx - 1) + newidx = pairset.load + i if !basis.is_redundant[i] && - !monom_is_gcd_const(ht.monoms[basis.monoms[i][1]], ht.monoms[new_lead]) - lcms[i] = hashtable_get_lcm!(basis.monoms[i][1], new_lead, ht, update_ht) + !monom_is_gcd_const(ht.monoms[basis.monoms[i][1]], ht.monoms[lead_idx]) + lcms[i] = hashtable_get_lcm!(basis.monoms[i][1], lead_idx, ht, update_ht) degs[newidx] = monom_totaldeg(update_ht.monoms[lcms[i]]) - ps[newidx] = CriticalPair(Int32(i), Int32(idx), lcms[i]) + pairs[newidx] = CriticalPair(Int32(i), Int32(idx), lcms[i]) else lcms[i] = CRITICAL_PAIR_REDUNDANT degs[newidx] = typemax(D) - ps[newidx] = CriticalPair(Int32(i), Int32(idx), CRITICAL_PAIR_REDUNDANT) + pairs[newidx] = CriticalPair(Int32(i), Int32(idx), CRITICAL_PAIR_REDUNDANT) end end - # Traverse existing pairs... - @inbounds for i in 1:pl - if ps[i].lcm == CRITICAL_PAIR_REDUNDANT - continue - end - - j = ps[i].poly1 - l = ps[i].poly2 - m = max(degs[pl + l], degs[pl + j]) - - # ...if an existing pair is divisible by the lead of the new poly and - # has a greater degree than newly generated critical pair, then mark the - # existing pair redundant - if degs[i] > m && hashtable_monom_is_divisible(ps[i].lcm, new_lead, ht) - ps[i] = CriticalPair(ps[i].poly1, ps[i].poly2, CRITICAL_PAIR_REDUNDANT) + # Criterion B_k(j, l). + # T(j, l) is divisible by T(k) + # max(deg T(j, k), deg T(l, k)) < deg T(j, l) + @inbounds for i in 1:pairset.load + (pairs[i].lcm == CRITICAL_PAIR_REDUNDANT) && continue + j = pairs[i].poly1 + l = pairs[i].poly2 + m = max(degs[pairset.load + l], degs[pairset.load + j]) + if degs[i] > m && hashtable_monom_is_divisible(pairs[i].lcm, lead_idx, ht) + pairs[i] = CriticalPair(pairs[i].poly1, pairs[i].poly2, CRITICAL_PAIR_REDUNDANT) end end # Traverse new pairs to move non-redundant ones first. - j = 1 - @inbounds for i in 1:(bl - 1) - if !basis.is_redundant[i] - ps[pl + j] = ps[pl + i] - degs[pl + j] = degs[pl + i] - j += 1 - end - end - - sort_pairset_by_degree!(pairset, pl + 1, j - 2) - - @inbounds for i in 1:(j - 1) - lcms[i] = ps[pl + i].lcm - end - @inbounds lcms[j] = CRITICAL_PAIR_REDUNDANT - pc = j - pc -= 1 - - # Mark redundancy of some pairs based on their lcms. - @inbounds for j in 1:pc - if !(lcms[j] == CRITICAL_PAIR_REDUNDANT) - hashtable_check_monomial_division_in_update(lcms, j + 1, pc, lcms[j], update_ht) + cnt = 1 + @inbounds for i in 1:(idx - 1) + (lcms[i] == CRITICAL_PAIR_REDUNDANT) && continue + # if !basis.is_redundant[i] + pairs[pairset.load + cnt] = pairs[pairset.load + i] + degs[pairset.load + cnt] = degs[pairset.load + i] + cnt += 1 + # end + end + + sort_pairset_by_degree!(pairset, pairset.load + 1, cnt - 2) + + @inbounds for i in 1:(cnt - 1) + lcms[i] = pairs[pairset.load + i].lcm + end + @inbounds lcms[cnt] = CRITICAL_PAIR_REDUNDANT + new_pairs = cnt - 1 + + # Criterion M(i, j). + # T(i, j) | T(k, j) + @inbounds for i in 1:new_pairs + (lcms[i] == CRITICAL_PAIR_REDUNDANT) && continue + divmask_i = update_ht.divmasks[lcms[i]] + lcm_i = update_ht.monoms[lcms[i]] + k = i + 1 + for k in (i+1):new_pairs + (lcms[k] == CRITICAL_PAIR_REDUNDANT) && continue + if update_ht.use_divmask && !divmask_is_probably_divisible(update_ht.divmasks[lcms[k]], divmask_i) + continue + end + ea = update_ht.monoms[lcms[k]] + if monom_is_divisible(update_ht.monoms[lcms[k]], lcm_i) + lcms[k] = CRITICAL_PAIR_REDUNDANT + end end end # Remove redundant pairs from the pairset. - j = 1 + cnt = 1 @inbounds for i in 1:(pairset.load) - (ps[i].lcm == CRITICAL_PAIR_REDUNDANT) && continue - ps[j] = ps[i] - degs[j] = degs[i] - j += 1 + (pairs[i].lcm == CRITICAL_PAIR_REDUNDANT) && continue + pairs[cnt] = pairs[i] + degs[cnt] = degs[i] + cnt += 1 end - hashtable_resize_if_needed!(ht, pc) - - # Add new lcm monomials to the basis hashtable - # (including index j and not including index pc). - insert_lcms_in_basis_hashtable!(pairset, pl, ht, update_ht, basis, lcms, j, pc + 1) + # Add new lcm monomials to the main hashtable. + hashtable_resize_if_needed!(ht, new_pairs) + @inbounds for i in 1:new_pairs + (lcms[i] == CRITICAL_PAIR_REDUNDANT) && continue + id = hashtable_insert!(ht, update_ht.monoms[lcms[i]]) + pair = pairs[pairset.load + i] + pairs[cnt] = CriticalPair(pair.poly1, pair.poly2, id) + degs[cnt] = degs[pairset.load + i] + cnt += 1 + end + pairset.load = cnt - 1 - # Mark redundant polynomials in basis. + # Mark redundant polynomials in the generating set. nonred = basis.nonredundant_indices - lml = basis.n_nonredundant - @inbounds for i in 1:lml - if !basis.is_redundant[nonred[i]] - if hashtable_monom_is_divisible(basis.monoms[nonred[i]][1], new_lead, ht) - basis.is_redundant[nonred[i]] = true - end + @inbounds for i in 1:basis.n_nonredundant + basis.is_redundant[nonred[i]] && continue + if hashtable_monom_is_divisible(basis.monoms[nonred[i]][1], lead_idx, ht) + basis.is_redundant[nonred[i]] = true end end - - nothing end function basis_update!(basis::Basis, ht::MonomialHashtable{M}) where {M <: Monom} @@ -573,24 +557,21 @@ function basis_update!(basis::Basis, ht::MonomialHashtable{M}) where {M <: Monom basis.n_processed = basis.n_filled end -function basis_is_new_polynomial_redundant!( +function basis_mark_redundant_elements!( pairset::Pairset, basis::Basis, ht::MonomialHashtable{M}, update_ht::MonomialHashtable{M}, idx::Int ) where {M <: Monom} - hashtable_resize_if_needed!(update_ht, 0) - - @inbounds lead_new = basis.monoms[idx][1] + @inbounds lead_idx = basis.monoms[idx][1] ps = pairset.pairs degs = pairset.degrees @inbounds for i in (idx + 1):(basis.n_filled) basis.is_redundant[i] && continue - lead_i = basis.monoms[i][1] - @invariant !monom_isless(ht.monoms[lead_i], ht.monoms[lead_new], ht.ord) - !hashtable_monom_is_divisible(lead_i, lead_new, ht) && continue + @invariant !monom_isless(ht.monoms[lead_i], ht.monoms[lead_idx], ht.ord) + !hashtable_monom_is_divisible(lead_i, lead_idx, ht) && continue # Add a new critical pair corresponding to Spoly(i, idx). pairset_resize_if_needed!(pairset, 1) @@ -600,8 +581,6 @@ function basis_is_new_polynomial_redundant!( basis.is_redundant[i] = true end - - false end function basis_fill_data!( @@ -627,26 +606,7 @@ function basis_fill_data!( basis.n_filled = ngens end -function basis_sweep_redundant!(basis::Basis, hashtable::MonomialHashtable) - # here -- assert that basis is in fact a Groebner basis. - # NOTE: maybe sort generators for more effective sweeping? - @inbounds for i in 1:(basis.n_processed) - for j in (i + 1):(basis.n_processed) - basis.is_redundant[i] && continue - basis.is_redundant[j] && continue - lead_i = basis.monoms[i][1] - lead_j = basis.monoms[j][1] - if hashtable_monom_is_divisible(lead_i, lead_j, hashtable) - basis.is_redundant[i] = true - elseif hashtable_monom_is_divisible(lead_j, lead_i, hashtable) - basis.is_redundant[j] = true - end - end - end - nothing -end - -function basis_mark_redundant_elements!(basis::Basis) +function basis_move_redundant_elements!(basis::Basis) j = 1 @inbounds for i in 1:(basis.n_nonredundant) if !basis.is_redundant[basis.nonredundant_indices[i]] @@ -692,8 +652,7 @@ end function basis_get_monoms_by_identifiers(basis::Basis, ht::MonomialHashtable{M}) where {M <: Monom} monoms = Vector{Vector{M}}(undef, basis.n_nonredundant) @inbounds for i in 1:(basis.n_nonredundant) - idx = basis.nonredundant_indices[i] - poly = basis.monoms[idx] + poly = basis.monoms[basis.nonredundant_indices[i]] monoms[i] = Vector{M}(undef, length(poly)) for j in 1:length(poly) monoms[i][j] = ht.monoms[poly[j]] @@ -705,8 +664,7 @@ end function basis_export_coeffs(basis::Basis{C}) where {C <: Coeff} coeffs = Vector{Vector{C}}(undef, basis.n_nonredundant) @inbounds for i in 1:(basis.n_nonredundant) - idx = basis.nonredundant_indices[i] - coeffs[i] = basis.coeffs[idx] + coeffs[i] = basis.coeffs[basis.nonredundant_indices[i]] end coeffs end @@ -716,87 +674,3 @@ function basis_export_data(basis::Basis{C}, ht::MonomialHashtable{M}) where {M < coeffs = basis_export_coeffs(basis) exps, coeffs end - -# For a given list of S-pairs and a list of indices `plcm` -# adds indices from plcm[ifirst:ilast] -# to the hashtable ht -function insert_lcms_in_basis_hashtable!( - pairset::Pairset, - off::Int, - ht::MonomialHashtable{M}, - update_ht::MonomialHashtable{M}, - basis::Basis, - plcm::Vector{MonomId}, - ifirst::Int, - ilast::Int -) where {M} - # including ifirst and not including ilast - - monoms = basis.monoms - ps = pairset.pairs - degs = pairset.degrees - - mod = MonomHash(ht.size - 1) - @invariant ispow2(mod + 1) - - m = ifirst - l = 1 - @label Letsgo - @inbounds while l < ilast - if plcm[l] == CRITICAL_PAIR_REDUNDANT - l += 1 - continue - end - - if hashtable_monom_is_gcd_const( - monoms[ps[off + l].poly1][1], - monoms[ps[off + 1].poly2][1], - ht - ) - l += 1 - continue - end - - ps[m] = ps[off + l] - degs[m] = degs[off + l] - - h = update_ht.hashvals[plcm[l]] - ht.monoms[ht.load + 1] = monom_copy(update_ht.monoms[plcm[l]]) - n = ht.monoms[ht.load + 1] - - k = h - i = MonomHash(0) - @inbounds while i <= ht.size - k = hashtable_next_lookup_index(h, i, mod) - hm = ht.hashtable[k] - - # if free - iszero(hm) && break - - if hashtable_is_hash_collision(ht, hm, n, h) - i += MonomHash(1) - continue - end - - ps[m] = CriticalPair(ps[m].poly1, ps[m].poly2, hm) - m += 1 - l += 1 - @goto Letsgo - end - - @invariant !ht.frozen - - ht.hashtable[k] = pos = ht.load + 1 - - ht.labels[ht.load + 1] = NON_PIVOT_COLUMN - ht.hashvals[ht.load + 1] = h - ht.divmasks[ht.load + 1] = update_ht.divmasks[plcm[l]] - - ht.load += 1 - ps[m] = CriticalPair(ps[m].poly1, ps[m].poly2, MonomId(pos)) - m += 1 - l += 1 - end - - pairset.load = m - 1 -end diff --git a/src/f4/f4.jl b/src/f4/f4.jl index 1fc697fe..0f9625f8 100644 --- a/src/f4/f4.jl +++ b/src/f4/f4.jl @@ -28,8 +28,8 @@ function f4_initialize_structs( end if sort_input - # The sorting of input polynomials is not deterministic across different - # Julia versions when sorting only w.r.t. the leading term. + # The order of polynomials is not deterministic across different Julia + # versions when sorting only w.r.t. the leading term. permutation = sort_polys_by_lead_increasing!(basis, hashtable, params.changematrix) else permutation = collect(1:(basis.n_filled)) @@ -69,30 +69,25 @@ function f4_update!( ht::MonomialHashtable, update_ht::MonomialHashtable ) - @invariant basis.n_filled >= basis.n_processed - @inbounds for i in (basis.n_processed + 1):(basis.n_filled) - basis_is_new_polynomial_redundant!(pairset, basis, ht, update_ht, i) && continue - pairset_resize_lcms_if_needed!(pairset, basis.n_filled) - pairset_resize_if_needed!(pairset, basis.n_filled) + @invariant basis.n_processed <= basis.n_filled + for i in (basis.n_processed + 1):(basis.n_filled) + basis_mark_redundant_elements!(pairset, basis, ht, update_ht, i) pairset_update!(pairset, basis, ht, update_ht, i) end basis_update!(basis, ht) end +# Monomials that represent the columns of the matrix are stored in the symbol_ht +# hashtable. +# We traverse monomials searching for a reducer for each monomial. The hashtable +# may grow as reducers are added to the matrix, and the loop accounts for that. function f4_symbolic_preprocessing!( basis::Basis, matrix::MacaulayMatrix, ht::MonomialHashtable, symbol_ht::MonomialHashtable ) - # Monomials that represent the columns of the matrix are stored in the - # symbol_ht hashtable. - ncols = matrix.ncols_left matrix_resize_upper_part_if_needed!(matrix, matrix.ncols_left + symbol_ht.load) - - # Traverse all monomials in symbol_ht and search for a polynomial reducer - # for each monomial. The hashtable grows as polynomials with new monomials - # are added to the matrix, and the loop accounts for that. i = symbol_ht.offset @inbounds while i <= symbol_ht.load if symbol_ht.labels[i] != NON_PIVOT_COLUMN @@ -106,8 +101,6 @@ function f4_symbolic_preprocessing!( f4_find_multiplied_reducer!(basis, matrix, ht, symbol_ht, MonomId(i)) i += 1 end - - nothing end function f4_autoreduce!( @@ -210,8 +203,6 @@ function f4_select_tobereduced!( basis.nonredundant_indices[i] = i basis.divmasks[i] = ht.divmasks[basis.monoms[i][1]] end - - nothing end function f4_find_lead_divisor_use_divmask(i, divmask, basis) @@ -297,49 +288,30 @@ function f4_find_multiplied_reducer!( end function f4_select_critical_pairs!( - pairset::Pairset{ExponentType}, + pairset::Pairset, basis::Basis, matrix::MacaulayMatrix, ht::MonomialHashtable, symbol_ht::MonomialHashtable; - maxpairs::Int=INT_INF, select_all::Bool=false -) where {ExponentType} - npairs::Int = pairset.load - if !select_all - npairs = pairset_lowest_degree_pairs!(pairset) +) + if select_all + npairs = pairset.load + else + npairs = pairset_partition_by_degree!(pairset) end - npairs = min(npairs, maxpairs) - @invariant npairs > 0 - ps = pairset.pairs - degs = pairset.degrees - @inbounds deg::ExponentType = degs[1] - sort_pairset_by_lcm!(pairset, npairs, ht) - # When there is a limit on the number of selected pairs, we still add pairs - # which have the same lcm as the selected ones. - if npairs > maxpairs - navailable = npairs - npairs = maxpairs - lastlcm = ps[npairs].lcm - while npairs < navailable && ps[npairs + 1].lcm == lastlcm - npairs += 1 - end - end - f4_add_critical_pairs_to_matrix!(pairset, npairs, basis, matrix, ht, symbol_ht) - # Remove selected parirs from the pairset. + # Remove selected pairs from the pairset. @inbounds for i in 1:(pairset.load - npairs) - ps[i] = ps[i + npairs] - degs[i] = degs[i + npairs] + pairset.pairs[i] = pairset.pairs[i + npairs] + pairset.degs[i] = pairset.degs[i + npairs] end pairset.load -= npairs - - deg, npairs end function f4_add_critical_pairs_to_matrix!( @@ -446,14 +418,7 @@ function f4!( f4_update!(pairset, basis, hashtable, update_ht) while !isempty(pairset) - f4_select_critical_pairs!( - pairset, - basis, - matrix, - hashtable, - symbol_ht, - maxpairs=params.maxpairs - ) + f4_select_critical_pairs!(pairset, basis, matrix, hashtable, symbol_ht) f4_symbolic_preprocessing!(basis, matrix, hashtable, symbol_ht) @@ -465,12 +430,8 @@ function f4!( hashtable_reinitialize!(symbol_ht) end - if params.sweep - basis_sweep_redundant!(basis, hashtable) - end - - basis_mark_redundant_elements!(basis) - + basis_move_redundant_elements!(basis) + if params.reduced f4_autoreduce!(ring, basis, matrix, hashtable, symbol_ht, params) end diff --git a/src/f4/hashtable.jl b/src/f4/hashtable.jl index 18026cbb..bd263071 100644 --- a/src/f4/hashtable.jl +++ b/src/f4/hashtable.jl @@ -255,22 +255,14 @@ end ### # Insertion of monomials -# Returns the next look-up position in the table. # Must be within 1 <= ... <= mod+1 function hashtable_next_lookup_index(h::MonomHash, j::MonomHash, mod::MonomHash) ((h + j) & mod) + MonomHash(1) end -# if hash collision happened function hashtable_is_hash_collision(ht::MonomialHashtable, vidx, e, he) - # if not free and not same hash - @inbounds if ht.hashvals[vidx] != he - return true - end - # if not free and not same monomial - @inbounds if !monom_is_equal(ht.monoms[vidx], e) - return true - end + @inbounds ht.hashvals[vidx] != he && return true + @inbounds !monom_is_equal(ht.monoms[vidx], e) && return true false end @@ -436,12 +428,6 @@ function hashtable_monom_is_divisible(h1::MonomId, h2::MonomId, ht::MonomialHash monom_is_divisible(e1, e2) end -function hashtable_monom_is_gcd_const(h1::MonomId, h2::MonomId, ht::MonomialHashtable) - @inbounds e1 = ht.monoms[h1] - @inbounds e2 = ht.monoms[h2] - monom_is_gcd_const(e1, e2) -end - function hashtable_get_lcm!( he1::MonomId, he2::MonomId, @@ -458,42 +444,6 @@ function hashtable_get_lcm!( end ### -# What are those once again?.. - -# compare pairwise divisibility of lcms from a[first:last] with lcm -function hashtable_check_monomial_division_in_update( - a::Vector{MonomId}, - first::Int, - last::Int, - lcm::MonomId, - ht::MonomialHashtable{M} -) where {M <: Monom} - # Pairs are sorted w.r.t. lcm, we only need to check entries above the - # starting point. - - @inbounds divmask = ht.divmasks[lcm] - @inbounds lcmexp = ht.monoms[lcm] - - j = first - @inbounds while j <= last - if a[j] == CRITICAL_PAIR_REDUNDANT - j += 1 - continue - end - if ht.use_divmask && !divmask_is_probably_divisible(ht.divmasks[a[j]], divmask) - j += 1 - continue - end - ea = ht.monoms[a[j]] - if !monom_is_divisible(ea, lcmexp) - j += 1 - continue - end - a[j] = CRITICAL_PAIR_REDUNDANT - end - - nothing -end # Inserts a multiple of the polynomial into symbolic hashtable. # Writes the resulting monomial identifiers to the given row. diff --git a/src/f4/learn_apply.jl b/src/f4/learn_apply.jl index b0839913..38fb0e71 100644 --- a/src/f4/learn_apply.jl +++ b/src/f4/learn_apply.jl @@ -171,14 +171,7 @@ function f4_learn!( pairset_size = f4_update!(pairset, basis, hashtable, update_ht) while !isempty(pairset) - degree_i, npairs_i = f4_select_critical_pairs!( - pairset, - basis, - matrix, - hashtable, - symbol_ht, - maxpairs=params.maxpairs - ) + degree_i, npairs_i = f4_select_critical_pairs!(pairset, basis, matrix, hashtable, symbol_ht) push!(trace.critical_pair_sequence, (degree_i, npairs_i)) f4_symbolic_preprocessing!(basis, matrix, hashtable, symbol_ht) @@ -191,11 +184,7 @@ function f4_learn!( hashtable_reinitialize!(symbol_ht) end - if params.sweep - basis_sweep_redundant!(basis, hashtable) - end - - basis_mark_redundant_elements!(basis) + basis_move_redundant_elements!(basis) if params.reduced f4_reducegb_learn!(trace, ring, basis, matrix, hashtable, symbol_ht, params) @@ -523,8 +512,6 @@ function f4_apply!( hashtable_reinitialize!(symbol_ht) end - # basis_mark_redundant_elements!(basis) - if params.reduced symbol_ht = hashtable_initialize_secondary(hashtable) flag = f4_autoreduce_apply!( diff --git a/src/f4/linalg/backend.jl b/src/f4/linalg/backend.jl index cf7bc381..85e665c7 100644 --- a/src/f4/linalg/backend.jl +++ b/src/f4/linalg/backend.jl @@ -716,25 +716,6 @@ function linalg_row_make_monic!( pinv end -function linalg_row_make_monic!( - row::Vector{T}, - arithmetic::Union{FloatingPointCompositeArithmeticZp{A, T}, FloatingPointArithmeticZp{A, T}}, - first_nnz_index::Int=1 -) where {A <: Union{CoeffZp, CompositeCoeffZp}, T <: Union{CoeffZp, CompositeCoeffZp}} - @invariant !iszero(row[first_nnz_index]) - - lead = row[first_nnz_index] - isone(lead) && return lead - - @inbounds pinv = inv_mod_p(A(lead), arithmetic) - @inbounds row[first_nnz_index] = one(T) - @inbounds for i in (first_nnz_index + 1):length(row) - row[i] = mod_p(A(row[i]) * A(pinv), arithmetic) - end - - pinv -end - function linalg_row_make_monic!( row::Vector{T}, arithmetic::AbstractArithmeticQQ{T}, @@ -777,46 +758,6 @@ function linalg_vector_addmul_sparsedense_mod_p!( mul end -# Linear combination of dense vector and sparse vector modulo a prime. -function linalg_vector_addmul_sparsedense_mod_p!( - row::Vector{A}, - indices::Vector{I}, - coeffs::Vector{T}, - arithmetic::FloatingPointArithmeticZp -) where {I, A <: CoeffZp, T <: CoeffZp} - @invariant isone(coeffs[1]) - @invariant length(indices) == length(coeffs) - @invariant !isempty(indices) - - @inbounds mul = divisor(arithmetic) - row[indices[1]] - @fastmath @inbounds for j in 1:length(indices) - idx = indices[j] - row[idx] = fma_mod_p(mul, coeffs[j], row[idx], arithmetic) - end - - mul -end - -# Linear combination of dense vector and sparse vector modulo a prime. -function linalg_vector_addmul_sparsedense_mod_p!( - row::Vector{A}, - indices::Vector{I}, - coeffs::Vector{T}, - arithmetic::FloatingPointCompositeArithmeticZp -) where {I, A <: CompositeCoeffZp, T <: CompositeCoeffZp} - @invariant isone(coeffs[1]) - @invariant length(indices) == length(coeffs) - @invariant !isempty(indices) - - @inbounds mul = divisor(arithmetic) - row[indices[1]] - @fastmath @inbounds for j in 1:length(indices) - idx = indices[j] - row[idx] = fma_mod_p(mul, coeffs[j], row[idx], arithmetic) - end - - mul -end - # Linear combination of dense vector and sparse vector. # The most generic version. function linalg_vector_addmul_sparsedense!( diff --git a/src/f4/matrix.jl b/src/f4/matrix.jl index bf270462..5ed73d7d 100644 --- a/src/f4/matrix.jl +++ b/src/f4/matrix.jl @@ -147,7 +147,6 @@ function matrix_resize_lower_part!(matrix::MacaulayMatrix, size::Int) nothing end -# statistics_refresh and partially initialize the matrix function matrix_reinitialize!(matrix::MacaulayMatrix, size::Int) new_size = size * 2 matrix_resize_upper_part!(matrix, new_size) diff --git a/src/f4/trace.jl b/src/f4/trace.jl index 840d5fb7..22e22e47 100644 --- a/src/f4/trace.jl +++ b/src/f4/trace.jl @@ -37,7 +37,6 @@ mutable struct Trace{C1 <: Coeff, C2 <: Coeff, M <: Monom, Ord1, Ord2} output_sort_indices::Vector{Int} params::AlgorithmParameters - sweep_output::Bool representation::PolynomialRepresentation homogenize::Bool @@ -91,7 +90,6 @@ function trace_initialize( Vector{Int}(), Vector{Int}(), params, - params.sweep, PolynomialRepresentation(ExponentVector{UInt64}, UInt64, false), params.homogenize, 0, @@ -131,7 +129,6 @@ function trace_deepcopy( copy(trace.nonredundant_indices_before_reduce), copy(trace.output_sort_indices), deepcopy(trace.params), - trace.sweep_output, PolynomialRepresentation( trace.representation.monomtype, trace.representation.coefftype, @@ -209,7 +206,6 @@ function trace_copy( trace.nonredundant_indices_before_reduce, trace.output_sort_indices, trace.params, - trace.sweep_output, new_representation, trace.homogenize, trace.napply, @@ -345,7 +341,6 @@ function Base.show(io::IO, ::MIME"text/plain", trace::Trace) """ input order : $(trace.original_ord) output order : $(trace.ring.ord) - sweep : $(trace.sweep_output) homogenize : $(trace.homogenize) permute : $(permute_input) monom. type : $(trace.representation.monomtype) diff --git a/src/groebner/groebner.jl b/src/groebner/groebner.jl index 7cfb813c..bb84222a 100644 --- a/src/groebner/groebner.jl +++ b/src/groebner/groebner.jl @@ -1,17 +1,5 @@ # This file is a part of Groebner.jl. License is GNU GPL v2. -# Groebner.jl interface goes in four levels of access: -# high level (0), intermediate level (1), low level (2), and very low level (3). -# -# Level (0) works with frontend polynomials (say, AbstractAlgebra.jl). -# Level (1) works with exponent vectors and coefficients. -# Level (2) works with internal representations of monomials and coefficients. -# Level (3) works with internal F4 data structures (matrices, hashtables, etc). -# -# User interface is levels (0) and (1). Useful work is done on levels (2) and -# (3). F4 runs on level (3). A series of conversions carries data across levels. -# Generally, input and output of a function must live on the same level. - ### # Backend for `groebner` diff --git a/src/groebner/parameters.jl b/src/groebner/parameters.jl index 96408db2..4987637b 100644 --- a/src/groebner/parameters.jl +++ b/src/groebner/parameters.jl @@ -139,10 +139,6 @@ function param_select_coefftype( tight_signed_type = tight_signed_int_type(char) - if arithmetic === :floating - return Float64, true - end - if arithmetic === :signed if typemax(Int32) < char < typemax(UInt32) || typemax(Int64) < char < typemax(UInt64) @info "Cannot use $(arithmetic) arithmetic with characteristic $char" @@ -195,11 +191,6 @@ mutable struct AlgorithmParameters{Arithmetic <: AbstractArithmetic} # If reduced Groebner basis is needed reduced::Bool - # Limit the number of critical pairs in the F4 matrix by this number - maxpairs::Int - - selection_strategy::Symbol - # Strategy for modular computation in groebner. This can be one of the # following: # - :classic_modular @@ -209,26 +200,16 @@ mutable struct AlgorithmParameters{Arithmetic <: AbstractArithmetic} # If learn & apply strategy can use apply in batches batched::Bool - # In modular computation of the basis, compute (at least!) this many bases - # modulo different primes until a consensus in is reached + # In multi-modular computation, compute (at least!) this many bases modulo + # different primes until a consensus in is reached majority_threshold::Int # Multi-threading threaded_f4::Symbol threaded_multimodular::Symbol - # Random number generator rng::Random.Xoshiro - # Internal option. - # At the end of F4, polynomials are interreduced. - # We can mark and sweep polynomials that are redundant prior to - # interreduction to speed things up a bit. This option specifies if such - # sweep should be done. - sweep::Bool - - statistics::Symbol - changematrix::Bool end @@ -300,12 +281,6 @@ function AlgorithmParameters(ring::PolyRing, kwargs::KeywordArguments; hint=:non end reduced = kwargs.reduced - maxpairs = kwargs.maxpairs - - selection_strategy = kwargs.selection - if selection_strategy === :auto - selection_strategy = :normal - end threaded = kwargs.threaded if !(_threaded[]) @@ -346,10 +321,6 @@ function AlgorithmParameters(ring::PolyRing, kwargs::KeywordArguments; hint=:non seed = kwargs.seed rng = Random.Xoshiro(seed) - sweep = kwargs.sweep - - statistics = kwargs.statistics - changematrix = kwargs.changematrix if changematrix if !(target_ord isa DegRevLex) @@ -369,16 +340,12 @@ function AlgorithmParameters(ring::PolyRing, kwargs::KeywordArguments; hint=:non arithmetic, representation, reduced, - maxpairs, - selection_strategy, modular_strategy, batched, majority_threshold, threaded_f4, threaded_multimodular, rng, - sweep, - statistics, changematrix ) end @@ -410,16 +377,12 @@ function params_mod_p( select_arithmetic(C, prime, :auto, is_wide_type_coeffs), representation, params.reduced, - params.maxpairs, - params.selection_strategy, params.modular_strategy, params.batched, params.majority_threshold, params.threaded_f4, params.threaded_multimodular, params.rng, - params.sweep, - params.statistics, params.changematrix ) end diff --git a/src/interface.jl b/src/interface.jl index 473dd494..5fdbe1b2 100644 --- a/src/interface.jl +++ b/src/interface.jl @@ -58,20 +58,11 @@ The `groebner` routine takes the following optional arguments: - `loglevel`: Logging level, one of `:all`, `:debug`, `:info`, `:warn`, `:no`, or an integer. An integer `0` is equivalent to `:info`, and lower integer values mean more verbose. Default is `:info`. -- `maxpairs`: The maximum number of critical pairs used at once in matrix - construction in the F4 algorithm. By default, this is not specified. Tweak - this option to try to lower the amount of RAM consumed. - `homogenize`: Controls the use of homogenization in the algorithm. Possible options are: - `:yes`, always homogenize the input ideal, - `:no`, never homogenize the input ideal, - `:auto`, for the automatic choice (default). -- `statistics`: After the function exits, print some runtime statistics. - Possible options are: - - `:no`, print nothing (default), - - `:timings`, print the table with timings and allocations. Note that - `Groebner.performance_counters_enabled()` must be set to `true` for this to - have effect. ## Example @@ -450,12 +441,6 @@ The `isgroebner` routine takes the following optional arguments: - `loglevel`: Logging level, one of `:all`, `:debug`, `:info`, `:warn`, `:no`, or an integer. An integer `0` is equivalent to `:info`, and lower integer values mean more verbose. Default is `:info`. -- `statistics`: After the function exits, print some runtime statistics. - Possible options are: - - `:no`, do not print anything (default), - - `:timings`, print the table with timings and allocations. Note that - `Groebner.performance_counters_enabled()` must be set to `true` for this to - have effect. ## Example @@ -519,12 +504,6 @@ The `normalform` routine takes the following optional arguments: - `loglevel`: Logging level, one of `:all`, `:debug`, `:info`, `:warn`, `:no`, or an integer. An integer `0` is equivalent to `:info`, and lower integer values mean more verbose. Default is `:info`. -- `statistics`: After the function exits, print some runtime statistics. - Possible options are: - - `:no`, do not print anything (default), - - `:timings`, print the table with timings and allocations. Note that - `Groebner.performance_counters_enabled()` must be set to `true` for this to - have effect. ## Example diff --git a/src/utils/keywords.jl b/src/utils/keywords.jl index 1bdbd71e..33d4b85f 100644 --- a/src/utils/keywords.jl +++ b/src/utils/keywords.jl @@ -18,13 +18,10 @@ const _supported_kw_args = ( monoms = :auto, arithmetic = :auto, seed = 42, - maxpairs = INT_INF, selection = :auto, modular = :auto, threaded = :auto, - sweep = false, homogenize = :auto, - statistics = :no, batched = true, changematrix = false ), @@ -32,14 +29,12 @@ const _supported_kw_args = ( check = false, ordering = InputOrdering(), monoms = :dense, - statistics = :no ), isgroebner = ( ordering = InputOrdering(), certify = true, seed = 42, - monoms = :dense, - statistics = :no + monoms = :dense ), groebner_learn = ( seed = 42, @@ -47,8 +42,6 @@ const _supported_kw_args = ( monoms = :auto, arithmetic = :auto, homogenize = :auto, - sweep = true, - statistics = :no, threaded = :auto, ), groebner_apply! = ( @@ -56,8 +49,6 @@ const _supported_kw_args = ( ordering = InputOrdering(), monoms = :auto, arithmetic = :auto, - sweep = true, - statistics = :no, threaded = :auto, ), groebner_with_change_matrix = ( @@ -68,13 +59,10 @@ const _supported_kw_args = ( monoms = :auto, arithmetic = :auto, seed = 42, - maxpairs = INT_INF, selection = :auto, modular = :auto, threaded = :auto, - sweep = false, homogenize = :auto, - statistics = :no, batched = true, changematrix = true ), @@ -102,14 +90,11 @@ mutable struct KeywordArguments monoms::Symbol arithmetic::Symbol seed::Int - maxpairs::Int selection::Symbol modular::Symbol batched::Bool check::Bool - sweep::Bool homogenize::Symbol - statistics::Symbol changematrix::Bool KeywordArguments(function_id::Symbol; passthrough...) = @@ -157,16 +142,13 @@ mutable struct KeywordArguments `:auto`, `:dense`, `:packed`""" arithmetic = get(kws, :arithmetic, get(default_kw_args, :arithmetic, :auto)) - @assert arithmetic in (:auto, :delayed, :signed, :basic, :floating) """ + @assert arithmetic in (:auto, :delayed, :signed, :basic) """ Not recognized arithmetic: $arithmetic Possible choices for keyword "arithmetic" are: - `:auto`, `:delayed`, `:signed`, `:basic`, `:floating`""" + `:auto`, `:delayed`, `:signed`, `:basic`""" seed = get(kws, :seed, get(default_kw_args, :seed, 42)) - maxpairs = get(kws, :maxpairs, get(default_kw_args, :maxpairs, INT_INF)) - @assert maxpairs > 0 "The limit on the number of critical pairs must be positive" - modular = get(kws, :modular, get(default_kw_args, :modular, :auto)) @assert modular in (:auto, :classic_modular, :learn_and_apply) """ Not recognized modular strategy: $modular @@ -179,19 +161,12 @@ mutable struct KeywordArguments @assert selection in (:auto, :normal, :sugar, :be_divided_and_perish) check = get(kws, :check, get(default_kw_args, :check, true)) - sweep = get(kws, :sweep, get(default_kw_args, :sweep, false)) homogenize = get(kws, :homogenize, get(default_kw_args, :homogenize, :auto)) @assert homogenize in (:auto, :no, :yes) """ Not recognized homogenization strategy: $homogenize Possible choices for keyword "homogenize" are: `:auto`, `:no`, `:yes`""" - statistics = get(kws, :statistics, get(default_kw_args, :statistics, :no)) - @assert statistics in (:no, :timings, :stats, :all) """ - Not recognized option for collecting statistics: $statistics - Possible choices for keyword "statistics" are: - `:no`, `:timings`, `:stats`, `:all`""" - changematrix = get(kws, :changematrix, get(default_kw_args, :changematrix, false)) new( @@ -204,14 +179,11 @@ mutable struct KeywordArguments monoms, arithmetic, seed, - maxpairs, selection, modular, batched, check, - sweep, homogenize, - statistics, changematrix ) end diff --git a/test/groebner/groebner.jl b/test/groebner/groebner.jl index 5f5784b5..11233bce 100644 --- a/test/groebner/groebner.jl +++ b/test/groebner/groebner.jl @@ -490,25 +490,6 @@ end @test Groebner.groebner(system) == Groebner.groebner(system, certify=true) end -@testset "groebner maxpairs" begin - # TODO: fix some bugs in maxpairs - s = Groebner.Examples.noonn(5, internal_ordering=:degrevlex, k=GF(2^31 - 1)) - gb = Groebner.groebner(s) - # @test gb == Groebner.groebner(s, maxpairs=100) - # @test gb == Groebner.groebner(s, maxpairs=10) - # @test gb == Groebner.groebner(s, maxpairs=2) - # @test gb == Groebner.groebner(s, maxpairs=1) - @test_throws AssertionError Groebner.groebner(s, maxpairs=0) - - s = Groebner.Examples.katsuran(5, internal_ordering=:degrevlex, k=QQ) - gb = Groebner.groebner(s) - # @test gb == Groebner.groebner(s, maxpairs=100) - # @test gb == Groebner.groebner(s, maxpairs=10) - # @test gb == Groebner.groebner(s, maxpairs=2) - # @test gb == Groebner.groebner(s, maxpairs=1) - @test_throws AssertionError Groebner.groebner(s, maxpairs=0) -end - @testset "groebner orderings" begin R, (x, y, z, w) = polynomial_ring(QQ, ["x", "y", "z", "w"], internal_ordering=:deglex) @@ -815,7 +796,7 @@ end @testset "groebner arithmetic" begin R, (x, y, z) = polynomial_ring(GF(10007), ["x", "y", "z"], internal_ordering=:degrevlex) - for arithmetic in [:auto, :delayed, :signed, :basic, :floating] + for arithmetic in [:auto, :delayed, :signed, :basic] @test Groebner.groebner([x, y], arithmetic=arithmetic, linalg=:deterministic) == Groebner.groebner([y, x]) == [y, x] diff --git a/test/runtests.jl b/test/runtests.jl index f27620f2..cf63548d 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -9,7 +9,6 @@ using Groebner # Check invariants during testing. Groebner.invariants_enabled() = true -Groebner.logging_enabled() = true function is_github_ci() return parse(Bool, get(ENV, "GITHUB_ACTIONS", "false")) From a47676d39566ac3aa492d67e2eeba304e11bbf48 Mon Sep 17 00:00:00 2001 From: Sasha Demin Date: Thu, 19 Sep 2024 02:05:04 +0200 Subject: [PATCH 2/6] add a couple more assertions for basis --- src/f4/basis.jl | 4 +++- src/f4/f4.jl | 6 +++++- 2 files changed, 8 insertions(+), 2 deletions(-) diff --git a/src/f4/basis.jl b/src/f4/basis.jl index e73d263d..ff9dabfb 100644 --- a/src/f4/basis.jl +++ b/src/f4/basis.jl @@ -164,7 +164,9 @@ function basis_well_formed(ring::PolyRing, basis::Basis, hashtable::MonomialHash !is_sorted_by_lead_increasing(basis, hashtable) && error("Basis elements must be sorted") for i in basis.n_filled isempty(basis.monoms[i]) && error("Zero polynomials are not allowed.") - !(length(basis.monoms[i]) == length(basis.coeffs[i])) && error("Beda!") + !(length(basis.monoms[i]) == length(basis.coeffs[i])) && error("Bad polynomial") + !allunique(basis.monoms[i]) && error("Bad polynomial") + !issorted(basis.monoms[i], lt=(j,k) -> monom_isless(hashtable.monoms[k], hashtable.monoms[j], ring.ord)) && error("Bad polynomial") for j in 1:length(basis.coeffs[i]) iszero(basis.coeffs[i][j]) && error("Coefficient is zero") (ring.ch > 0) && diff --git a/src/f4/f4.jl b/src/f4/f4.jl index 0f9625f8..c209aa4c 100644 --- a/src/f4/f4.jl +++ b/src/f4/f4.jl @@ -306,12 +306,16 @@ function f4_select_critical_pairs!( f4_add_critical_pairs_to_matrix!(pairset, npairs, basis, matrix, ht, symbol_ht) + deg = pairset.degrees[1] + # Remove selected pairs from the pairset. @inbounds for i in 1:(pairset.load - npairs) pairset.pairs[i] = pairset.pairs[i + npairs] - pairset.degs[i] = pairset.degs[i + npairs] + pairset.degrees[i] = pairset.degrees[i + npairs] end pairset.load -= npairs + + deg, npairs end function f4_add_critical_pairs_to_matrix!( From 534a8f4fa2124c2e9b6d511904255a2840752fc7 Mon Sep 17 00:00:00 2001 From: Sasha Demin Date: Fri, 20 Sep 2024 22:58:46 +0200 Subject: [PATCH 3/6] guess lucky prime in multi-modular --- src/Groebner.jl | 3 +- src/f4/basis.jl | 24 +-- src/f4/f4.jl | 5 +- src/f4/learn_apply.jl | 31 ++-- src/groebner/groebner.jl | 189 ++++++++++---------- src/groebner/groebner_with_change_matrix.jl | 48 ++--- src/groebner/isgroebner.jl | 4 +- src/groebner/modular.jl | 94 +++++++--- src/groebner/primes.jl | 88 --------- src/reconstruction/crt.jl | 51 +----- src/reconstruction/ratrec.jl | 9 +- test/groebner/groebner.jl | 33 ++-- 12 files changed, 248 insertions(+), 331 deletions(-) delete mode 100644 src/groebner/primes.jl diff --git a/src/Groebner.jl b/src/Groebner.jl index b1717774..8f264cfb 100644 --- a/src/Groebner.jl +++ b/src/Groebner.jl @@ -33,7 +33,7 @@ invariants_enabled() = false # Groebner does not provide a polynomial implementation of its own but relies on # existing symbolic computation packages in Julia for communicating with the # user. Groebner accepts as its input polynomials from the Julia packages -# AbstractAlgebra.jl, Nemo.jl (Oscar.jl), and MultivariatePolynomials.jl. +# AbstractAlgebra.jl, Nemo.jl, and MultivariatePolynomials.jl. import AbstractAlgebra import AbstractAlgebra: base_ring, elem_type @@ -137,7 +137,6 @@ include("reconstruction/crt.jl") include("reconstruction/ratrec.jl") #= more high level functions =# -include("groebner/primes.jl") include("groebner/modular.jl") include("groebner/groebner.jl") include("groebner/groebner_with_change_matrix.jl") diff --git a/src/f4/basis.jl b/src/f4/basis.jl index ff9dabfb..c8e1628b 100644 --- a/src/f4/basis.jl +++ b/src/f4/basis.jl @@ -166,7 +166,10 @@ function basis_well_formed(ring::PolyRing, basis::Basis, hashtable::MonomialHash isempty(basis.monoms[i]) && error("Zero polynomials are not allowed.") !(length(basis.monoms[i]) == length(basis.coeffs[i])) && error("Bad polynomial") !allunique(basis.monoms[i]) && error("Bad polynomial") - !issorted(basis.monoms[i], lt=(j,k) -> monom_isless(hashtable.monoms[k], hashtable.monoms[j], ring.ord)) && error("Bad polynomial") + !issorted( + basis.monoms[i], + lt=(j, k) -> monom_isless(hashtable.monoms[k], hashtable.monoms[j], ring.ord) + ) && error("Bad polynomial") for j in 1:length(basis.coeffs[i]) iszero(basis.coeffs[i][j]) && error("Coefficient is zero") (ring.ch > 0) && @@ -455,7 +458,7 @@ function pairset_update!( # Criterion B_k(j, l). # T(j, l) is divisible by T(k) # max(deg T(j, k), deg T(l, k)) < deg T(j, l) - @inbounds for i in 1:pairset.load + @inbounds for i in 1:(pairset.load) (pairs[i].lcm == CRITICAL_PAIR_REDUNDANT) && continue j = pairs[i].poly1 l = pairs[i].poly2 @@ -470,9 +473,9 @@ function pairset_update!( @inbounds for i in 1:(idx - 1) (lcms[i] == CRITICAL_PAIR_REDUNDANT) && continue # if !basis.is_redundant[i] - pairs[pairset.load + cnt] = pairs[pairset.load + i] - degs[pairset.load + cnt] = degs[pairset.load + i] - cnt += 1 + pairs[pairset.load + cnt] = pairs[pairset.load + i] + degs[pairset.load + cnt] = degs[pairset.load + i] + cnt += 1 # end end @@ -491,9 +494,10 @@ function pairset_update!( divmask_i = update_ht.divmasks[lcms[i]] lcm_i = update_ht.monoms[lcms[i]] k = i + 1 - for k in (i+1):new_pairs + for k in (i + 1):new_pairs (lcms[k] == CRITICAL_PAIR_REDUNDANT) && continue - if update_ht.use_divmask && !divmask_is_probably_divisible(update_ht.divmasks[lcms[k]], divmask_i) + if update_ht.use_divmask && + !divmask_is_probably_divisible(update_ht.divmasks[lcms[k]], divmask_i) continue end ea = update_ht.monoms[lcms[k]] @@ -522,11 +526,11 @@ function pairset_update!( degs[cnt] = degs[pairset.load + i] cnt += 1 end - pairset.load = cnt - 1 + pairset.load = cnt - 1 # Mark redundant polynomials in the generating set. nonred = basis.nonredundant_indices - @inbounds for i in 1:basis.n_nonredundant + @inbounds for i in 1:(basis.n_nonredundant) basis.is_redundant[nonred[i]] && continue if hashtable_monom_is_divisible(basis.monoms[nonred[i]][1], lead_idx, ht) basis.is_redundant[nonred[i]] = true @@ -608,7 +612,7 @@ function basis_fill_data!( basis.n_filled = ngens end -function basis_move_redundant_elements!(basis::Basis) +function basis_discard_redundant_elements!(basis::Basis) j = 1 @inbounds for i in 1:(basis.n_nonredundant) if !basis.is_redundant[basis.nonredundant_indices[i]] diff --git a/src/f4/f4.jl b/src/f4/f4.jl index c209aa4c..bd17ee13 100644 --- a/src/f4/f4.jl +++ b/src/f4/f4.jl @@ -430,12 +430,11 @@ function f4!( f4_update!(pairset, basis, hashtable, update_ht) - matrix_reinitialize!(matrix, 0) hashtable_reinitialize!(symbol_ht) end - basis_move_redundant_elements!(basis) - + basis_discard_redundant_elements!(basis) + if params.reduced f4_autoreduce!(ring, basis, matrix, hashtable, symbol_ht, params) end diff --git a/src/f4/learn_apply.jl b/src/f4/learn_apply.jl index 38fb0e71..8b258c30 100644 --- a/src/f4/learn_apply.jl +++ b/src/f4/learn_apply.jl @@ -26,9 +26,7 @@ function f4_initialize_structs_with_trace( make_monic=make_monic, sort_input=sort_input ) - trace = trace_initialize(ring, basis_deepcopy(basis), basis, hashtable, permutation, params) - trace, basis, pairset, hashtable, permutation end @@ -129,7 +127,7 @@ function f4_reducegb_learn!( i = 1 @label Letsgo @inbounds while i <= basis.n_processed - @inbounds for j in 1:k + for j in 1:k if hashtable_monom_is_divisible( basis.monoms[basis.n_filled - i + 1][1], basis.monoms[basis.nonredundant_indices[j]][1], @@ -158,7 +156,6 @@ function f4_learn!( params::AlgorithmParameters ) where {M <: Monom, C <: Coeff} @invariant basis_well_formed(ring, basis, hashtable) - @invariant basis == trace.gb_basis @invariant params.reduced @@ -168,7 +165,7 @@ function f4_learn!( update_ht = hashtable_initialize_secondary(hashtable) symbol_ht = hashtable_initialize_secondary(hashtable) - pairset_size = f4_update!(pairset, basis, hashtable, update_ht) + f4_update!(pairset, basis, hashtable, update_ht) while !isempty(pairset) degree_i, npairs_i = f4_select_critical_pairs!(pairset, basis, matrix, hashtable, symbol_ht) @@ -178,13 +175,12 @@ function f4_learn!( f4_reduction_learn!(trace, basis, matrix, hashtable, symbol_ht, params) - pairset_size = f4_update!(pairset, basis, hashtable, update_ht) + f4_update!(pairset, basis, hashtable, update_ht) - matrix_reinitialize!(matrix, 0) hashtable_reinitialize!(symbol_ht) end - basis_move_redundant_elements!(basis) + basis_discard_redundant_elements!(basis) if params.reduced f4_reducegb_learn!(trace, ring, basis, matrix, hashtable, symbol_ht, params) @@ -469,16 +465,15 @@ function f4_apply!( basis::Basis{C}, params::AlgorithmParameters ) where {C <: Coeff} - ring = trace.ring @invariant basis_well_formed(ring, basis, trace.hashtable) @invariant params.reduced - basis_make_monic!(basis, params.arithmetic, params.changematrix) - iters_total = length(trace.matrix_infos) - 1 iters = 0 hashtable = trace.hashtable + basis_make_monic!(basis, params.arithmetic, params.changematrix) + symbol_ht = hashtable_initialize_secondary(hashtable) matrix = matrix_initialize(ring, C) @@ -501,10 +496,7 @@ function f4_apply!( cache_column_order, params ) - if !flag - # Unlucky cancellation of basis coefficients may have happened - return false - end + !flag && return false cache_column_order = new_cache_column_order && cache_column_order basis_update!(basis, hashtable) @@ -524,17 +516,14 @@ function f4_apply!( cache_column_order, params ) - if !flag - return false - end + !flag && return false end f4_standardize_basis_in_apply!(ring, trace, params.arithmetic) - basis = trace.gb_basis - - @invariant basis_well_formed(ring, basis, hashtable) trace.napply += 1 + @invariant basis_well_formed(ring, trace.gb_basis, hashtable) + true end diff --git a/src/groebner/groebner.jl b/src/groebner/groebner.jl index bb84222a..377dccf1 100644 --- a/src/groebner/groebner.jl +++ b/src/groebner/groebner.jl @@ -141,23 +141,58 @@ function get_next_batchsize(primes_used::Int, prev_batchsize::Int, batchsize_sca max(new_batchsize, prev_batchsize) end +function _groebner_guess_lucky_prime( + state::ModularState, + ring::PolyRing, + basis_zz::Basis, + pairset::Pairset, + hashtable::MonomialHashtable{M}, + params::AlgorithmParameters +) where {M <: Monom} + prime_1 = modular_random_prime(state, params.rng) + ring_ff_1, basis_ff_1 = modular_reduce_mod_p!(ring, basis_zz, prime_1, deepcopy=true) + params_zp = params_mod_p(params, prime_1) + f4!(ring_ff_1, basis_ff_1, pairset, hashtable, params_zp) + + prime_2 = modular_random_prime(state, params.rng) + ring_ff_2, basis_ff_2 = modular_reduce_mod_p!(ring, basis_zz, prime_2, deepcopy=true) + params_zp = params_mod_p(params, prime_2) + f4!(ring_ff_2, basis_ff_2, pairset, hashtable, params_zp) + + if basis_ff_1.monoms == basis_ff_2.monoms + return prime_1 + end + + prime_3 = modular_random_prime(state, params.rng) + ring_ff_3, basis_ff_3 = modular_reduce_mod_p!(ring, basis_zz, prime_3, deepcopy=true) + params_zp = params_mod_p(params, prime_3) + f4!(ring_ff_3, basis_ff_3, pairset, hashtable, params_zp) + + if basis_ff_1.monoms == basis_ff_3.monoms + return prime_1 + else + @assert basis_ff_2.monoms == basis_ff_3.monoms + return prime_2 + end + + prime_1 +end + function _groebner_learn_and_apply( ring::PolyRing, monoms::Vector{Vector{M}}, coeffs::Vector{Vector{C}}, params::AlgorithmParameters ) where {M <: Monom, C <: CoeffQQ} - # Initialize supporting structs - state = ModularState{BigInt, C, Int32}() basis, pairset, hashtable, permutation = f4_initialize_structs(ring, monoms, coeffs, params, make_monic=false) basis_zz = clear_denominators!(basis, deepcopy=false) - # Handler for lucky primes - lucky = LuckyPrimes(basis_zz.coeffs) - prime = primes_next_lucky_prime!(lucky) + state = ModularState{BigInt, C, Int32}(basis_zz.coeffs) + + prime = _groebner_guess_lucky_prime(state, ring, basis_zz, pairset, hashtable, params) ring_ff, basis_ff = modular_reduce_mod_p!(ring, basis_zz, prime, deepcopy=true) @@ -174,33 +209,26 @@ function _groebner_learn_and_apply( f4_learn!(trace, ring_ff, trace.gb_basis, pairset, hashtable, params_zp) # TODO: no need to deepcopy! + push!(state.used_primes, prime) push!(state.gb_coeffs_ff_all, deepcopy(trace.gb_basis.coeffs)) # Reconstruct coefficients and write results to the accumulator. modular_prepare!(state) crt_vec_full!( state.gb_coeffs_zz, - lucky.modulo, + state.modulo, state.gb_coeffs_ff_all, - lucky.used_primes, + state.used_primes, state.crt_mask ) success_reconstruct = - ratrec_vec_full!(state.gb_coeffs_qq, state.gb_coeffs_zz, lucky.modulo, state.ratrec_mask) + ratrec_vec_full!(state.gb_coeffs_qq, state.gb_coeffs_zz, state.modulo, state.ratrec_mask) correct_basis = false if success_reconstruct - correct_basis = modular_lift_check!( - state, - lucky, - ring_ff, - basis, - basis_zz, - trace.gb_basis, - hashtable, - params - ) + correct_basis = + modular_lift_check!(state, ring_ff, basis, basis_zz, trace.gb_basis, hashtable, params) # At this point, the constructed basis is deemed correct, we return it. if correct_basis gb_monoms, _ = basis_export_data(trace.gb_basis, hashtable) @@ -224,7 +252,7 @@ function _groebner_learn_and_apply( while !correct_basis if iszero(batchsize % 4) && params.batched for j in 1:4:batchsize - prime_4x = ntuple(_ -> Int32(primes_next_lucky_prime!(lucky)), 4) + prime_4x = ntuple(i -> Int32(modular_next_prime!(state)), 4) # Perform reduction modulo primes and store result in basis_ff_4x ring_ff_4x, basis_ff_4x = modular_reduce_mod_p_in_batch!(ring, basis_zz, prime_4x) @@ -241,6 +269,7 @@ function _groebner_learn_and_apply( ir_unpack_composite_coefficients(trace_4x.gb_basis.coeffs) # TODO: This causes unnecessary conversions of arrays. + append!(state.used_primes, prime_4x) push!(state.gb_coeffs_ff_all, gb_coeffs_1) push!(state.gb_coeffs_ff_all, gb_coeffs_2) push!(state.gb_coeffs_ff_all, gb_coeffs_3) @@ -249,7 +278,7 @@ function _groebner_learn_and_apply( end else for j in 1:batchsize - prime = primes_next_lucky_prime!(lucky) + prime = modular_next_prime!(state) ring_ff, basis_ff = modular_reduce_mod_p!(ring, basis_zz, prime, deepcopy=true) params_zp = params_mod_p(params, prime) @@ -259,6 +288,7 @@ function _groebner_learn_and_apply( f4_apply!(trace, ring_ff, trace.buf_basis, params_zp) + push!(state.used_primes, prime) push!(state.gb_coeffs_ff_all, deepcopy(trace.gb_basis.coeffs)) if !modular_majority_vote!(state, trace.gb_basis, params) @@ -270,16 +300,16 @@ function _groebner_learn_and_apply( crt_vec_partial!( state.gb_coeffs_zz, - lucky.modulo, + state.modulo, state.gb_coeffs_ff_all, - lucky.used_primes, + state.used_primes, witness_set, state.crt_mask ) success_reconstruct = ratrec_vec_partial!( state.gb_coeffs_qq, state.gb_coeffs_zz, - lucky.modulo, + state.modulo, witness_set, state.ratrec_mask ) @@ -292,7 +322,7 @@ function _groebner_learn_and_apply( if params.heuristic_check success_check = - modular_lift_heuristic_check_partial(state.gb_coeffs_qq, lucky.modulo, witness_set) + modular_lift_heuristic_check_partial(state.gb_coeffs_qq, state.modulo, witness_set) if !success_check iters += 1 batchsize = get_next_batchsize(primes_used, batchsize, batchsize_scaling) @@ -303,16 +333,16 @@ function _groebner_learn_and_apply( # Perform full reconstruction crt_vec_full!( state.gb_coeffs_zz, - lucky.modulo, + state.modulo, state.gb_coeffs_ff_all, - lucky.used_primes, + state.used_primes, state.crt_mask ) success_reconstruct = ratrec_vec_full!( state.gb_coeffs_qq, state.gb_coeffs_zz, - lucky.modulo, + state.modulo, state.ratrec_mask ) @@ -322,16 +352,8 @@ function _groebner_learn_and_apply( continue end - correct_basis = modular_lift_check!( - state, - lucky, - ring_ff, - basis, - basis_zz, - trace.gb_basis, - hashtable, - params - ) + correct_basis = + modular_lift_check!(state, ring_ff, basis, basis_zz, trace.gb_basis, hashtable, params) iters += 1 batchsize = get_next_batchsize(primes_used, batchsize, batchsize_scaling) @@ -358,15 +380,14 @@ function _groebner_learn_and_apply_threaded( end # Initialize supporting structs - state = ModularState{BigInt, C, Int32}() basis, pairset, hashtable, permutation = f4_initialize_structs(ring, monoms, coeffs, params, make_monic=false) basis_zz = clear_denominators!(basis, deepcopy=false) - # Handler for lucky primes - lucky = LuckyPrimes(basis_zz.coeffs) - prime = primes_next_lucky_prime!(lucky) + state = ModularState{BigInt, C, Int32}(basis_zz.coeffs) + + prime = _groebner_guess_lucky_prime(state, ring, basis_zz, pairset, hashtable, params) ring_ff, basis_ff = modular_reduce_mod_p!(ring, basis_zz, prime, deepcopy=true) @@ -383,33 +404,26 @@ function _groebner_learn_and_apply_threaded( f4_learn!(trace, ring_ff, trace.gb_basis, pairset, hashtable, params_zp) # TODO: no need to deepcopy! + push!(state.used_primes, prime) push!(state.gb_coeffs_ff_all, deepcopy(trace.gb_basis.coeffs)) # Reconstruct coefficients and write results to the accumulator. modular_prepare!(state) crt_vec_full!( state.gb_coeffs_zz, - lucky.modulo, + state.modulo, state.gb_coeffs_ff_all, - lucky.used_primes, + state.used_primes, state.crt_mask ) success_reconstruct = - ratrec_vec_full!(state.gb_coeffs_qq, state.gb_coeffs_zz, lucky.modulo, state.ratrec_mask) + ratrec_vec_full!(state.gb_coeffs_qq, state.gb_coeffs_zz, state.modulo, state.ratrec_mask) correct_basis = false if success_reconstruct - correct_basis = modular_lift_check!( - state, - lucky, - ring_ff, - basis, - basis_zz, - trace.gb_basis, - hashtable, - params - ) + correct_basis = + modular_lift_check!(state, ring_ff, basis, basis_zz, trace.gb_basis, hashtable, params) if correct_basis gb_monoms, _ = basis_export_data(trace.gb_basis, hashtable) gb_coeffs_qq = state.gb_coeffs_qq @@ -441,7 +455,7 @@ function _groebner_learn_and_apply_threaded( while !correct_basis @invariant iszero(batchsize % 4) - threadbuf_primes = map(_ -> Int32(primes_next_lucky_prime!(lucky)), 1:batchsize) + threadbuf_primes = ntuple(_ -> Int32(modular_next_prime!(state)), batchsize) for i in 1:nthreads() empty!(threadbuf_gb_coeffs[i]) end @@ -483,22 +497,23 @@ function _groebner_learn_and_apply_threaded( threadbuf_gb_coeffs_union = reduce(vcat, threadbuf_gb_coeffs) sort!(threadbuf_gb_coeffs_union, by=first, rev=true) - for (_, coeffs_ff_) in threadbuf_gb_coeffs_union + for (prime_, coeffs_ff_) in threadbuf_gb_coeffs_union + push!(state.used_primes, prime_) push!(state.gb_coeffs_ff_all, coeffs_ff_) end crt_vec_partial!( state.gb_coeffs_zz, - lucky.modulo, + state.modulo, state.gb_coeffs_ff_all, - lucky.used_primes, + state.used_primes, witness_set, state.crt_mask ) success_reconstruct = ratrec_vec_partial!( state.gb_coeffs_qq, state.gb_coeffs_zz, - lucky.modulo, + state.modulo, witness_set, state.ratrec_mask ) @@ -511,7 +526,7 @@ function _groebner_learn_and_apply_threaded( if params.heuristic_check success_check = - modular_lift_heuristic_check_partial(state.gb_coeffs_qq, lucky.modulo, witness_set) + modular_lift_heuristic_check_partial(state.gb_coeffs_qq, state.modulo, witness_set) if !success_check iters += 1 batchsize = get_next_batchsize(primes_used, batchsize, batchsize_scaling) @@ -522,15 +537,15 @@ function _groebner_learn_and_apply_threaded( # Perform full reconstruction crt_vec_full!( state.gb_coeffs_zz, - lucky.modulo, + state.modulo, state.gb_coeffs_ff_all, - lucky.used_primes, + state.used_primes, state.crt_mask ) success_reconstruct = ratrec_vec_full!( state.gb_coeffs_qq, state.gb_coeffs_zz, - lucky.modulo, + state.modulo, state.ratrec_mask ) @@ -541,16 +556,8 @@ function _groebner_learn_and_apply_threaded( continue end - correct_basis = modular_lift_check!( - state, - lucky, - ring_ff, - basis, - basis_zz, - trace.gb_basis, - hashtable, - params - ) + correct_basis = + modular_lift_check!(state, ring_ff, basis, basis_zz, trace.gb_basis, hashtable, params) iters += 1 batchsize = get_next_batchsize(primes_used, batchsize, batchsize_scaling) @@ -571,17 +578,15 @@ function _groebner_classic_modular( coeffs::Vector{Vector{C}}, params::AlgorithmParameters ) where {M <: Monom, C <: CoeffQQ} - # Initialize supporting structs - state = ModularState{BigInt, C, CoeffModular}() basis, pairset, hashtable = f4_initialize_structs(ring, monoms, coeffs, params, make_monic=false) basis_zz = clear_denominators!(basis, deepcopy=false) - # Handler for lucky primes - lucky = LuckyPrimes(basis_zz.coeffs) - prime = primes_next_lucky_prime!(lucky) + state = ModularState{BigInt, C, CoeffModular}(basis_zz.coeffs) + + prime = _groebner_guess_lucky_prime(state, ring, basis_zz, pairset, hashtable, params) ring_ff, basis_ff = modular_reduce_mod_p!(ring, basis_zz, prime, deepcopy=true) @@ -590,25 +595,26 @@ function _groebner_classic_modular( # NOTE: basis_ff may not own its coefficients, one should not mutate it # directly further in the code + push!(state.used_primes, prime) push!(state.gb_coeffs_ff_all, basis_ff.coeffs) # Reconstruct coefficients and write results to the accumulator. modular_prepare!(state) crt_vec_full!( state.gb_coeffs_zz, - lucky.modulo, + state.modulo, state.gb_coeffs_ff_all, - lucky.used_primes, + state.used_primes, state.crt_mask ) success_reconstruct = - ratrec_vec_full!(state.gb_coeffs_qq, state.gb_coeffs_zz, lucky.modulo, state.ratrec_mask) + ratrec_vec_full!(state.gb_coeffs_qq, state.gb_coeffs_zz, state.modulo, state.ratrec_mask) correct_basis = false if success_reconstruct correct_basis = - modular_lift_check!(state, lucky, ring_ff, basis, basis_zz, basis_ff, hashtable, params) + modular_lift_check!(state, ring_ff, basis, basis_zz, basis_ff, hashtable, params) if correct_basis gb_monoms, _ = basis_export_data(basis_ff, hashtable) gb_coeffs_qq = state.gb_coeffs_qq @@ -628,13 +634,14 @@ function _groebner_classic_modular( iters = 0 while !correct_basis for j in 1:batchsize - prime = primes_next_lucky_prime!(lucky) + prime = modular_next_prime!(state) ring_ff, basis_ff = modular_reduce_mod_p!(ring, basis_zz, prime, deepcopy=true) params_zp = params_mod_p(params, prime) f4!(ring_ff, basis_ff, pairset, hashtable, params_zp) + push!(state.used_primes, prime) push!(state.gb_coeffs_ff_all, basis_ff.coeffs) if !modular_majority_vote!(state, basis_ff, params) @@ -645,16 +652,16 @@ function _groebner_classic_modular( crt_vec_partial!( state.gb_coeffs_zz, - lucky.modulo, + state.modulo, state.gb_coeffs_ff_all, - lucky.used_primes, + state.used_primes, witness_set, state.crt_mask ) success_reconstruct = ratrec_vec_partial!( state.gb_coeffs_qq, state.gb_coeffs_zz, - lucky.modulo, + state.modulo, witness_set, state.ratrec_mask ) @@ -667,7 +674,7 @@ function _groebner_classic_modular( if params.heuristic_check success_check = - modular_lift_heuristic_check_partial(state.gb_coeffs_qq, lucky.modulo, witness_set) + modular_lift_heuristic_check_partial(state.gb_coeffs_qq, state.modulo, witness_set) if !success_check iters += 1 batchsize = get_next_batchsize(primes_used, batchsize, batchsize_scaling) @@ -678,15 +685,15 @@ function _groebner_classic_modular( # Perform full reconstruction crt_vec_full!( state.gb_coeffs_zz, - lucky.modulo, + state.modulo, state.gb_coeffs_ff_all, - lucky.used_primes, + state.used_primes, state.crt_mask ) success_reconstruct = ratrec_vec_full!( state.gb_coeffs_qq, state.gb_coeffs_zz, - lucky.modulo, + state.modulo, state.ratrec_mask ) @@ -697,7 +704,7 @@ function _groebner_classic_modular( end correct_basis = - modular_lift_check!(state, lucky, ring_ff, basis, basis_zz, basis_ff, hashtable, params) + modular_lift_check!(state, ring_ff, basis, basis_zz, basis_ff, hashtable, params) iters += 1 batchsize = get_next_batchsize(primes_used, batchsize, batchsize_scaling) diff --git a/src/groebner/groebner_with_change_matrix.jl b/src/groebner/groebner_with_change_matrix.jl index 2396e054..0939eec7 100644 --- a/src/groebner/groebner_with_change_matrix.jl +++ b/src/groebner/groebner_with_change_matrix.jl @@ -139,17 +139,15 @@ function _groebner_with_change_classic_modular( coeffs::Vector{Vector{C}}, params::AlgorithmParameters ) where {M <: Monom, C <: CoeffQQ} - # Initialize supporting structs - state = ModularState{BigInt, C, CoeffModular}() basis, pairset, hashtable = f4_initialize_structs(ring, monoms, coeffs, params, make_monic=false) basis_zz = clear_denominators!(basis, deepcopy=false) - # Handler for lucky primes - lucky = LuckyPrimes(basis_zz.coeffs) - prime = primes_next_lucky_prime!(lucky) + state = ModularState{BigInt, C, CoeffModular}(basis_zz.coeffs) + + prime = modular_next_prime!(state) # Perform reduction modulo prime and store result in basis_ff ring_ff, basis_ff = modular_reduce_mod_p!(ring, basis_zz, prime, deepcopy=true) @@ -159,6 +157,7 @@ function _groebner_with_change_classic_modular( # NOTE: basis_ff may not own its coefficients, one should not mutate it # directly further in the code + push!(state.used_primes, prime) push!(state.gb_coeffs_ff_all, basis_ff.coeffs) changematrix_monoms, changematrix_coeffs = @@ -169,22 +168,22 @@ function _groebner_with_change_classic_modular( modular_prepare!(state) crt_vec_full!( state.gb_coeffs_zz, - lucky.modulo, + state.modulo, state.gb_coeffs_ff_all, - lucky.used_primes, + state.used_primes, state.crt_mask ) - modular_crt_full_changematrix!(state, lucky) + modular_crt_full_changematrix!(state) success_reconstruct = - ratrec_vec_full!(state.gb_coeffs_qq, state.gb_coeffs_zz, lucky.modulo, state.ratrec_mask) + ratrec_vec_full!(state.gb_coeffs_qq, state.gb_coeffs_zz, state.modulo, state.ratrec_mask) - changematrix_success_reconstruct = modular_ratrec_full_changematrix!(state, lucky) + changematrix_success_reconstruct = modular_ratrec_full_changematrix!(state) correct_basis = false if success_reconstruct && changematrix_success_reconstruct correct_basis = - modular_lift_check!(state, lucky, ring_ff, basis, basis_zz, basis_ff, hashtable, params) + modular_lift_check!(state, ring_ff, basis, basis_zz, basis_ff, hashtable, params) # At this point, if the constructed basis is correct, we return it. if correct_basis gb_monoms, _ = basis_export_data(basis_ff, hashtable) @@ -204,9 +203,9 @@ function _groebner_with_change_classic_modular( crt_vec_partial!( state.gb_coeffs_zz, - lucky.modulo, + state.modulo, state.gb_coeffs_ff_all, - lucky.used_primes, + state.used_primes, witness_set, state.crt_mask ) @@ -214,13 +213,14 @@ function _groebner_with_change_classic_modular( iters = 0 while !correct_basis for j in 1:batchsize - prime = primes_next_lucky_prime!(lucky) + prime = modular_next_prime!(state) ring_ff, basis_ff = modular_reduce_mod_p!(ring, basis_zz, prime, deepcopy=true) params_zp = params_mod_p(params, prime) f4!(ring_ff, basis_ff, pairset, hashtable, params_zp) + push!(state.used_primes, prime) push!(state.gb_coeffs_ff_all, basis_ff.coeffs) _, changematrix_coeffs = basis_changematrix_export(basis_ff, hashtable, length(monoms)) push!(state.changematrix_coeffs_ff_all, changematrix_coeffs) @@ -234,9 +234,9 @@ function _groebner_with_change_classic_modular( crt_vec_partial!( state.gb_coeffs_zz, - lucky.modulo, + state.modulo, state.gb_coeffs_ff_all, - lucky.used_primes, + state.used_primes, witness_set, state.crt_mask ) @@ -244,7 +244,7 @@ function _groebner_with_change_classic_modular( success_reconstruct = ratrec_vec_partial!( state.gb_coeffs_qq, state.gb_coeffs_zz, - lucky.modulo, + state.modulo, witness_set, state.ratrec_mask ) @@ -257,7 +257,7 @@ function _groebner_with_change_classic_modular( if params.heuristic_check success_check = - modular_lift_heuristic_check_partial(state.gb_coeffs_qq, lucky.modulo, witness_set) + modular_lift_heuristic_check_partial(state.gb_coeffs_qq, state.modulo, witness_set) if !success_check iters += 1 batchsize = get_next_batchsize(primes_used, batchsize, batchsize_scaling) @@ -268,16 +268,16 @@ function _groebner_with_change_classic_modular( # Perform full reconstruction crt_vec_full!( state.gb_coeffs_zz, - lucky.modulo, + state.modulo, state.gb_coeffs_ff_all, - lucky.used_primes, + state.used_primes, state.crt_mask ) success_reconstruct = ratrec_vec_full!( state.gb_coeffs_qq, state.gb_coeffs_zz, - lucky.modulo, + state.modulo, state.ratrec_mask ) @@ -287,8 +287,8 @@ function _groebner_with_change_classic_modular( continue end - modular_crt_full_changematrix!(state, lucky) - success_reconstruct_changematrix = modular_ratrec_full_changematrix!(state, lucky) + modular_crt_full_changematrix!(state) + success_reconstruct_changematrix = modular_ratrec_full_changematrix!(state) if !success_reconstruct_changematrix iters += 1 batchsize = get_next_batchsize(primes_used, batchsize, batchsize_scaling) @@ -296,7 +296,7 @@ function _groebner_with_change_classic_modular( end correct_basis = - modular_lift_check!(state, lucky, ring_ff, basis, basis_zz, basis_ff, hashtable, params) + modular_lift_check!(state, ring_ff, basis, basis_zz, basis_ff, hashtable, params) iters += 1 batchsize = get_next_batchsize(primes_used, batchsize, batchsize_scaling) diff --git a/src/groebner/isgroebner.jl b/src/groebner/isgroebner.jl index 53c729e1..d23e6309 100644 --- a/src/groebner/isgroebner.jl +++ b/src/groebner/isgroebner.jl @@ -72,8 +72,8 @@ function _isgroebner2( end # Otherwise, check modulo a prime basis_zz = clear_denominators!(basis, deepcopy=false) - luckyprimes = LuckyPrimes(basis_zz.coeffs) - prime = primes_next_aux_prime!(luckyprimes) + state = ModularState{BigInt, Rational{BigInt}, UInt32}(basis_zz.coeffs) + prime = modular_random_prime(state, params.rng) ring_ff, basis_ff = modular_reduce_mod_p!(ring, basis_zz, prime, deepcopy=true) arithmetic = select_arithmetic(CoeffModular, prime, :auto, false) flag = f4_isgroebner!(ring_ff, basis_ff, pairset, hashtable, arithmetic) diff --git a/src/groebner/modular.jl b/src/groebner/modular.jl index 3683a607..b646e7ff 100644 --- a/src/groebner/modular.jl +++ b/src/groebner/modular.jl @@ -1,10 +1,23 @@ # This file is a part of Groebner.jl. License is GNU GPL v2. +# The sequence of lucky prime candidates is decreasing and deterministic. +# Make sure this number is at most 31 bits to be able to use signed ints. +const FIRST_LUCKY_PRIME_CANDIDATE = 2^31 - 1 + +const RANDOM_PRIME_LB = 2^29 +const RANDOM_PRIME_UB = 2^30 + mutable struct ModularState{T1 <: CoeffZZ, T2 <: CoeffQQ, T3} gb_coeffs_zz::Vector{Vector{T1}} gb_coeffs_qq::Vector{Vector{T2}} gb_coeffs_ff_all::Vector{Vector{Vector{T3}}} + gen_coeffs_zz::Vector{Vector{T1}} + + modulo::BigInt + current_prime::UInt64 + used_primes::Vector{UInt64} + crt_mask::Vector{BitVector} ratrec_mask::Vector{BitVector} @@ -12,11 +25,17 @@ mutable struct ModularState{T1 <: CoeffZZ, T2 <: CoeffQQ, T3} changematrix_coeffs_zz::Vector{Vector{Vector{T1}}} changematrix_coeffs_qq::Vector{Vector{Vector{T2}}} - function ModularState{T1, T2, T3}() where {T1 <: CoeffZZ, T2 <: CoeffQQ, T3 <: CoeffZp} + function ModularState{T1, T2, T3}( + coeffs_zz::Vector{Vector{T1}} + ) where {T1 <: CoeffZZ, T2 <: CoeffQQ, T3 <: CoeffZp} new( Vector{Vector{T1}}(), Vector{Vector{T2}}(), Vector{Vector{Vector{T3}}}(), + coeffs_zz, + BigInt(), + UInt64(FIRST_LUCKY_PRIME_CANDIDATE), + Vector{UInt64}(), Vector{BitVector}(), Vector{BitVector}(), Vector{Vector{Vector{Vector{T3}}}}(), @@ -26,6 +45,47 @@ mutable struct ModularState{T1 <: CoeffZZ, T2 <: CoeffQQ, T3} end end +### +# Prime number selection + +function modular_accept_prime(state::ModularState, prime::UInt64) + p = BigInt(prime) + buf = BigInt() + @inbounds for coeffs in state.gen_coeffs_zz + c = coeffs[1] + if Base.GMP.MPZ.cmp_ui(Base.GMP.MPZ.tdiv_r!(buf, c, p), 0) == 0 + return false + end + c = coeffs[end] + if Base.GMP.MPZ.cmp_ui(Base.GMP.MPZ.tdiv_r!(buf, c, p), 0) == 0 + return false + end + end + true +end + +function modular_next_prime!(state::ModularState) + prime = Primes.prevprime(state.current_prime - 1) + while !modular_accept_prime(state, prime) + prime = Primes.prevprime(prime - 1) + end + state.current_prime = prime + prime +end + +function modular_random_prime(state::ModularState, rng::AbstractRNG) + lb, ub = RANDOM_PRIME_LB, RANDOM_PRIME_UB + while true + candidate = rand(rng, lb:ub) % UInt64 + !modular_accept_prime(state, candidate) && continue + Primes.isprime(candidate) && return UInt64(candidate) + end + UInt64(0) +end + +### +# Other stuff + function modular_witness_set( gb_coeffs::Vector{Vector{C}}, params::AlgorithmParameters @@ -183,12 +243,12 @@ end ### # Lifting change matrix -function modular_ratrec_full_changematrix!(state::ModularState, lucky::LuckyPrimes) +function modular_ratrec_full_changematrix!(state::ModularState) @inbounds for i in 1:length(state.changematrix_coeffs_zz) flag = ratrec_vec_full!( state.changematrix_coeffs_qq[i], state.changematrix_coeffs_zz[i], - lucky.modulo, + state.modulo, [ falses(length(state.changematrix_coeffs_qq[i][j])) for j in 1:length(state.changematrix_coeffs_qq[i]) @@ -199,7 +259,7 @@ function modular_ratrec_full_changematrix!(state::ModularState, lucky::LuckyPrim true end -function modular_crt_full_changematrix!(state::ModularState, lucky::LuckyPrimes) +function modular_crt_full_changematrix!(state::ModularState) if isempty(state.changematrix_coeffs_zz) coeffs_ff = state.changematrix_coeffs_ff_all[1] resize!(state.changematrix_coeffs_zz, length(coeffs_ff)) @@ -221,16 +281,14 @@ function modular_crt_full_changematrix!(state::ModularState, lucky::LuckyPrimes) crt_vec_full!( state.changematrix_coeffs_zz[i], modulo, - [state.changematrix_coeffs_ff_all[ell][i] for ell in 1:length(lucky.used_primes)], - lucky.used_primes, + [state.changematrix_coeffs_ff_all[ell][i] for ell in 1:length(state.used_primes)], + state.used_primes, [ falses(length(state.changematrix_coeffs_ff_all[1][i][j])) for j in 1:length(state.changematrix_coeffs_ff_all[1][i]) ] ) end - - nothing end ### @@ -238,7 +296,6 @@ end function modular_lift_check!( state::ModularState, - lucky::LuckyPrimes, ring::PolyRing, basis_qq::Basis, basis_zz::Basis, @@ -248,22 +305,14 @@ function modular_lift_check!( ) # First we check the size of the coefficients with a heuristic if params.heuristic_check - if !modular_lift_heuristic_check(state.gb_coeffs_qq, lucky.modulo) + if !modular_lift_heuristic_check(state.gb_coeffs_qq, state.modulo) return false end end # Then check that a basis is also a basis modulo a prime if params.randomized_check - if !modular_lift_randomized_check!( - state, - ring, - basis_zz, - basis_ff, - lucky, - hashtable, - params - ) + if !modular_lift_randomized_check!(state, ring, basis_zz, basis_ff, hashtable, params) return false end end @@ -318,12 +367,11 @@ function modular_lift_randomized_check!( ring::PolyRing, input_zz::Basis, gb_ff::Basis, - lucky::LuckyPrimes, hashtable::MonomialHashtable, params::AlgorithmParameters ) # !!! note that this function may modify the given hashtable! - prime = primes_next_aux_prime!(lucky) + prime = modular_random_prime(state, params.rng) ring_ff, input_ff = modular_reduce_mod_p!(ring, input_zz, prime, deepcopy=true) # TODO: do we really need to re-scale things to be fraction-free? gb_coeffs_zz = _clear_denominators!(state.gb_coeffs_qq) @@ -357,14 +405,14 @@ function modular_lift_certify_check!( gb_qq = basis_deep_copy_with_new_coeffs(gb_ff, state.gb_coeffs_qq) ring_qq = PolyRing(ring.nvars, ring.ord, 0) input_qq = basis_deepcopy(input_qq) - # Check that some polynomial is not reduced to zero modulo a prime + # Check that some polynomial is not reduced to zero f4_normalform!(ring_qq, gb_qq, input_qq, hashtable, params.arithmetic) for i in 1:(input_qq.n_processed) if !isempty(input_qq.coeffs[i]) return false end end - # Check that the basis is a groebner basis modulo a prime + # Check that the basis is a groebner basis pairset = pairset_initialize(UInt64) if !f4_isgroebner!(ring_qq, gb_qq, pairset, hashtable, params.arithmetic) return false diff --git a/src/groebner/primes.jl b/src/groebner/primes.jl deleted file mode 100644 index 9f1f5c94..00000000 --- a/src/groebner/primes.jl +++ /dev/null @@ -1,88 +0,0 @@ -# This file is a part of Groebner.jl. License is GNU GPL v2. - -### -# Lucky primes - -# The sequence of lucky prime candidates is decreasing and deterministic. -# Groebner basis is computed modulo a lucky prime, and the correctness of the -# computed basis is checked modulo some another prime(s). -# -# We use 31 bit primes for modular computation, since it allows the use of -# signed integers for implementing arithmetic in the ground field. - -const FIRST_LUCKY_PRIME_CANDIDATE = 2^31 - 1 -const FIRST_AUX_PRIME_CANDIDATE = 2^30 + 3 - -# Keeps track of used prime numbers and helps selecting new ones -mutable struct LuckyPrimes - # Integer coefficients of input polynomials - coeffs::Vector{Vector{BigInt}} - buf::BigInt - lucky_prime::UInt64 - aux_prime::UInt64 - used_primes::Vector{UInt64} - # product of used lucky primes - modulo::BigInt - - function LuckyPrimes(coeffs::Vector{Vector{BigInt}}) - new( - coeffs, - BigInt(), - FIRST_LUCKY_PRIME_CANDIDATE, - FIRST_AUX_PRIME_CANDIDATE, - UInt64[], - BigInt(1) - ) - end -end - -# Check if the prime is lucky, that is, if the prime does not divide any of the -# leading coefficients of the input. -function primes_is_lucky_prime(lucky::LuckyPrimes, prime::UInt64) - buf = lucky.buf - p = BigInt(prime) - @inbounds for coeffs in lucky.coeffs - @invariant !isempty(coeffs) - # TODO(Sasha): I am too greedy to check all coefficients. - # for c in coeffs - # if Base.GMP.MPZ.cmp_ui(Base.GMP.MPZ.tdiv_r!(buf, c, p), 0) == 0 - # return false - # end - # end - c = coeffs[1] - if Base.GMP.MPZ.cmp_ui(Base.GMP.MPZ.tdiv_r!(buf, c, p), 0) == 0 - return false - end - c = coeffs[end] - if Base.GMP.MPZ.cmp_ui(Base.GMP.MPZ.tdiv_r!(buf, c, p), 0) == 0 - return false - end - end - true -end - -function primes_next_lucky_prime!(lucky::LuckyPrimes) - prime = lucky.lucky_prime - while !primes_is_lucky_prime(lucky, prime) - prime = Primes.prevprime(prime - 1) - end - lucky.lucky_prime = Primes.prevprime(prime - 1) - @invariant !(prime in lucky.used_primes) - push!(lucky.used_primes, prime) - prime -end - -function primes_next_aux_prime!(lucky::LuckyPrimes) - prime = lucky.aux_prime - while !primes_is_lucky_prime(lucky, prime) - prime = nextprime(prime + 1) - prime >= FIRST_LUCKY_PRIME_CANDIDATE && throw( - DomainError( - "Something went wrong in Groebner.jl. Please consider submitting a GitHub issue." - ) - ) - end - lucky.aux_prime = nextprime(prime + 1) - @invariant !(prime in lucky.used_primes) - prime -end diff --git a/src/reconstruction/crt.jl b/src/reconstruction/crt.jl index c99c1a00..05f6ffb7 100644 --- a/src/reconstruction/crt.jl +++ b/src/reconstruction/crt.jl @@ -70,10 +70,8 @@ end ### # Element-wise CRT -# table of big integers must be initialized -# -# A mod M, a_1 mod m_1, ..., a_n mod m_n => B mod M m_1 ... m_n. # Reconstructs only the witness set. +# Table of big integers must be initialized. function crt_vec_partial!( table_zz::Vector{Vector{BigInt}}, modulo::BigInt, @@ -82,23 +80,8 @@ function crt_vec_partial!( witness_set::Vector{Tuple{Int, Int}}, mask::Vector{BitVector} ) where {T <: Integer, U <: Integer} - @invariant isbitstype(T) @invariant length(tables_ff) == length(moduli) - @invariant all(<(typemax(UInt64)), moduli) - - # Base case - if length(moduli) == 1 - table_ff = tables_ff[1] - Base.GMP.MPZ.set_ui!(modulo, UInt64(moduli[1])) - @invariant length(table_zz) == length(table_ff) - @inbounds for k in 1:length(witness_set) - i, j = witness_set[k] - @invariant 1 <= i <= length(table_ff) && 1 <= j <= length(table_ff[i]) - rem_ij = UInt64(table_ff[i][j]) - Base.GMP.MPZ.set_ui!(table_zz[i][j], rem_ij) - end - return nothing - end + @invariant length(table_zz) == length(mask) # Precompute CRT multipliers buf, n1, n2 = BigInt(), BigInt(), BigInt() @@ -106,12 +89,12 @@ function crt_vec_partial!( for i in 1:length(mults) mults[i] = BigInt(0) end + crt_precompute!(modulo, n1, n2, mults, map(UInt64, moduli)) rems = Vector{UInt64}(undef, length(moduli)) @inbounds for k in 1:length(witness_set) i, j = witness_set[k] - @invariant 1 <= i <= length(tables_ff[1]) && 1 <= j <= length(tables_ff[1][i]) for t in 1:length(moduli) rems[t] = UInt64(tables_ff[t][i][j]) @@ -127,7 +110,7 @@ function crt_vec_partial!( nothing end -# A mod M, a_1 mod m_1, ..., a_n mod m_n => B mod M m_1 ... m_n. +# Table of big integers must be initialized. function crt_vec_full!( table_zz::Vector{Vector{BigInt}}, modulo::BigInt, @@ -135,41 +118,23 @@ function crt_vec_full!( moduli::Vector{U}, mask::Vector{BitVector} ) where {T <: Integer, U <: Integer} - @invariant isbitstype(T) @invariant length(tables_ff) == length(moduli) - @invariant all(<(typemax(UInt64)), moduli) - - # Base case - if length(moduli) == 1 - table_ff = tables_ff[1] - Base.GMP.MPZ.set_ui!(modulo, UInt64(moduli[1])) - @inbounds for i in 1:length(table_zz) - @invariant length(table_zz[i]) == length(table_ff[i]) - for j in 1:length(table_zz[i]) - rem_ij = UInt64(table_ff[i][j]) - @invariant 0 <= rem_ij < moduli[1] - Base.GMP.MPZ.set_ui!(table_zz[i][j], rem_ij) - end - end - return nothing - end + @invariant length(table_zz) == length(mask) buf, n1, n2 = BigInt(), BigInt(), BigInt() mults = Vector{BigInt}(undef, length(moduli)) for i in 1:length(mults) mults[i] = BigInt(0) end + crt_precompute!(modulo, n1, n2, mults, map(UInt64, moduli)) rems = Vector{UInt64}(undef, length(moduli)) @inbounds for i in 1:length(table_zz) for j in 1:length(table_zz[i]) - if mask[i][j] - continue - end + mask[i][j] && continue for k in 1:length(moduli) - @invariant length(table_zz[i]) == length(tables_ff[k][i]) @invariant 0 <= tables_ff[k][i][j] < moduli[k] rems[k] = UInt64(tables_ff[k][i][j]) end @@ -179,6 +144,4 @@ function crt_vec_full!( Base.GMP.MPZ.set!(table_zz[i][j], buf) end end - - nothing end diff --git a/src/reconstruction/ratrec.jl b/src/reconstruction/ratrec.jl index c6d89f21..0b34b2ba 100644 --- a/src/reconstruction/ratrec.jl +++ b/src/reconstruction/ratrec.jl @@ -16,7 +16,8 @@ end ### # Element-wise rational reconstruction -# table of rationals need not be initialized +# Table of rationals need not be initialized. +# Reconstructs only the witness set. function ratrec_vec_partial!( table_qq::Vector{Vector{Rational{BigInt}}}, table_zz::Vector{Vector{BigInt}}, @@ -25,6 +26,7 @@ function ratrec_vec_partial!( mask::Vector{BitVector} ) @invariant length(table_qq) == length(table_zz) + @invariant modulo > 1 bound = ratrec_reconstruction_bound(modulo) nemo_bound = Nemo.ZZRingElem(bound) @@ -47,6 +49,7 @@ function ratrec_vec_partial!( true end +# Table of rationals need not be initialized. function ratrec_vec_full!( table_qq::Vector{Vector{Rational{BigInt}}}, table_zz::Vector{Vector{BigInt}}, @@ -63,9 +66,7 @@ function ratrec_vec_full!( @inbounds for i in 1:length(table_zz) @invariant length(table_zz[i]) == length(table_qq[i]) for j in 1:length(table_zz[i]) - if mask[i][j] - continue - end + mask[i][j] && continue rem_nemo = Nemo.ZZRingElem(table_zz[i][j]) @invariant 0 <= rem_nemo < modulo diff --git a/test/groebner/groebner.jl b/test/groebner/groebner.jl index 11233bce..47e39c1f 100644 --- a/test/groebner/groebner.jl +++ b/test/groebner/groebner.jl @@ -193,10 +193,9 @@ end gb = Groebner.groebner(noon) @test Groebner.isgroebner(gb) && all(iszero, Groebner.normalform(gb, noon)) - # TODO - # trace, gb00 = Groebner.groebner_learn(gb) - # flag, gb01 = Groebner.groebner_apply!(trace, gb) - # @test gb == gb00 == gb01 && flag + trace, gb00 = Groebner.groebner_learn(gb) + flag, gb01 = Groebner.groebner_apply!(trace, gb) + @test gb == gb00 == gb01 && flag end # Larger fields @@ -211,9 +210,9 @@ end gb = Groebner.groebner(noon) @test Groebner.isgroebner(gb) && all(iszero, Groebner.normalform(gb, noon)) - # trace, gb00 = Groebner.groebner_learn(gb) - # flag, gb01 = Groebner.groebner_apply!(trace, gb) - # @test gb == gb00 == gb01 && flag + trace, gb00 = Groebner.groebner_learn(gb) + flag, gb01 = Groebner.groebner_apply!(trace, gb) + @test gb == gb00 == gb01 && flag end # Smaller fields @@ -227,9 +226,9 @@ end gb = Groebner.groebner(noon) @test Groebner.isgroebner(gb) && all(iszero, Groebner.normalform(gb, noon)) - # trace, gb00 = Groebner.groebner_learn(gb) - # flag, gb01 = Groebner.groebner_apply!(trace, gb) - # @test gb == gb00 == gb01 && flag + trace, gb00 = Groebner.groebner_learn(gb) + flag, gb01 = Groebner.groebner_apply!(trace, gb) + @test gb == gb00 == gb01 && flag end # Tiny fields @@ -243,9 +242,9 @@ end gb = Groebner.groebner(noon) @test Groebner.isgroebner(gb) && all(iszero, Groebner.normalform(gb, noon)) - # trace, gb00 = Groebner.groebner_learn(gb) - # flag, gb01 = Groebner.groebner_apply!(trace, gb) - # @test gb == gb00 == gb01 && flag + trace, gb00 = Groebner.groebner_learn(gb) + flag, gb01 = Groebner.groebner_apply!(trace, gb) + @test gb == gb00 == gb01 && flag end end @@ -363,9 +362,6 @@ end fs = [(12345678 // 12347)x, (222222221111123 // 2131232232097)y + z] G = Groebner.groebner(fs) @test G == [y + 2131232232097 // 222222221111123 * z, x] - - # TODO TODO TODO: infinite loop - # groebner([x^2 + (2^31 - 1)*x + 1, y]) end @testset "groebner output sorted" begin @@ -907,9 +903,8 @@ end gb = Groebner.groebner(system) @test gb == result - # TODO: Sasha is too greedy to support this case. - # system = [x1 - (2^31 - 1) * x2 - (2^30 + 3) * x3] - # @test Groebner.groebner(system) == system + system = [x1 - (2^31 - 1) * x2 - (2^30 + 3) * x3] + @test Groebner.groebner(system) == system for start_of_range in [2^10, 2^20, 2^30, 2^40] for size_of_range in [10, 100, 1000] From 472053a8f34c9dcb609a27c4def95a9f0c3ccab7 Mon Sep 17 00:00:00 2001 From: Sasha Demin Date: Fri, 20 Sep 2024 23:42:01 +0200 Subject: [PATCH 4/6] cleanup test directory a bit --- src/Groebner.jl | 2 +- src/groebner/groebner.jl | 70 ++-- src/groebner/groebner_with_change_matrix.jl | 2 +- src/groebner/isgroebner.jl | 22 +- src/groebner/parameters.jl | 8 - test/{arithmetic/Zp.jl => arithmetic.jl} | 0 test/{groebner => }/groebner.jl | 314 ++++++++++++++++++ test/groebner/groebner_stress.jl | 105 ------ test/groebner/groebner_with_change_matrix.jl | 41 --- test/groebner/homogenization.jl | 100 ------ test/groebner/multi_threading.jl | 74 ----- test/{isgroebner => }/isgroebner.jl | 0 test/{learn_and_apply => }/learn_and_apply.jl | 89 +++++ test/learn_and_apply/apply_in_batches.jl | 91 ----- .../normalform_stress.jl => normalform.jl} | 92 ++++- test/normalform/normalform.jl | 89 ----- test/runtests.jl | 27 +- 17 files changed, 547 insertions(+), 579 deletions(-) rename test/{arithmetic/Zp.jl => arithmetic.jl} (100%) rename test/{groebner => }/groebner.jl (75%) delete mode 100644 test/groebner/groebner_stress.jl delete mode 100644 test/groebner/groebner_with_change_matrix.jl delete mode 100644 test/groebner/homogenization.jl delete mode 100644 test/groebner/multi_threading.jl rename test/{isgroebner => }/isgroebner.jl (100%) rename test/{learn_and_apply => }/learn_and_apply.jl (84%) delete mode 100644 test/learn_and_apply/apply_in_batches.jl rename test/{normalform/normalform_stress.jl => normalform.jl} (54%) delete mode 100644 test/normalform/normalform.jl diff --git a/src/Groebner.jl b/src/Groebner.jl index 8f264cfb..177ec1c9 100644 --- a/src/Groebner.jl +++ b/src/Groebner.jl @@ -150,7 +150,7 @@ include("groebner/homogenization.jl") include("interface.jl") using PrecompileTools -# include("precompile.jl") +include("precompile.jl") ### # Exports diff --git a/src/groebner/groebner.jl b/src/groebner/groebner.jl index 377dccf1..76b99bfb 100644 --- a/src/groebner/groebner.jl +++ b/src/groebner/groebner.jl @@ -240,7 +240,7 @@ function _groebner_learn_and_apply( # At this point, either the reconstruction or the correctness check failed. # Continue to compute Groebner bases modulo different primes in batches. primes_used = 1 - batchsize = 1 + batchsize = 4 batchsize_scaling = 0.10 witness_set = modular_witness_set(state.gb_coeffs_zz, params) @@ -250,52 +250,30 @@ function _groebner_learn_and_apply( iters = 0 while !correct_basis - if iszero(batchsize % 4) && params.batched - for j in 1:4:batchsize - prime_4x = ntuple(i -> Int32(modular_next_prime!(state)), 4) - - # Perform reduction modulo primes and store result in basis_ff_4x - ring_ff_4x, basis_ff_4x = modular_reduce_mod_p_in_batch!(ring, basis_zz, prime_4x) - params_zp_4x = params_mod_p( - params, - CompositeNumber{4, Int32}(prime_4x), - using_wide_type_for_coeffs=false - ) - trace_4x.buf_basis = basis_ff_4x - trace_4x.ring = ring_ff_4x - - f4_apply!(trace_4x, ring_ff_4x, trace_4x.buf_basis, params_zp_4x) - gb_coeffs_1, gb_coeffs_2, gb_coeffs_3, gb_coeffs_4 = - ir_unpack_composite_coefficients(trace_4x.gb_basis.coeffs) - - # TODO: This causes unnecessary conversions of arrays. - append!(state.used_primes, prime_4x) - push!(state.gb_coeffs_ff_all, gb_coeffs_1) - push!(state.gb_coeffs_ff_all, gb_coeffs_2) - push!(state.gb_coeffs_ff_all, gb_coeffs_3) - push!(state.gb_coeffs_ff_all, gb_coeffs_4) - primes_used += 4 - end - else - for j in 1:batchsize - prime = modular_next_prime!(state) - - ring_ff, basis_ff = modular_reduce_mod_p!(ring, basis_zz, prime, deepcopy=true) - params_zp = params_mod_p(params, prime) - - trace.buf_basis = basis_ff - trace.ring = ring_ff - - f4_apply!(trace, ring_ff, trace.buf_basis, params_zp) - - push!(state.used_primes, prime) - push!(state.gb_coeffs_ff_all, deepcopy(trace.gb_basis.coeffs)) + for j in 1:4:batchsize + prime_4x = ntuple(i -> Int32(modular_next_prime!(state)), 4) + + # Perform reduction modulo primes and store result in basis_ff_4x + ring_ff_4x, basis_ff_4x = modular_reduce_mod_p_in_batch!(ring, basis_zz, prime_4x) + params_zp_4x = params_mod_p( + params, + CompositeNumber{4, Int32}(prime_4x), + using_wide_type_for_coeffs=false + ) + trace_4x.buf_basis = basis_ff_4x + trace_4x.ring = ring_ff_4x - if !modular_majority_vote!(state, trace.gb_basis, params) - continue - end - primes_used += 1 - end + f4_apply!(trace_4x, ring_ff_4x, trace_4x.buf_basis, params_zp_4x) + gb_coeffs_1, gb_coeffs_2, gb_coeffs_3, gb_coeffs_4 = + ir_unpack_composite_coefficients(trace_4x.gb_basis.coeffs) + + # TODO: This causes unnecessary conversions of arrays. + append!(state.used_primes, prime_4x) + push!(state.gb_coeffs_ff_all, gb_coeffs_1) + push!(state.gb_coeffs_ff_all, gb_coeffs_2) + push!(state.gb_coeffs_ff_all, gb_coeffs_3) + push!(state.gb_coeffs_ff_all, gb_coeffs_4) + primes_used += 4 end crt_vec_partial!( diff --git a/src/groebner/groebner_with_change_matrix.jl b/src/groebner/groebner_with_change_matrix.jl index 0939eec7..ff0b4ca5 100644 --- a/src/groebner/groebner_with_change_matrix.jl +++ b/src/groebner/groebner_with_change_matrix.jl @@ -147,7 +147,7 @@ function _groebner_with_change_classic_modular( state = ModularState{BigInt, C, CoeffModular}(basis_zz.coeffs) - prime = modular_next_prime!(state) + prime = _groebner_guess_lucky_prime(state, ring, basis_zz, pairset, hashtable, params) # Perform reduction modulo prime and store result in basis_ff ring_ff, basis_ff = modular_reduce_mod_p!(ring, basis_zz, prime, deepcopy=true) diff --git a/src/groebner/isgroebner.jl b/src/groebner/isgroebner.jl index d23e6309..6ee5ff5e 100644 --- a/src/groebner/isgroebner.jl +++ b/src/groebner/isgroebner.jl @@ -70,12 +70,28 @@ function _isgroebner2( flag = f4_isgroebner!(ring, basis, pairset, hashtable, params.arithmetic) return flag end - # Otherwise, check modulo a prime + # Otherwise, check modulo primes basis_zz = clear_denominators!(basis, deepcopy=false) state = ModularState{BigInt, Rational{BigInt}, UInt32}(basis_zz.coeffs) prime = modular_random_prime(state, params.rng) ring_ff, basis_ff = modular_reduce_mod_p!(ring, basis_zz, prime, deepcopy=true) arithmetic = select_arithmetic(CoeffModular, prime, :auto, false) - flag = f4_isgroebner!(ring_ff, basis_ff, pairset, hashtable, arithmetic) - flag + flag1 = f4_isgroebner!(ring_ff, basis_ff, pairset, hashtable, arithmetic) + prime = modular_random_prime(state, params.rng) + ring_ff, basis_ff = modular_reduce_mod_p!(ring, basis_zz, prime, deepcopy=true) + arithmetic = select_arithmetic(CoeffModular, prime, :auto, false) + flag2 = f4_isgroebner!(ring_ff, basis_ff, pairset, hashtable, arithmetic) + if flag1 == flag2 + return flag1 + end + prime = modular_random_prime(state, params.rng) + ring_ff, basis_ff = modular_reduce_mod_p!(ring, basis_zz, prime, deepcopy=true) + arithmetic = select_arithmetic(CoeffModular, prime, :auto, false) + flag3 = f4_isgroebner!(ring_ff, basis_ff, pairset, hashtable, arithmetic) + if flag1 == flag3 + return flag1 + else + return flag2 + end + flag1 end diff --git a/src/groebner/parameters.jl b/src/groebner/parameters.jl index 4987637b..f4a55ae4 100644 --- a/src/groebner/parameters.jl +++ b/src/groebner/parameters.jl @@ -200,10 +200,6 @@ mutable struct AlgorithmParameters{Arithmetic <: AbstractArithmetic} # If learn & apply strategy can use apply in batches batched::Bool - # In multi-modular computation, compute (at least!) this many bases modulo - # different primes until a consensus in is reached - majority_threshold::Int - # Multi-threading threaded_f4::Symbol threaded_multimodular::Symbol @@ -316,8 +312,6 @@ function AlgorithmParameters(ring::PolyRing, kwargs::KeywordArguments; hint=:non end batched = kwargs.batched - majority_threshold = 1 - seed = kwargs.seed rng = Random.Xoshiro(seed) @@ -342,7 +336,6 @@ function AlgorithmParameters(ring::PolyRing, kwargs::KeywordArguments; hint=:non reduced, modular_strategy, batched, - majority_threshold, threaded_f4, threaded_multimodular, rng, @@ -379,7 +372,6 @@ function params_mod_p( params.reduced, params.modular_strategy, params.batched, - params.majority_threshold, params.threaded_f4, params.threaded_multimodular, params.rng, diff --git a/test/arithmetic/Zp.jl b/test/arithmetic.jl similarity index 100% rename from test/arithmetic/Zp.jl rename to test/arithmetic.jl diff --git a/test/groebner/groebner.jl b/test/groebner.jl similarity index 75% rename from test/groebner/groebner.jl rename to test/groebner.jl index 47e39c1f..75ff6e87 100644 --- a/test/groebner/groebner.jl +++ b/test/groebner.jl @@ -1,4 +1,5 @@ using AbstractAlgebra +using Base.Threads using Combinatorics, Primes, Random using Test, TestSetExtensions @@ -1032,3 +1033,316 @@ end gb = Groebner.groebner(f) end end + +@testset "homogenization, basic" begin + # TODO: broken for char. 113 + for field in [GF(2^10 + 7), GF(2^62 + 135), QQ] + R, x = polynomial_ring(field, "x") + @test Groebner.groebner([x^2 - 1, (x + 1) * x^2], homogenize=:yes) == [x + 1] + @test Groebner.groebner([x^2 - 1, (x + 1) * x^2], homogenize=:no) == [x + 1] + @test Groebner.groebner([x^2 - 1, (x + 1) * x^2], homogenize=:auto) == [x + 1] + + for ordering in [:lex, :deglex, :degrevlex] + R, (x, y) = polynomial_ring(field, ["x", "y"], internal_ordering=ordering) + + @test Groebner.groebner([R(3)], homogenize=:yes) == [R(1)] + @test Groebner.groebner([R(0)], homogenize=:yes) == [R(0)] + @test Groebner.groebner([x], homogenize=:yes) == [x] + @test Groebner.groebner([x + 11y], homogenize=:yes) == [x + 11y] + @test Groebner.groebner([x, y], homogenize=:yes) == [y, x] + @test Groebner.isgroebner(Groebner.groebner([x + 1, y + 2], homogenize=:yes)) + + if ordering === :lex + gb = Groebner.groebner([x + y^2, x^2 - 1], homogenize=:yes) + @test (y^4 - 1) in gb && (x + y^2) in gb + end + + for case in [ + Groebner.Examples.katsuran(2, k=field, internal_ordering=ordering), + Groebner.Examples.katsuran(3, k=field, internal_ordering=ordering), + Groebner.Examples.katsuran(4, k=field, internal_ordering=ordering), + Groebner.Examples.noonn(2, k=field, internal_ordering=ordering), + Groebner.Examples.noonn(3, k=field, internal_ordering=ordering) + ] + gb = Groebner.groebner(case, homogenize=:yes) + @test Groebner.isgroebner(gb) + end + end + end +end + +@testset "homogenization, orderings" begin + for field in [GF(113), GF(2^62 + 135), QQ] + # Test that the basis obtained with the use of homogenization + # *coincides* with the one obtained without it + R, (x, y, z) = polynomial_ring(field, ["x", "y", "z"]) + for case in [ + (system=[R(1)], ord=Groebner.Lex()), + (system=[R(5)], ord=Groebner.DegRevLex()), + (system=[x, y, z], ord=Groebner.Lex()), + (system=[x, y, z], ord=Groebner.Lex(z) * Groebner.Lex(x, y)), + (system=[x^2 + 1, x * y + 2, y * z + 3], ord=Groebner.DegLex()), + ( + system=[x^20 - x^15 - x^5 + y * z, x * y^10 + x^3 * y^3 + x * y], + ord=Groebner.DegRevLex() + ), + (system=[x + 5y, x + 7z, y + 11z], ord=Groebner.Lex(z) * Groebner.Lex(x, y)), + (system=[x + 5y + 11z], ord=Groebner.DegRevLex(y) * Groebner.Lex(z, x)), + (system=[x + 5y + 11z], ord=Groebner.DegRevLex(z) * Groebner.DegRevLex(y, x)), + (system=[x * y^2 + x + 1, y * z^2 + y + 1, z^4 - z^2 - 1], ord=Groebner.Lex()), + (system=[x * y^2 + x + 1, y * z^2 + y + 1, z^4 - z^2 - 1], ord=Groebner.DegRevLex()), + ( + system=[x * y^2 + x + 1, y * z^2 + y + 1, z^4 - z^2 - 1], + ord=Groebner.DegRevLex(z) * Groebner.DegRevLex(x, y) + ), + ( + system=[x * y^2 + x + 1, y * z^2 + y + 1, x * y * z^4 - x * z^2 - y + z], + ord=Groebner.DegRevLex(y) * Groebner.DegRevLex(z, x) + ), + (system=Groebner.Examples.katsuran(4, k=field), ord=Groebner.DegRevLex()), + (system=Groebner.Examples.rootn(5, k=field), ord=Groebner.Lex()), + (system=Groebner.Examples.reimern(5, k=field), ord=Groebner.DegRevLex()) + ] + ord = case.ord + system = case.system + @debug "" system ord field + + gb1 = Groebner.groebner(system, ordering=ord, homogenize=:no) + gb2 = Groebner.groebner(system, ordering=ord, homogenize=:yes) + + @test Groebner.isgroebner(gb1, ordering=ord) + + @test gb1 == gb2 + + # Also test learn / apply + field == QQ && continue + + context3, gb3 = + Groebner.groebner_learn(system, ordering=ord, homogenize=:yes) + context4, gb4 = + Groebner.groebner_learn(system, ordering=ord, homogenize=:no) + for _ in 1:4 + flag3, gb33 = Groebner.groebner_apply!(context3, system) + flag4, gb44 = Groebner.groebner_apply!(context4, system) + @test flag3 && flag4 + @test gb1 == gb3 == gb4 == gb33 == gb44 + end + end + end +end + +@testset "groebner, change matrix" begin + R, (x, y, z) = polynomial_ring(GF(2^30 + 3), ["x", "y", "z"], internal_ordering=:degrevlex) + f = [x * y * z - 1, x * y + x * z + y * z, x + y + z] + g, m = Groebner.groebner_with_change_matrix(f) + @test m * f == g + @test size(m) == (length(g), length(f)) + @test Groebner.isgroebner(g) + + for ground in [GF(2^30 + 3), QQ] + R, (x, y) = polynomial_ring(ground, ["x", "y"], internal_ordering=:degrevlex) + + cases = [ + [R(0), R(0), R(0)], + [R(0), R(1), R(0)], + [x], + [y], + [x, x, x, x, R(0), R(0), x, x, x, x], + [x + y, R(1), x^5 - y^5], + [x^200 + y^250, x^100 * y^200 + 1], + Groebner.Examples.katsuran(4, k=ground), + Groebner.Examples.noonn(4, k=ground), + Groebner.Examples.cyclicn(4, k=ground), + [x^i * y + x for i in 1:500] + ] + + for f in cases + g, m = Groebner.groebner_with_change_matrix(f) + @test m * f == g + @test size(m) == (length(g), length(f)) + @test Groebner.isgroebner(g) + end + end + + f = Groebner.Examples.katsuran(4, k=QQ, internal_ordering=:lex) + g, m = Groebner.groebner_with_change_matrix(f, ordering=Groebner.DegRevLex()) + @test m * f == g + @test_throws DomainError Groebner.groebner_with_change_matrix(f, ordering=Groebner.DegLex()) + @test_throws DomainError Groebner.groebner_with_change_matrix(f) +end + +@testset "groebner, multi-threading, Zp" begin + @info "Testing multi-threading over Zp using $(nthreads()) threads" + + R, (x, y, z) = GF(2^31 - 1)["x", "y", "z"] + + s = [R(1), R(2), R(4)] + @test Groebner.groebner(s, threaded=:yes) == + Groebner.groebner(s, threaded=:no) == + Groebner.groebner(s, threaded=:auto) == + [R(1)] + + s = [x^(2^10) + 1, y^3000 + 2, z^1000 + 3] + @test Groebner.groebner(s, threaded=:yes) == + Groebner.groebner(s, threaded=:no) == + Groebner.groebner(s, threaded=:auto) == + [z^1000 + 3, y^3000 + 2, x^(2^10) + 1] + + linalg = [:deterministic, :randomized] + threaded = [:yes, :no, :auto] + grid = [(linalg=l, threaded=t) for l in linalg for t in threaded] + for system in [ + Groebner.Examples.katsuran(3, internal_ordering=:lex, k=GF(2^31 - 1)), + Groebner.Examples.katsuran(4, internal_ordering=:lex, k=GF(2^20 + 7)), + Groebner.Examples.katsuran(7, internal_ordering=:degrevlex, k=GF(2^31 - 1)), + Groebner.Examples.eco5(internal_ordering=:deglex, k=GF(2^20 + 7)), + Groebner.Examples.eco5(internal_ordering=:degrevlex, k=GF(2^20 + 7)), + Groebner.Examples.cyclicn(5, internal_ordering=:degrevlex, k=GF(2^20 + 7)), + Groebner.Examples.cyclicn(6, internal_ordering=:degrevlex, k=GF(2^20 + 7)), + Groebner.Examples.ojika4(internal_ordering=:lex, k=GF(2^20 + 7)), + Groebner.Examples.henrion5(internal_ordering=:degrevlex, k=GF(2^20 + 7)), + Groebner.Examples.noonn(6, internal_ordering=:degrevlex, k=GF(2^10 + 7)), + Groebner.Examples.ku10(internal_ordering=:degrevlex, k=GF(2^10 + 7)) + ] + results = Array{Any}(undef, size(grid)) + for (i, kw) in enumerate(grid) + gb = Groebner.groebner(system; kw...) + results[i] = gb + end + @test length(unique(results)) == 1 + @test Groebner.isgroebner(results[1]) + end +end + +@testset "groebner, multi-threading, QQ" begin + @info "Testing multi-threading over QQ using $(nthreads()) threads" + + R, (x, y) = QQ["x", "y"] + + threaded = [:yes, :no, :auto] + grid = [(threaded=t,) for t in threaded] + for system in [ + [x - 1, y - 2], + [x + (BigInt(2)^1000 + 1) // 2^61, x * y + BigInt(2)^(2^10)], + Groebner.Examples.katsuran(3, internal_ordering=:lex, k=QQ), + Groebner.Examples.katsuran(4, internal_ordering=:lex, k=QQ), + Groebner.Examples.eco5(internal_ordering=:degrevlex, k=QQ), + Groebner.Examples.cyclicn(5, internal_ordering=:degrevlex, k=QQ), + Groebner.Examples.ojika4(internal_ordering=:lex, k=QQ), + Groebner.Examples.noonn(6, internal_ordering=:degrevlex, k=QQ), + Groebner.Examples.ku10(internal_ordering=:degrevlex, k=QQ), + [x + BigInt(2)^1_000 // (BigInt(2)^1_000 - 1), y - BigInt(2)^1_000], + [x + BigInt(2)^10_000 // (BigInt(2)^10_000 - 1), y^2 - BigInt(2)^10_000 * y] + ] + results = Array{Any}(undef, size(grid)) + for (i, kw) in enumerate(grid) + gb = Groebner.groebner(system; kw...) + results[i] = gb + end + @test length(unique(results)) == 1 + @test Groebner.isgroebner(results[1]) + end +end + +function test_params( + rng, + nvariables, + maxdegs, + nterms, + npolys, + grounds, + coeffssize, + orderings, + linalgs, + monoms, + homogenizes +) + boot = 1 + for _ in 1:boot + for nv in nvariables + for md in maxdegs + for nt in nterms + for np in npolys + for k in grounds + for ord in orderings + for cf in coeffssize + set = Groebner.Examples.random_generating_set( + rng, + k, + ord, + nv, + md, + nt, + np, + cf + ) + isempty(set) && continue + + for linalg in linalgs + for monom in monoms + for homogenize in homogenizes + try + gb = Groebner.groebner( + set, + linalg=linalg, + monoms=monom, + homogenize=homogenize + ) + flag = Groebner.isgroebner(gb) + if !flag + @error "Beda!" nv md nt np k ord monom + println("Rng:\n", rng) + println("Set:\n", set) + println("Gb:\n", gb) + end + @test flag + catch err + @error "Beda!" nv md nt np k ord monom + println(err) + println("Rng:\n", rng) + println("Set:\n", set) + rethrow(err) + end + end + end + end + end + end + end + end + end + end + end + end +end + +@testset "groebner random stress tests" begin + rng = Random.MersenneTwister(42) + + nvariables = [2, 3] + maxdegs = [2, 4] + nterms = [1, 2, 4] + npolys = [1, 4, 100] + grounds = [GF(1031), GF(2^50 + 55), AbstractAlgebra.QQ] + coeffssize = [3, 1000, 2^31 - 1, BigInt(2)^100] + orderings = [:degrevlex, :lex, :deglex] + linalgs = [:deterministic, :randomized] + monoms = [:auto, :dense, :packed] + homogenize = [:yes, :auto] + p = prod(map(length, (nvariables, maxdegs, nterms, npolys, grounds, orderings, coeffssize, linalgs, monoms, homogenize))) + @info "Producing $p small random tests for groebner. This may take a minute" + + test_params( + rng, + nvariables, + maxdegs, + nterms, + npolys, + grounds, + coeffssize, + orderings, + linalgs, + monoms, + homogenize + ) +end diff --git a/test/groebner/groebner_stress.jl b/test/groebner/groebner_stress.jl deleted file mode 100644 index ad541d4c..00000000 --- a/test/groebner/groebner_stress.jl +++ /dev/null @@ -1,105 +0,0 @@ -import Random -using AbstractAlgebra - -function test_params( - rng, - nvariables, - maxdegs, - nterms, - npolys, - grounds, - coeffssize, - orderings, - linalgs, - monoms, - homogenizes -) - boot = 1 - for _ in 1:boot - for nv in nvariables - for md in maxdegs - for nt in nterms - for np in npolys - for k in grounds - for ord in orderings - for cf in coeffssize - set = Groebner.Examples.random_generating_set( - rng, - k, - ord, - nv, - md, - nt, - np, - cf - ) - isempty(set) && continue - - for linalg in linalgs - for monom in monoms - for homogenize in homogenizes - try - gb = Groebner.groebner( - set, - linalg=linalg, - monoms=monom, - homogenize=homogenize - ) - flag = Groebner.isgroebner(gb) - if !flag - @error "Beda!" nv md nt np k ord monom - println("Rng:\n", rng) - println("Set:\n", set) - println("Gb:\n", gb) - end - @test flag - catch err - @error "Beda!" nv md nt np k ord monom - println(err) - println("Rng:\n", rng) - println("Set:\n", set) - rethrow(err) - end - end - end - end - end - end - end - end - end - end - end - end -end - -@testset "groebner random stress tests" begin - rng = Random.MersenneTwister(42) - - nvariables = [2, 3] - maxdegs = [2, 4] - nterms = [1, 2, 4] - npolys = [1, 4, 100] - grounds = [GF(1031), GF(2^50 + 55), AbstractAlgebra.QQ] - coeffssize = [3, 1000, 2^31 - 1, BigInt(2)^100] - orderings = [:degrevlex, :lex, :deglex] - linalgs = [:deterministic, :randomized] - monoms = [:auto, :dense, :packed] - homogenize = [:yes, :auto] - p = prod(map(length, (nvariables, maxdegs, nterms, npolys, grounds, orderings, coeffssize, linalgs, monoms, homogenize))) - @info "Producing $p small random tests for groebner. This may take a minute" - - test_params( - rng, - nvariables, - maxdegs, - nterms, - npolys, - grounds, - coeffssize, - orderings, - linalgs, - monoms, - homogenize - ) -end diff --git a/test/groebner/groebner_with_change_matrix.jl b/test/groebner/groebner_with_change_matrix.jl deleted file mode 100644 index 2ea8758f..00000000 --- a/test/groebner/groebner_with_change_matrix.jl +++ /dev/null @@ -1,41 +0,0 @@ -using AbstractAlgebra - -@testset "groebner, change matrix" begin - R, (x, y, z) = polynomial_ring(GF(2^30 + 3), ["x", "y", "z"], internal_ordering=:degrevlex) - f = [x * y * z - 1, x * y + x * z + y * z, x + y + z] - g, m = Groebner.groebner_with_change_matrix(f) - @test m * f == g - @test size(m) == (length(g), length(f)) - @test Groebner.isgroebner(g) - - for ground in [GF(2^30 + 3), QQ] - R, (x, y) = polynomial_ring(ground, ["x", "y"], internal_ordering=:degrevlex) - - cases = [ - [R(0), R(0), R(0)], - [R(0), R(1), R(0)], - [x], - [y], - [x, x, x, x, R(0), R(0), x, x, x, x], - [x + y, R(1), x^5 - y^5], - [x^200 + y^250, x^100 * y^200 + 1], - Groebner.Examples.katsuran(4, k=ground), - Groebner.Examples.noonn(4, k=ground), - Groebner.Examples.cyclicn(4, k=ground), - [x^i * y + x for i in 1:500] - ] - - for f in cases - g, m = Groebner.groebner_with_change_matrix(f) - @test m * f == g - @test size(m) == (length(g), length(f)) - @test Groebner.isgroebner(g) - end - end - - f = Groebner.Examples.katsuran(4, k=QQ, internal_ordering=:lex) - g, m = Groebner.groebner_with_change_matrix(f, ordering=Groebner.DegRevLex()) - @test m * f == g - @test_throws DomainError Groebner.groebner_with_change_matrix(f, ordering=Groebner.DegLex()) - @test_throws DomainError Groebner.groebner_with_change_matrix(f) -end diff --git a/test/groebner/homogenization.jl b/test/groebner/homogenization.jl deleted file mode 100644 index 3d858261..00000000 --- a/test/groebner/homogenization.jl +++ /dev/null @@ -1,100 +0,0 @@ -using AbstractAlgebra -using Combinatorics -import Random - -@testset "homogenization, basic" begin - # TODO: broken for char. 113 - for field in [GF(2^10 + 7), GF(2^62 + 135), QQ] - R, x = polynomial_ring(field, "x") - @test Groebner.groebner([x^2 - 1, (x + 1) * x^2], homogenize=:yes) == [x + 1] - @test Groebner.groebner([x^2 - 1, (x + 1) * x^2], homogenize=:no) == [x + 1] - @test Groebner.groebner([x^2 - 1, (x + 1) * x^2], homogenize=:auto) == [x + 1] - - for ordering in [:lex, :deglex, :degrevlex] - R, (x, y) = polynomial_ring(field, ["x", "y"], internal_ordering=ordering) - - @test Groebner.groebner([R(3)], homogenize=:yes) == [R(1)] - @test Groebner.groebner([R(0)], homogenize=:yes) == [R(0)] - @test Groebner.groebner([x], homogenize=:yes) == [x] - @test Groebner.groebner([x + 11y], homogenize=:yes) == [x + 11y] - @test Groebner.groebner([x, y], homogenize=:yes) == [y, x] - @test Groebner.isgroebner(Groebner.groebner([x + 1, y + 2], homogenize=:yes)) - - if ordering === :lex - gb = Groebner.groebner([x + y^2, x^2 - 1], homogenize=:yes) - @test (y^4 - 1) in gb && (x + y^2) in gb - end - - for case in [ - Groebner.Examples.katsuran(2, k=field, internal_ordering=ordering), - Groebner.Examples.katsuran(3, k=field, internal_ordering=ordering), - Groebner.Examples.katsuran(4, k=field, internal_ordering=ordering), - Groebner.Examples.noonn(2, k=field, internal_ordering=ordering), - Groebner.Examples.noonn(3, k=field, internal_ordering=ordering) - ] - gb = Groebner.groebner(case, homogenize=:yes) - @test Groebner.isgroebner(gb) - end - end - end -end - -@testset "homogenization, orderings" begin - for field in [GF(113), GF(2^62 + 135), QQ] - # Test that the basis obtained with the use of homogenization - # *coincides* with the one obtained without it - R, (x, y, z) = polynomial_ring(field, ["x", "y", "z"]) - for case in [ - (system=[R(1)], ord=Groebner.Lex()), - (system=[R(5)], ord=Groebner.DegRevLex()), - (system=[x, y, z], ord=Groebner.Lex()), - (system=[x, y, z], ord=Groebner.Lex(z) * Groebner.Lex(x, y)), - (system=[x^2 + 1, x * y + 2, y * z + 3], ord=Groebner.DegLex()), - ( - system=[x^20 - x^15 - x^5 + y * z, x * y^10 + x^3 * y^3 + x * y], - ord=Groebner.DegRevLex() - ), - (system=[x + 5y, x + 7z, y + 11z], ord=Groebner.Lex(z) * Groebner.Lex(x, y)), - (system=[x + 5y + 11z], ord=Groebner.DegRevLex(y) * Groebner.Lex(z, x)), - (system=[x + 5y + 11z], ord=Groebner.DegRevLex(z) * Groebner.DegRevLex(y, x)), - (system=[x * y^2 + x + 1, y * z^2 + y + 1, z^4 - z^2 - 1], ord=Groebner.Lex()), - (system=[x * y^2 + x + 1, y * z^2 + y + 1, z^4 - z^2 - 1], ord=Groebner.DegRevLex()), - ( - system=[x * y^2 + x + 1, y * z^2 + y + 1, z^4 - z^2 - 1], - ord=Groebner.DegRevLex(z) * Groebner.DegRevLex(x, y) - ), - ( - system=[x * y^2 + x + 1, y * z^2 + y + 1, x * y * z^4 - x * z^2 - y + z], - ord=Groebner.DegRevLex(y) * Groebner.DegRevLex(z, x) - ), - (system=Groebner.Examples.katsuran(4, k=field), ord=Groebner.DegRevLex()), - (system=Groebner.Examples.rootn(5, k=field), ord=Groebner.Lex()), - (system=Groebner.Examples.reimern(5, k=field), ord=Groebner.DegRevLex()) - ] - ord = case.ord - system = case.system - @debug "" system ord field - - gb1 = Groebner.groebner(system, ordering=ord, homogenize=:no) - gb2 = Groebner.groebner(system, ordering=ord, homogenize=:yes) - - @test Groebner.isgroebner(gb1, ordering=ord) - - @test gb1 == gb2 - - # Also test learn / apply - # if field != QQ - # context3, gb3 = - # Groebner.groebner_learn(system, ordering=ord, homogenize=:yes) - # context4, gb4 = - # Groebner.groebner_learn(system, ordering=ord, homogenize=:no) - # for _ in 1:4 - # flag3, gb33 = Groebner.groebner_apply!(context3, system) - # flag4, gb44 = Groebner.groebner_apply!(context4, system) - # @test flag3 && flag4 - # @test gb1 == gb3 == gb4 == gb33 == gb44 - # end - # end - end - end -end diff --git a/test/groebner/multi_threading.jl b/test/groebner/multi_threading.jl deleted file mode 100644 index e592e3ef..00000000 --- a/test/groebner/multi_threading.jl +++ /dev/null @@ -1,74 +0,0 @@ -using Base.Threads - -@testset "groebner, multi-threading, Zp" begin - @info "Testing multi-threading over Zp using $(nthreads()) threads" - - R, (x, y, z) = GF(2^31 - 1)["x", "y", "z"] - - s = [R(1), R(2), R(4)] - @test Groebner.groebner(s, threaded=:yes) == - Groebner.groebner(s, threaded=:no) == - Groebner.groebner(s, threaded=:auto) == - [R(1)] - - s = [x^(2^10) + 1, y^3000 + 2, z^1000 + 3] - @test Groebner.groebner(s, threaded=:yes) == - Groebner.groebner(s, threaded=:no) == - Groebner.groebner(s, threaded=:auto) == - [z^1000 + 3, y^3000 + 2, x^(2^10) + 1] - - linalg = [:deterministic, :randomized] - threaded = [:yes, :no, :auto] - grid = [(linalg=l, threaded=t) for l in linalg for t in threaded] - for system in [ - Groebner.Examples.katsuran(3, internal_ordering=:lex, k=GF(2^31 - 1)), - Groebner.Examples.katsuran(4, internal_ordering=:lex, k=GF(2^20 + 7)), - Groebner.Examples.katsuran(7, internal_ordering=:degrevlex, k=GF(2^31 - 1)), - Groebner.Examples.eco5(internal_ordering=:deglex, k=GF(2^20 + 7)), - Groebner.Examples.eco5(internal_ordering=:degrevlex, k=GF(2^20 + 7)), - Groebner.Examples.cyclicn(5, internal_ordering=:degrevlex, k=GF(2^20 + 7)), - Groebner.Examples.cyclicn(6, internal_ordering=:degrevlex, k=GF(2^20 + 7)), - Groebner.Examples.ojika4(internal_ordering=:lex, k=GF(2^20 + 7)), - Groebner.Examples.henrion5(internal_ordering=:degrevlex, k=GF(2^20 + 7)), - Groebner.Examples.noonn(6, internal_ordering=:degrevlex, k=GF(2^10 + 7)), - Groebner.Examples.ku10(internal_ordering=:degrevlex, k=GF(2^10 + 7)) - ] - results = Array{Any}(undef, size(grid)) - for (i, kw) in enumerate(grid) - gb = Groebner.groebner(system; kw...) - results[i] = gb - end - @test length(unique(results)) == 1 - @test Groebner.isgroebner(results[1]) - end -end - -@testset "groebner, multi-threading, QQ" begin - @info "Testing multi-threading over QQ using $(nthreads()) threads" - - R, (x, y) = QQ["x", "y"] - - threaded = [:yes, :no, :auto] - grid = [(threaded=t,) for t in threaded] - for system in [ - [x - 1, y - 2], - [x + (BigInt(2)^1000 + 1) // 2^61, x * y + BigInt(2)^(2^10)], - Groebner.Examples.katsuran(3, internal_ordering=:lex, k=QQ), - Groebner.Examples.katsuran(4, internal_ordering=:lex, k=QQ), - Groebner.Examples.eco5(internal_ordering=:degrevlex, k=QQ), - Groebner.Examples.cyclicn(5, internal_ordering=:degrevlex, k=QQ), - Groebner.Examples.ojika4(internal_ordering=:lex, k=QQ), - Groebner.Examples.noonn(6, internal_ordering=:degrevlex, k=QQ), - Groebner.Examples.ku10(internal_ordering=:degrevlex, k=QQ), - [x + BigInt(2)^1_000 // (BigInt(2)^1_000 - 1), y - BigInt(2)^1_000], - [x + BigInt(2)^10_000 // (BigInt(2)^10_000 - 1), y^2 - BigInt(2)^10_000 * y] - ] - results = Array{Any}(undef, size(grid)) - for (i, kw) in enumerate(grid) - gb = Groebner.groebner(system; kw...) - results[i] = gb - end - @test length(unique(results)) == 1 - @test Groebner.isgroebner(results[1]) - end -end diff --git a/test/isgroebner/isgroebner.jl b/test/isgroebner.jl similarity index 100% rename from test/isgroebner/isgroebner.jl rename to test/isgroebner.jl diff --git a/test/learn_and_apply/learn_and_apply.jl b/test/learn_and_apply.jl similarity index 84% rename from test/learn_and_apply/learn_and_apply.jl rename to test/learn_and_apply.jl index 425df76f..2216742f 100644 --- a/test/learn_and_apply/learn_and_apply.jl +++ b/test/learn_and_apply.jl @@ -494,3 +494,92 @@ end @info "Apply (expectedly) failed in $A / $B cases." end end + +@testset "learn & apply, in batches" begin + # Simple sanity checks + ks = map(GF, Primes.nextprimes(2^30, 40)) + + R, (x, y) = polynomial_ring(AbstractAlgebra.ZZ, ["x", "y"], internal_ordering=:degrevlex) + sys_qq = [44x^2 + x + 2^50, y^10 - 10 * y^5 - 99] + sys_gf = map( + j -> + map(poly -> AbstractAlgebra.map_coefficients(c -> (k = ks[j]; k(c)), poly), sys_qq), + 1:length(ks) + ) + + trace, gb = Groebner.groebner_learn(sys_gf[1]) + + flag1, gb1 = Groebner.groebner_apply!(trace, sys_gf[1]) + flag2, gb2 = Groebner.groebner_apply!(trace, sys_gf[2]) + flag3, (gb3, gb4) = Groebner.groebner_apply!(trace, (sys_gf[3], sys_gf[4])) + flag4, (gb5, gb6, gb7, gb8) = + Groebner.groebner_apply!(trace, (sys_gf[3], sys_gf[4], sys_gf[5], sys_gf[6])) + flag5, (gb9, gb10, gb11, gb12, gb13, gb14, gb15, gb16) = Groebner.groebner_apply!( + trace, + (sys_gf[3], sys_gf[4], sys_gf[5], sys_gf[6], sys_gf[7], sys_gf[8], sys_gf[9], sys_gf[10]) + ) + # he-he + flag6, gbs = Groebner.groebner_apply!(trace, (sys_gf[3:18]...,)) + + @test flag1 && flag2 && flag3 && flag4 && flag5 && flag6 + + true_gb = map(Groebner.groebner, sys_gf) + @test [gb1, gb2, gb3, gb4] == true_gb[1:4] + @test [gb5, gb6, gb7, gb8] == true_gb[3:6] + @test [gb9, gb10, gb11, gb12, gb13, gb14, gb15, gb16] == true_gb[3:10] + @test collect(gbs) == true_gb[3:18] + + function test_learn_apply(system, primes) + @assert length(primes) > 4 + systems_zp = + map(p -> map(f -> AbstractAlgebra.map_coefficients(c -> GF(p)(c), f), system), primes) + + trace, gb_0 = Groebner.groebner_learn(systems_zp[1]) + + flag, gb = Groebner.groebner_apply!(trace, systems_zp[2]) + @test gb == Groebner.groebner(systems_zp[2]) + + flag, (gb,) = Groebner.groebner_apply!(trace, (systems_zp[2],)) + @test gb == Groebner.groebner(systems_zp[2]) + + flag, gbs = Groebner.groebner_apply!(trace, (systems_zp[2:5]...,)) + @test collect(gbs) == map(Groebner.groebner, systems_zp[2:5]) + end + + kat = Groebner.Examples.katsuran(8, k=ZZ, internal_ordering=:degrevlex) + cyc = Groebner.Examples.cyclicn(6, k=ZZ, internal_ordering=:degrevlex) + cyc_lex = Groebner.Examples.cyclicn(5, k=ZZ, internal_ordering=:lex) + + ps1 = Primes.nextprimes(2^30, 5) + ps2 = [2^31 - 1, 2^30 + 3, 2^32 + 15, 2^40 + 15, 2^30 + 3] + + test_learn_apply(kat, ps1) + test_learn_apply(kat, ps2) + test_learn_apply(cyc, ps2) + test_learn_apply(cyc_lex, ps2) +end + +@testset "learn & apply low level, in batches" begin + ring0 = Groebner.PolyRing(2, Groebner.DegLex(), 5) + ring1 = Groebner.PolyRing(2, Groebner.DegLex(), 7) + ring2 = Groebner.PolyRing(2, Groebner.DegLex(), 11) + ring3 = Groebner.PolyRing(2, Groebner.DegLex(), 13) + ring4 = Groebner.PolyRing(2, Groebner.Lex(1) * Groebner.Lex(2), 17) + trace, (gb0...) = Groebner.groebner_learn(ring0, [[[0, 0], [1, 1]]], [[1, -1]]) + @test gb0 == ([[[1, 1], [0, 0]]], [[1, 4]]) + flag1, gb1 = Groebner.groebner_apply!(trace, ring1, [[[0, 0], [1, 1]]], [[1, -1]]) + flag2, gb2 = Groebner.groebner_apply!(trace, ring2, [[[1, 1], [0, 0]]], [[-1, 1]]) + flag3, gb3 = Groebner.groebner_apply!(trace, ring3, [[[0, 0], [0, 1], [1, 1]]], [[1, 0, -1]]) + flag4, gb4 = Groebner.groebner_apply!(trace, ring4, [[[0, 0], [1, 1]]], [[1, -1]]) + @test flag1 && flag2 && flag3 && flag4 + @test gb1 == ([[1, 6]]) + @test gb2 == ([[1, 10]]) + @test gb3 == ([[1, 12]]) + @test gb4 == ([[1, 16]]) + + flag1, gb1 = Groebner.groebner_apply!(trace, ring1, [[[0, 0]]], [[1]]) + flag2, gb2 = Groebner.groebner_apply!(trace, ring2, [[[1, 1], [0, 0]]], [[-1, 0]]) + flag3, gb3 = Groebner.groebner_apply!(trace, ring3, [[[0, 0], [0, 1], [1, 1]]], [[1, 0, -1]]) + flag4, (gb4...) = Groebner.groebner_apply!(trace, ring4, [[[0, 0], [1, 1]]], [[1, -1]]) + @test !flag1 && !flag2 && flag3 && flag4 +end diff --git a/test/learn_and_apply/apply_in_batches.jl b/test/learn_and_apply/apply_in_batches.jl deleted file mode 100644 index 53074312..00000000 --- a/test/learn_and_apply/apply_in_batches.jl +++ /dev/null @@ -1,91 +0,0 @@ -import Primes -using AbstractAlgebra, Test - -@testset "learn & apply, in batches" begin - # Simple sanity checks - ks = map(GF, Primes.nextprimes(2^30, 40)) - - R, (x, y) = polynomial_ring(AbstractAlgebra.ZZ, ["x", "y"], internal_ordering=:degrevlex) - sys_qq = [44x^2 + x + 2^50, y^10 - 10 * y^5 - 99] - sys_gf = map( - j -> - map(poly -> AbstractAlgebra.map_coefficients(c -> (k = ks[j]; k(c)), poly), sys_qq), - 1:length(ks) - ) - - trace, gb = Groebner.groebner_learn(sys_gf[1]) - - flag1, gb1 = Groebner.groebner_apply!(trace, sys_gf[1]) - flag2, gb2 = Groebner.groebner_apply!(trace, sys_gf[2]) - flag3, (gb3, gb4) = Groebner.groebner_apply!(trace, (sys_gf[3], sys_gf[4])) - flag4, (gb5, gb6, gb7, gb8) = - Groebner.groebner_apply!(trace, (sys_gf[3], sys_gf[4], sys_gf[5], sys_gf[6])) - flag5, (gb9, gb10, gb11, gb12, gb13, gb14, gb15, gb16) = Groebner.groebner_apply!( - trace, - (sys_gf[3], sys_gf[4], sys_gf[5], sys_gf[6], sys_gf[7], sys_gf[8], sys_gf[9], sys_gf[10]) - ) - # he-he - flag6, gbs = Groebner.groebner_apply!(trace, (sys_gf[3:18]...,)) - - @test flag1 && flag2 && flag3 && flag4 && flag5 && flag6 - - true_gb = map(Groebner.groebner, sys_gf) - @test [gb1, gb2, gb3, gb4] == true_gb[1:4] - @test [gb5, gb6, gb7, gb8] == true_gb[3:6] - @test [gb9, gb10, gb11, gb12, gb13, gb14, gb15, gb16] == true_gb[3:10] - @test collect(gbs) == true_gb[3:18] - - function test_learn_apply(system, primes) - @assert length(primes) > 4 - systems_zp = - map(p -> map(f -> AbstractAlgebra.map_coefficients(c -> GF(p)(c), f), system), primes) - - trace, gb_0 = Groebner.groebner_learn(systems_zp[1]) - - flag, gb = Groebner.groebner_apply!(trace, systems_zp[2]) - @test gb == Groebner.groebner(systems_zp[2]) - - flag, (gb,) = Groebner.groebner_apply!(trace, (systems_zp[2],)) - @test gb == Groebner.groebner(systems_zp[2]) - - flag, gbs = Groebner.groebner_apply!(trace, (systems_zp[2:5]...,)) - @test collect(gbs) == map(Groebner.groebner, systems_zp[2:5]) - end - - kat = Groebner.Examples.katsuran(8, k=ZZ, internal_ordering=:degrevlex) - cyc = Groebner.Examples.cyclicn(6, k=ZZ, internal_ordering=:degrevlex) - cyc_lex = Groebner.Examples.cyclicn(5, k=ZZ, internal_ordering=:lex) - - ps1 = Primes.nextprimes(2^30, 5) - ps2 = [2^31 - 1, 2^30 + 3, 2^32 + 15, 2^40 + 15, 2^30 + 3] - - test_learn_apply(kat, ps1) - test_learn_apply(kat, ps2) - test_learn_apply(cyc, ps2) - test_learn_apply(cyc_lex, ps2) -end - -@testset "learn & apply low level, in batches" begin - ring0 = Groebner.PolyRing(2, Groebner.DegLex(), 5) - ring1 = Groebner.PolyRing(2, Groebner.DegLex(), 7) - ring2 = Groebner.PolyRing(2, Groebner.DegLex(), 11) - ring3 = Groebner.PolyRing(2, Groebner.DegLex(), 13) - ring4 = Groebner.PolyRing(2, Groebner.Lex(1) * Groebner.Lex(2), 17) - trace, (gb0...) = Groebner.groebner_learn(ring0, [[[0, 0], [1, 1]]], [[1, -1]]) - @test gb0 == ([[[1, 1], [0, 0]]], [[1, 4]]) - flag1, gb1 = Groebner.groebner_apply!(trace, ring1, [[[0, 0], [1, 1]]], [[1, -1]]) - flag2, gb2 = Groebner.groebner_apply!(trace, ring2, [[[1, 1], [0, 0]]], [[-1, 1]]) - flag3, gb3 = Groebner.groebner_apply!(trace, ring3, [[[0, 0], [0, 1], [1, 1]]], [[1, 0, -1]]) - flag4, gb4 = Groebner.groebner_apply!(trace, ring4, [[[0, 0], [1, 1]]], [[1, -1]]) - @test flag1 && flag2 && flag3 && flag4 - @test gb1 == ([[1, 6]]) - @test gb2 == ([[1, 10]]) - @test gb3 == ([[1, 12]]) - @test gb4 == ([[1, 16]]) - - flag1, gb1 = Groebner.groebner_apply!(trace, ring1, [[[0, 0]]], [[1]]) - flag2, gb2 = Groebner.groebner_apply!(trace, ring2, [[[1, 1], [0, 0]]], [[-1, 0]]) - flag3, gb3 = Groebner.groebner_apply!(trace, ring3, [[[0, 0], [0, 1], [1, 1]]], [[1, 0, -1]]) - flag4, (gb4...) = Groebner.groebner_apply!(trace, ring4, [[[0, 0], [1, 1]]], [[1, -1]]) - @test !flag1 && !flag2 && flag3 && flag4 -end diff --git a/test/normalform/normalform_stress.jl b/test/normalform.jl similarity index 54% rename from test/normalform/normalform_stress.jl rename to test/normalform.jl index 428a63c1..a83666a8 100644 --- a/test/normalform/normalform_stress.jl +++ b/test/normalform.jl @@ -1,5 +1,93 @@ -import Random -using AbstractAlgebra +using AbstractAlgebra, Random + +@testset "normalform" begin + R, x = polynomial_ring(GF(2^31 - 1), "x") + + @test Groebner.normalform([x], R(0)) == R(0) + @test Groebner.normalform([x], R(1)) == R(1) + @test Groebner.normalform([x], [R(0)]) == [R(0)] + @test Groebner.normalform([x], [R(1)]) == [R(1)] + @test Groebner.normalform([x], [x, R(5), x + 1]) == [R(0), R(5), R(1)] + + R, (x, y, z) = polynomial_ring(GF(2^31 - 1), ["x", "y", "z"]) + + @test Groebner.normalform([x], R(0)) == R(0) + @test Groebner.normalform([x], R(1)) == R(1) + @test Groebner.normalform([x], [R(0)]) == [R(0)] + @test Groebner.normalform([x], [R(1)]) == [R(1)] + + Gs = [[x], [x, y], [x, y, z]] + + @test Groebner.normalform(Gs[1], x) == 0 + @test Groebner.normalform(Gs[1], y) == y + @test Groebner.normalform(Gs[3], x + y^2) == 0 + @test Groebner.normalform(Gs[2], 5x^2 - y^3 - 1) == -1 + @test Groebner.normalform(Gs[2], x + y + 3z) == 3z + + G = [x^2 + y * x + 1, y^2 - y] + + @test Groebner.normalform(G, x + y + z) == x + y + z + @test Groebner.normalform(G, x^2 + y^2) == y - y * x - 1 + @test Groebner.normalform(G, y^3) == Groebner.normalform(G, y^4) == y + + R, (x, y, z) = polynomial_ring(GF(2^31 - 1), ["x", "y", "z"], internal_ordering=:degrevlex) + G = [x^2 + y, y^2 + x] + @test Groebner.normalform(G, x^2 + y^2) == -x - y + + @test_nowarn Groebner.normalform([x, x * y], x) +end + +@testset "normalform many variables" begin + R, xs = polynomial_ring(QQ, ["x$i" for i in 1:63]) + GB = Groebner.groebner(xs) + @test GB == reverse(xs) + @test Groebner.normalform(GB, xs[1]) == R(0) + + R, xs = polynomial_ring(QQ, ["x$i" for i in 1:220]) + GB = Groebner.groebner(xs) + @test GB == reverse(xs) + @test Groebner.normalform(GB, xs[1]) == R(0) + + R, xs = polynomial_ring(QQ, ["x$i" for i in 1:1025]) + GB = Groebner.groebner(xs) + @test GB == reverse(xs) + @test Groebner.normalform(GB, xs) == zeros(R, length(xs)) +end + +@testset "normalform orderings" begin + R, (x, y, z) = polynomial_ring(QQ, ["x", "y", "z"]) + + @test Groebner.normalform([x + y], [x, y]) == [-y, y] + @test Groebner.normalform([x + y^2], [x, y], ordering=Groebner.DegLex()) == [x, y] + @test Groebner.normalform([x + y^2], [x, y], ordering=Groebner.DegRevLex()) == [x, y] + @test Groebner.normalform([x + y], [x, y], ordering=Groebner.Lex(y, x, z)) == [x, -x] +end + +@testset "normalform of an array" begin + for field in [GF(17), GF(2^31 - 1), QQ] + R, (x, y) = polynomial_ring(field, ["x", "y"]) + gb = [x, y] + @test Groebner.normalform(gb, [x, y + 1]) == [R(0), R(1)] + @test Groebner.normalform(gb, [y + 1, x]) == [R(1), R(0)] + @test Groebner.normalform(gb, [x, 3y, 4x + 5y, 6x * y, 7x + 1, y + 2, x^2 + y^2 + 3]) == + [R(0), R(0), R(0), R(0), R(1), R(2), R(3)] + + gb = [x] + @test Groebner.normalform(gb, [x, 3y, 4x + 5y, 6x * y, 7x + 1, y + 2, x^2 + y^2 + 3]) == + [R(0), 3y, 5y, R(0), R(1), y + 2, y^2 + 3] + + gb = [x^2 + y^2, y^3 - y] + @test Groebner.normalform( + gb, + [R(1), R(5), R(16), 3x, y, 4x * y + 2, y^2, x * y^2 - 8x + y] + ) == [R(1), R(5), R(16), 3x, y, 4x * y + 2, y^2, x * y^2 - 8x + y] + + gb = [x, R(0), y, R(0)] + @test Groebner.normalform(gb, [R(0)]) == [R(0)] + @test Groebner.normalform(gb, [R(0), R(1), R(0)]) == [R(0), R(1), R(0)] + @test Groebner.normalform([R(0), R(0)], [R(0), R(1), R(0)]) == [R(0), R(1), R(0)] + end +end function test_params_nf( rng, diff --git a/test/normalform/normalform.jl b/test/normalform/normalform.jl deleted file mode 100644 index fc8b2ebd..00000000 --- a/test/normalform/normalform.jl +++ /dev/null @@ -1,89 +0,0 @@ - -@testset "normalform" begin - R, x = polynomial_ring(GF(2^31 - 1), "x") - - @test Groebner.normalform([x], R(0)) == R(0) - @test Groebner.normalform([x], R(1)) == R(1) - @test Groebner.normalform([x], [R(0)]) == [R(0)] - @test Groebner.normalform([x], [R(1)]) == [R(1)] - @test Groebner.normalform([x], [x, R(5), x + 1]) == [R(0), R(5), R(1)] - - R, (x, y, z) = polynomial_ring(GF(2^31 - 1), ["x", "y", "z"]) - - @test Groebner.normalform([x], R(0)) == R(0) - @test Groebner.normalform([x], R(1)) == R(1) - @test Groebner.normalform([x], [R(0)]) == [R(0)] - @test Groebner.normalform([x], [R(1)]) == [R(1)] - - Gs = [[x], [x, y], [x, y, z]] - - @test Groebner.normalform(Gs[1], x) == 0 - @test Groebner.normalform(Gs[1], y) == y - @test Groebner.normalform(Gs[3], x + y^2) == 0 - @test Groebner.normalform(Gs[2], 5x^2 - y^3 - 1) == -1 - @test Groebner.normalform(Gs[2], x + y + 3z) == 3z - - G = [x^2 + y * x + 1, y^2 - y] - - @test Groebner.normalform(G, x + y + z) == x + y + z - @test Groebner.normalform(G, x^2 + y^2) == y - y * x - 1 - @test Groebner.normalform(G, y^3) == Groebner.normalform(G, y^4) == y - - R, (x, y, z) = polynomial_ring(GF(2^31 - 1), ["x", "y", "z"], internal_ordering=:degrevlex) - G = [x^2 + y, y^2 + x] - @test Groebner.normalform(G, x^2 + y^2) == -x - y - - @test_nowarn Groebner.normalform([x, x * y], x) -end - -@testset "normalform many variables" begin - R, xs = polynomial_ring(QQ, ["x$i" for i in 1:63]) - GB = Groebner.groebner(xs) - @test GB == reverse(xs) - @test Groebner.normalform(GB, xs[1]) == R(0) - - R, xs = polynomial_ring(QQ, ["x$i" for i in 1:220]) - GB = Groebner.groebner(xs) - @test GB == reverse(xs) - @test Groebner.normalform(GB, xs[1]) == R(0) - - R, xs = polynomial_ring(QQ, ["x$i" for i in 1:1025]) - GB = Groebner.groebner(xs) - @test GB == reverse(xs) - @test Groebner.normalform(GB, xs) == zeros(R, length(xs)) -end - -@testset "normalform orderings" begin - R, (x, y, z) = polynomial_ring(QQ, ["x", "y", "z"]) - - @test Groebner.normalform([x + y], [x, y]) == [-y, y] - @test Groebner.normalform([x + y^2], [x, y], ordering=Groebner.DegLex()) == [x, y] - @test Groebner.normalform([x + y^2], [x, y], ordering=Groebner.DegRevLex()) == [x, y] - @test Groebner.normalform([x + y], [x, y], ordering=Groebner.Lex(y, x, z)) == [x, -x] -end - -@testset "normalform of an array" begin - for field in [GF(17), GF(2^31 - 1), QQ] - R, (x, y) = polynomial_ring(field, ["x", "y"]) - gb = [x, y] - @test Groebner.normalform(gb, [x, y + 1]) == [R(0), R(1)] - @test Groebner.normalform(gb, [y + 1, x]) == [R(1), R(0)] - @test Groebner.normalform(gb, [x, 3y, 4x + 5y, 6x * y, 7x + 1, y + 2, x^2 + y^2 + 3]) == - [R(0), R(0), R(0), R(0), R(1), R(2), R(3)] - - gb = [x] - @test Groebner.normalform(gb, [x, 3y, 4x + 5y, 6x * y, 7x + 1, y + 2, x^2 + y^2 + 3]) == - [R(0), 3y, 5y, R(0), R(1), y + 2, y^2 + 3] - - gb = [x^2 + y^2, y^3 - y] - @test Groebner.normalform( - gb, - [R(1), R(5), R(16), 3x, y, 4x * y + 2, y^2, x * y^2 - 8x + y] - ) == [R(1), R(5), R(16), 3x, y, 4x * y + 2, y^2, x * y^2 - 8x + y] - - gb = [x, R(0), y, R(0)] - @test Groebner.normalform(gb, [R(0)]) == [R(0)] - @test Groebner.normalform(gb, [R(0), R(1), R(0)]) == [R(0), R(1), R(0)] - @test Groebner.normalform([R(0), R(0)], [R(0), R(1), R(0)]) == [R(0), R(1), R(0)] - end -end diff --git a/test/runtests.jl b/test/runtests.jl index cf63548d..163d645f 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -34,36 +34,27 @@ end @test isempty(Test.detect_ambiguities(Groebner)) @time @testset "All tests" verbose = true begin - # Basic tests for addition in Zp - @time @includetests ["arithmetic/Zp"] + @time @includetests ["arithmetic"] # Different implementations of a monomial @time @includetests ["monoms/exponentvector", "monoms/packedtuples"] - # High-level monomial arithmetic and term orders @time @includetests ["monoms/monom_arithmetic", "monoms/monom_orders"] - # Consistency of input-output - @time @includetests ["input_output/AbstractAlgebra"] - @time @includetests [ - "groebner/groebner", - "groebner/groebner_stress", - "groebner/homogenization", - "groebner/multi_threading" + "groebner" ] - @time @includetests ["groebner/groebner_with_change_matrix"] + @time @includetests ["learn_and_apply"] - @time @includetests ["learn_and_apply/learn_and_apply", "learn_and_apply/apply_in_batches"] + @time @includetests ["isgroebner"] - @time @includetests ["isgroebner/isgroebner"] - - @time @includetests ["normalform/normalform", "normalform/normalform_stress"] + @time @includetests ["normalform"] # Test for different frontends: - # - AbstractAlgebra.jl (AbstractAlgebra.Generic.MPoly{T}) - # - Nemo.jl (Nemo.fmpq_mpoly, Nemo.gfp_mpoly) - # - DynamicPolynomials.jl (DynamicPolynomials.Polynomial{true, T}) + # - AbstractAlgebra.jl + # - Nemo.jl + # - DynamicPolynomials.jl + @time @includetests ["input_output/AbstractAlgebra"] if try_import(:DynamicPolynomials) @time @includetests ["input_output/GroebnerDynamicPolynomialsExt"] end From 8f99a9e823616e4ed2148a98c259d17887b4269b Mon Sep 17 00:00:00 2001 From: Sasha Demin Date: Sat, 21 Sep 2024 00:05:39 +0200 Subject: [PATCH 5/6] save 1 prime --- src/groebner/groebner.jl | 2 ++ test/groebner.jl | 6 ++---- test/runtests.jl | 4 +--- 3 files changed, 5 insertions(+), 7 deletions(-) diff --git a/src/groebner/groebner.jl b/src/groebner/groebner.jl index 76b99bfb..54991c22 100644 --- a/src/groebner/groebner.jl +++ b/src/groebner/groebner.jl @@ -160,6 +160,8 @@ function _groebner_guess_lucky_prime( f4!(ring_ff_2, basis_ff_2, pairset, hashtable, params_zp) if basis_ff_1.monoms == basis_ff_2.monoms + push!(state.used_primes, prime_2) + push!(state.gb_coeffs_ff_all, basis_ff_2.coeffs) return prime_1 end diff --git a/test/groebner.jl b/test/groebner.jl index 75ff6e87..1e8f5e22 100644 --- a/test/groebner.jl +++ b/test/groebner.jl @@ -1117,10 +1117,8 @@ end # Also test learn / apply field == QQ && continue - context3, gb3 = - Groebner.groebner_learn(system, ordering=ord, homogenize=:yes) - context4, gb4 = - Groebner.groebner_learn(system, ordering=ord, homogenize=:no) + context3, gb3 = Groebner.groebner_learn(system, ordering=ord, homogenize=:yes) + context4, gb4 = Groebner.groebner_learn(system, ordering=ord, homogenize=:no) for _ in 1:4 flag3, gb33 = Groebner.groebner_apply!(context3, system) flag4, gb44 = Groebner.groebner_apply!(context4, system) diff --git a/test/runtests.jl b/test/runtests.jl index 163d645f..3727e93c 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -40,9 +40,7 @@ end @time @includetests ["monoms/exponentvector", "monoms/packedtuples"] @time @includetests ["monoms/monom_arithmetic", "monoms/monom_orders"] - @time @includetests [ - "groebner" - ] + @time @includetests ["groebner"] @time @includetests ["learn_and_apply"] From efb76ee9d7b1596bac4ca231d631b2cf49938870 Mon Sep 17 00:00:00 2001 From: Sasha Demin Date: Sat, 21 Sep 2024 00:19:20 +0200 Subject: [PATCH 6/6] revert 'save 1 prime' --- src/groebner/groebner.jl | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/groebner/groebner.jl b/src/groebner/groebner.jl index 54991c22..76b99bfb 100644 --- a/src/groebner/groebner.jl +++ b/src/groebner/groebner.jl @@ -160,8 +160,6 @@ function _groebner_guess_lucky_prime( f4!(ring_ff_2, basis_ff_2, pairset, hashtable, params_zp) if basis_ff_1.monoms == basis_ff_2.monoms - push!(state.used_primes, prime_2) - push!(state.gb_coeffs_ff_all, basis_ff_2.coeffs) return prime_1 end