From 0ec642c32a25fb101a0882da1450cdfdbf55b5ad Mon Sep 17 00:00:00 2001 From: Alexander Demin Date: Mon, 29 Jan 2024 21:06:15 +0300 Subject: [PATCH] update eco10 example --- src/f4/linalg/backend_threaded.jl | 145 +-------------- src/monomials/packedtuples.jl | 298 ++++++++++++++++++++++++++++-- src/monomials/packedutils.jl | 9 +- src/utils/examples.jl | 2 +- test/monoms/monom_arithmetic.jl | 14 ++ test/monoms/monom_orders.jl | 5 +- test/monoms/packedtuples.jl | 51 ++++- 7 files changed, 360 insertions(+), 164 deletions(-) diff --git a/src/f4/linalg/backend_threaded.jl b/src/f4/linalg/backend_threaded.jl index a12c71f3..45f73e0f 100644 --- a/src/f4/linalg/backend_threaded.jl +++ b/src/f4/linalg/backend_threaded.jl @@ -19,13 +19,7 @@ function linalg_deterministic_sparse_threaded!( @log level = -3 "linalg_deterministic_sparse!" @log level = -3 matrix_string_repr(matrix) # Reduce CD with AB - if true - linalg_reduce_matrix_lower_part_threaded_cas!(matrix, basis, arithmetic) - elseif false - linalg_reduce_matrix_lower_part_threaded_task_lock_free!(matrix, basis, arithmetic) - else - linalg_reduce_matrix_lower_part_threaded_task_cas!(matrix, basis, arithmetic) - end + linalg_reduce_matrix_lower_part_threaded_cas!(matrix, basis, arithmetic) # Interreduce CD linalg_interreduce_matrix_pivots!(matrix, basis, arithmetic) true @@ -42,8 +36,6 @@ end # A B # 0 D' # The new pivots in the D' block are known, but possibly not fully interreduced. -# -# NOTE: This function is incorrect. function linalg_reduce_matrix_lower_part_threaded_cas!( matrix::MacaulayMatrix{CoeffType}, basis::Basis{CoeffType}, @@ -62,137 +54,12 @@ function linalg_reduce_matrix_lower_part_threaded_cas!( resize!(matrix.some_coeffs, nlow) resize!(matrix.sentinels, ncols) - # If there is a pivot with the leading column i, then sentinels[i] = 1. - # Otherwise, sentinels[i] = 0. - # Once sentinels[i] becomes 1, this equality is unchanged. - # NOTE: for all valid indices, isassigned(sentinels, i) must hold. - sentinels = matrix.sentinels - @inbounds for i in 1:ncols - sentinels[i] = 0 - end - @inbounds for i in 1:nup - sentinels[matrix.upper_rows[i][1]] = 1 - end - - # Allocate thread-local buffers - buffers_row = map(_ -> zeros(AccumType, ncols), 1:nthreads()) - - # NOTE: by default, @threads uses the :dynamic execution schedule, which - # does not guarantee that threadid() is constant within one iteration - @inbounds Base.Threads.@threads for i in 1:nlow - t_id = threadid() - t_local_row = buffers_row[t_id] - new_sparse_row_support, new_sparse_row_coeffs = - linalg_new_empty_sparse_row(CoeffType) - - # Select a row to be reduced from the lower part of the matrix - # NOTE: no copy of coefficients is needed - sparse_row_support = matrix.lower_rows[i] - sparse_row_coeffs = basis.coeffs[row_index_to_coeffs[i]] - @invariant length(sparse_row_support) == length(sparse_row_coeffs) - - # Load the coefficients into a dense array - linalg_load_sparse_row!(t_local_row, sparse_row_support, sparse_row_coeffs) - - first_nnz_column = sparse_row_support[1] - - # In a somewhat ideal scenario, the while loop below executes just once - success = false - @inbounds while !success - @invariant 1 <= first_nnz_column <= length(t_local_row) - - # Note that this call may only mutate the first three arguments, and - # the arguments do not overlap in memory with each other - zeroed = linalg_reduce_dense_row_by_pivots_sparse!( - new_sparse_row_support, - new_sparse_row_coeffs, - t_local_row, - matrix, - basis, - pivots, - first_nnz_column, - ncols, - arithmetic, - tmp_pos=-1 - ) - - # If the row is fully reduced - zeroed && break - - # At this point, we have discovered a row that is not reduced to - # zero (yet). Thus, two outcomes are possible: - # 1. sentinels[...] is 0, meaning that this row is potentially a new - # pivot row. We try to set sentinels[...] to 1, update the pivots - # and the matrix, and break out of the loop. - # If assigning sentinels[...] to 1 fails, go to 2. - # 2. sentinels[...] is already 1, meaning that this row can be - # further reduced on the following iteration - - # Sync point. Everything before this point becomes visible to - # all other threads once they reach this point. - # NOTE: Note the absense of a total ordering on atomic operations - old, success = Atomix.replace!( - Atomix.IndexableRef(sentinels, (Int(new_sparse_row_support[1]),)), - Int8(0), - Int8(1), - Atomix.acquire_release, - Atomix.acquire - ) - - @invariant length(new_sparse_row_support) == length(new_sparse_row_coeffs) - - if success - # new pivot is formed, wake up from the inner loop - @invariant iszero(old) - - linalg_normalize_row!(new_sparse_row_coeffs, arithmetic) - - matrix.some_coeffs[i] = new_sparse_row_coeffs - matrix.lower_to_coeffs[new_sparse_row_support[1]] = i - pivots[new_sparse_row_support[1]] = new_sparse_row_support - else - # go for another iteration - first_nnz_column = new_sparse_row_support[1] - end - end - - new_sparse_row_support, new_sparse_row_coeffs = - linalg_new_empty_sparse_row(CoeffType) - end - - true -end - -# Given a matrix of the following form, where A is in REF, -# A B -# C D -# reduces the lower part CD with respect to the pivots in the upper part AB. -# As a result, a matrix of the following form is produced: -# A B -# 0 D' -# The new pivots in the D' block are known, but possibly not fully interreduced. -function linalg_reduce_matrix_lower_part_threaded_cas_maybe_correct!( - matrix::MacaulayMatrix{CoeffType}, - basis::Basis{CoeffType}, - arithmetic::AbstractArithmetic{AccumType, CoeffType} -) where {CoeffType <: Coeff, AccumType <: Coeff} - if nthreads() == 1 - @log level = 0 """ - Using multi-threaded linear algebra with nthreads() == 1. - Something probably went wrong.""" - end - _, ncols = size(matrix) - nup, nlow = matrix_nrows_filled(matrix) - - # Prepare the matrix - pivots, row_index_to_coeffs = linalg_prepare_matrix_pivots!(matrix) - resize!(matrix.some_coeffs, nlow) - resize!(matrix.sentinels, ncols) - - # If there is a pivot with the leading column i, then sentinels[i] = 1. + # If there is a pivot with the leading column i, then we set + # sentinels[i] = 1. + # If sentinels[i] = 1, and the pivot is synced, then we set + # sentinels[i] = 2. # Otherwise, sentinels[i] = 0. - # Once sentinels[i] becomes 1, this equality is unchanged. - # NOTE: for all valid indices, isassigned(sentinels, i) must hold. + # Once sentinels[i] becomes 1 or 2, this equality is unchanged. sentinels = matrix.sentinels @inbounds for i in 1:ncols sentinels[i] = 0 diff --git a/src/monomials/packedtuples.jl b/src/monomials/packedtuples.jl index 838c17ad..bbb74657 100644 --- a/src/monomials/packedtuples.jl +++ b/src/monomials/packedtuples.jl @@ -1,30 +1,31 @@ # This file is a part of Groebner.jl. License is GNU GPL v2. -# PackedTupleN - +### # PackedTupleN implements the monomial interface. -# PackedTupleN implements the interface of a vector of small nonegative integers -# with O(1) sum queries. # -# PackedTupleN always stores the sum of the vector explicitly. -# PackedTupleN always stores in degree-reverse-lex favorable order. +# PackedTupleN packs exponents in integers in a degrevlex-favorable ordering. +# Currently implemented packed monomials support up to 31 variables. abstract type AbstractPackedTuple{T <: Unsigned, B <: Unsigned} end struct PackedTuple1{T <: Unsigned, B <: Unsigned} <: AbstractPackedTuple{T, B} a1::T end - struct PackedTuple2{T <: Unsigned, B <: Unsigned} <: AbstractPackedTuple{T, B} a1::T a2::T end - struct PackedTuple3{T <: Unsigned, B <: Unsigned} <: AbstractPackedTuple{T, B} a1::T a2::T a3::T end +struct PackedTuple4{T <: Unsigned, B <: Unsigned} <: AbstractPackedTuple{T, B} + a1::T + a2::T + a3::T + a4::T +end # a `p` object can store monom_max_vars(p) integers at max. monom_max_vars(p::AbstractPackedTuple) = monom_max_vars(typeof(p)) @@ -35,7 +36,8 @@ function _monom_overflow_check(e::AbstractPackedTuple{T, B}) where {T, B} _monom_overflow_check(monom_totaldeg(e), B) end -const _defined_packed_pairs = ((:PackedTuple1, 1), (:PackedTuple2, 2), (:PackedTuple3, 3)) +const _defined_packed_pairs = + ((:PackedTuple1, 1), (:PackedTuple2, 2), (:PackedTuple3, 3), (:PackedTuple4, 4)) # for each PackedTupleI define something.. for (op, n) in _defined_packed_pairs @@ -72,13 +74,10 @@ for (op, n) in _defined_packed_pairs end end -monom_copy(pv::PackedTuple1{T, B}) where {T <: UInt64, B} = pv -monom_copy(pv::PackedTuple2{T, B}) where {T <: UInt64, B} = pv -monom_copy(pv::PackedTuple3{T, B}) where {T <: UInt64, B} = pv - -monom_copy(pv::PackedTuple1{T, B}) where {T, B} = PackedTuple1{T, B}(pv.a1) -monom_copy(pv::PackedTuple2{T, B}) where {T, B} = PackedTuple2{T, B}(pv.a1, pv.a2) -monom_copy(pv::PackedTuple3{T, B}) where {T, B} = PackedTuple3{T, B}(pv.a1, pv.a2, pv.a3) +monom_copy(pv::PackedTuple1{T, B}) where {T, B} = pv +monom_copy(pv::PackedTuple2{T, B}) where {T, B} = pv +monom_copy(pv::PackedTuple3{T, B}) where {T, B} = pv +monom_copy(pv::PackedTuple4{T, B}) where {T, B} = pv # Creates a packed monomial of the given type from regular vector `ev` function monom_construct_from_vector( @@ -165,6 +164,42 @@ function monom_construct_from_vector( a1 |= s << (indent * 8) PackedTuple3{T, B}(a1, a2, a3) end +function monom_construct_from_vector( + ::Type{PackedTuple4{T, B}}, + ev::Vector{U} +) where {T, B, U} + n = length(ev) + epc = packed_elperchunk(T, B) + @invariant n < 4 * epc + if n < 3 * epc + small = monom_construct_from_vector(PackedTuple3{T, B}, ev) + return PackedTuple4{T, B}(small.a1, small.a2, small.a3, zero(T)) + end + indent = sizeof(T) - packed_degsize(T, B, n) + a1, a2, a3, a4 = zero(T), zero(T), zero(T), zero(T) + s = zero(T) + @inbounds for i in n:-1:1 + _monom_overflow_check(ev[i], B) + d = T(ev[i]) + if div(i - 1, epc) == 3 + a1 = a1 << (sizeof(B) * 8) + a1 = a1 | d + elseif div(i - 1, epc) == 2 + a2 = a2 << (sizeof(B) * 8) + a2 = a2 | d + elseif div(i - 1, epc) == 1 + a3 = a3 << (sizeof(B) * 8) + a3 = a3 | d + else + a4 = a4 << (sizeof(B) * 8) + a4 = a4 | d + end + _monom_overflow_check(s, B) + s += d + end + a1 |= s << (indent * 8) + PackedTuple4{T, B}(a1, a2, a3, a4) +end # Creates a constant packed monomial of the given type of length n function monom_construct_const_monom(::Type{PackedTuple1{T, B}}, n::Integer) where {T, B} @@ -179,6 +214,10 @@ function monom_construct_const_monom(::Type{PackedTuple3{T, B}}, n::Integer) whe @invariant n < 3 * packed_elperchunk(T, B) PackedTuple3{T, B}(zero(T), zero(T), zero(T)) end +function monom_construct_const_monom(::Type{PackedTuple4{T, B}}, n::Integer) where {T, B} + @invariant n < 4 * packed_elperchunk(T, B) + PackedTuple4{T, B}(zero(T), zero(T), zero(T), zero(T)) +end # Hash of a packed monomial function monom_hash(x::PackedTuple1{T, B}, b::Vector{MH}) where {T, B, MH} @@ -210,6 +249,20 @@ function monom_hash(x::PackedTuple3{T, B}, b::Vector{MH}) where {T, B, MH} ) mod(h, MonomHash) end +function monom_hash(x::PackedTuple4{T, B}, b::Vector{MH}) where {T, B, MH} + epc = packed_elperchunk(T, B) + h = packed_dot_product(x.a4, b, B, 0) + h = packed_dot_product(x.a3, view(b, (epc + 1):(2 * epc)), B, 0) + h = h + packed_dot_product(x.a2, view(b, (2 * epc + 1):(3 * epc)), B, 0) + h = + h + packed_dot_product( + x.a1, + view(b, (3 * epc + 1):length(b)), + B, + epc - max(epc - 1, length(b) - 3 * epc) + ) + mod(h, MonomHash) +end # Creates a regular vector from a packed monomial and writes result to `tmp` function monom_to_vector!(tmp::Vector{I}, pv::PackedTuple1{T, B}) where {I, T, B} @@ -239,6 +292,20 @@ function monom_to_vector!(tmp::Vector{I}, pv::PackedTuple3{T, B}) where {I, T, B packed_unpack!(view(tmp, (2 * epc + 1):length(tmp)), pv.a1, B, indent) tmp end +function monom_to_vector!(tmp::Vector{I}, pv::PackedTuple4{T, B}) where {I, T, B} + epc = packed_elperchunk(T, B) + (length(tmp) < 3 * epc) && + return monom_to_vector!(tmp, PackedTuple3{T, B}(pv.a1, pv.a2, pv.a3)) + indent = 0 + packed_unpack!(tmp, pv.a4, B, indent) + indent = 0 + packed_unpack!(view(tmp, (epc + 1):(2 * epc)), pv.a3, B, indent) + indent = 0 + packed_unpack!(view(tmp, (2 * epc + 1):(3 * epc)), pv.a2, B, indent) + indent = epc - min(epc - 1, length(tmp) - 3 * epc) + packed_unpack!(view(tmp, (3 * epc + 1):length(tmp)), pv.a1, B, indent) + tmp +end ### # Monomial orderings for the `PackedTupleI` monomial implementation. @@ -315,6 +382,34 @@ function monom_isless( end end +function monom_isless( + ea::PackedTuple4{T, B}, + eb::PackedTuple4{T, B}, + ::_DegRevLex{true} +) where {T, B} + da, db = monom_totaldeg(ea), monom_totaldeg(eb) + if da < db + return true + end + if da > db + return false + end + + if ea.a1 == eb.a1 + if ea.a2 == eb.a2 + if ea.a3 == eb.a3 + return !(ea.a4 <= eb.a4) + else + return !(ea.a3 <= eb.a3) + end + else + return !(ea.a2 <= eb.a2) + end + else + return !(ea.a1 <= eb.a1) + end +end + function monom_isless( ea::PackedTuple1{T, B}, eb::PackedTuple1{T, B}, @@ -349,6 +444,17 @@ function monom_isless( b = monom_construct_from_vector(ExponentVector{T}, monom_to_vector!(tmp2, eb)) monom_isless(a, b, ord) end +function monom_isless( + ea::PackedTuple4{T, B}, + eb::PackedTuple4{T, B}, + ord::Ord +) where {T, B, Ord <: AbstractInternalOrdering} + s = 4 * div(sizeof(T), sizeof(B)) - 1 + tmp1, tmp2 = Vector{T}(undef, s), Vector{T}(undef, s) + a = monom_construct_from_vector(ExponentVector{T}, monom_to_vector!(tmp1, ea)) + b = monom_construct_from_vector(ExponentVector{T}, monom_to_vector!(tmp2, eb)) + monom_isless(a, b, ord) +end ### # Monomial-Monomial arithmetic. @@ -389,6 +495,20 @@ function monom_lcm!( _monom_overflow_check(ans) ans end +function monom_lcm!( + ec::PackedTuple4{T, B}, + ea::PackedTuple4{T, B}, + eb::PackedTuple4{T, B} +) where {T, B} + x1, si1 = packed_max(ea.a1, eb.a1, B, Val(1)) + x2, si2 = packed_max(ea.a2, eb.a2, B, Val(0)) + x3, si3 = packed_max(ea.a3, eb.a3, B, Val(0)) + x4, si4 = packed_max(ea.a4, eb.a4, B, Val(0)) + x1 = x1 + ((si1 + si2 + si3 + si4) << ((sizeof(T) - sizeof(B)) * 8)) + ans = PackedTuple4{T, B}(x1, x2, x3, x4) + _monom_overflow_check(ans) + ans +end function monom_is_gcd_const(ea::PackedTuple1{T, B}, eb::PackedTuple1{T, B}) where {T, B} if !packed_is_zero_dot_product(ea.a1, eb.a1, B, Val(1)) @@ -417,6 +537,21 @@ function monom_is_gcd_const(ea::PackedTuple3{T, B}, eb::PackedTuple3{T, B}) wher end return true end +function monom_is_gcd_const(ea::PackedTuple4{T, B}, eb::PackedTuple4{T, B}) where {T, B} + if !packed_is_zero_dot_product(ea.a1, eb.a1, B, Val(1)) + return false + end + if !packed_is_zero_dot_product(ea.a2, eb.a2, B, Val(0)) + return false + end + if !packed_is_zero_dot_product(ea.a3, eb.a3, B, Val(0)) + return false + end + if !packed_is_zero_dot_product(ea.a4, eb.a4, B, Val(0)) + return false + end + return true +end function monom_product!( ec::PackedTuple1{T, B}, @@ -451,6 +586,19 @@ function monom_product!( _monom_overflow_check(ans) ans end +function monom_product!( + ec::PackedTuple4{T, B}, + ea::PackedTuple4{T, B}, + eb::PackedTuple4{T, B} +) where {T, B} + x1 = ea.a1 + eb.a1 + x2 = ea.a2 + eb.a2 + x3 = ea.a3 + eb.a3 + x4 = ea.a4 + eb.a4 + ans = PackedTuple4{T, B}(x1, x2, x3, x4) + _monom_overflow_check(ans) + ans +end function monom_division!( ec::PackedTuple1{T, B}, @@ -482,6 +630,18 @@ function monom_division!( ans = PackedTuple3{T, B}(x1, x2, x3) ans end +function monom_division!( + ec::PackedTuple4{T, B}, + ea::PackedTuple4{T, B}, + eb::PackedTuple4{T, B} +) where {T, B} + x1 = ea.a1 - eb.a1 + x2 = ea.a2 - eb.a2 + x3 = ea.a3 - eb.a3 + x4 = ea.a4 - eb.a4 + ans = PackedTuple4{T, B}(x1, x2, x3, x4) + ans +end function monom_is_divisible(ea::PackedTuple1{T, B}, eb::PackedTuple1{T, B}) where {T, B} if !packed_ge(ea.a1, eb.a1, B, Val(1)) @@ -510,6 +670,21 @@ function monom_is_divisible(ea::PackedTuple3{T, B}, eb::PackedTuple3{T, B}) wher end return true end +function monom_is_divisible(ea::PackedTuple4{T, B}, eb::PackedTuple4{T, B}) where {T, B} + if !packed_ge(ea.a1, eb.a1, B, Val(1)) + return false + end + if !packed_ge(ea.a2, eb.a2, B, Val(0)) + return false + end + if !packed_ge(ea.a3, eb.a3, B, Val(0)) + return false + end + if !packed_ge(ea.a4, eb.a4, B, Val(0)) + return false + end + return true +end function monom_is_divisible!( ec::PackedTuple1{T, B}, @@ -541,6 +716,16 @@ function monom_is_divisible!( ans && (e = monom_division!(ec, ea, eb)) ans, e end +function monom_is_divisible!( + ec::PackedTuple4{T, B}, + ea::PackedTuple4{T, B}, + eb::PackedTuple4{T, B} +) where {T, B} + ans = monom_is_divisible(ea, eb) + e = ec + ans && (e = monom_division!(ec, ea, eb)) + ans, e +end function monom_is_equal(ea::PackedTuple1{T, B}, eb::PackedTuple1{T, B}) where {T, B} ea.a1 == eb.a1 @@ -551,6 +736,9 @@ end function monom_is_equal(ea::PackedTuple3{T, B}, eb::PackedTuple3{T, B}) where {T, B} ea.a1 == eb.a1 && ea.a2 == eb.a2 && ea.a3 == eb.a3 end +function monom_is_equal(ea::PackedTuple4{T, B}, eb::PackedTuple4{T, B}) where {T, B} + ea.a1 == eb.a1 && ea.a2 == eb.a2 && ea.a3 == eb.a3 && ea.a4 == eb.a4 +end ### # Monomial division masks. @@ -701,3 +889,81 @@ function monom_create_divmask( res end + +function monom_create_divmask( + e::PackedTuple4{T, B}, + DM::Type{Mask}, + ndivvars, + divmap, + ndivbits, + compressed +) where {T, B, Mask} + @invariant !compressed + + epc = packed_elperchunk(T, B) + + if ndivvars < 3 * epc + return monom_create_divmask( + PackedTuple3{T, B}(e.a1, e.a2, e.a3), + DM, + ndivvars, + divmap, + ndivbits, + compressed + ) + end + + ctr = one(Mask) + res = zero(Mask) + o = one(Mask) + + a4 = e.a4 + for i in 1:epc + ei = mod(a3, B) + a3 = a3 >> (sizeof(B) * 8) + for j in 1:ndivbits + @inbounds if ei >= divmap[ctr] + res |= o << (ctr - 1) + end + ctr += o + end + end + + a3 = e.a3 + for i in (epc + 1):(2 * epc) + ei = mod(a3, B) + a3 = a3 >> (sizeof(B) * 8) + for j in 1:ndivbits + @inbounds if ei >= divmap[ctr] + res |= o << (ctr - 1) + end + ctr += o + end + end + + a2 = e.a2 + for i in (2 * epc + 1):(3 * epc) + ei = mod(a2, B) + a2 = a2 >> (sizeof(B) * 8) + for j in 1:ndivbits + @inbounds if ei >= divmap[ctr] + res |= o << (ctr - 1) + end + ctr += o + end + end + + a1 = e.a1 + for i in (3 * epc + 1):min(4 * epc - 1, ndivvars) + ei = mod(a1, B) + a1 = a1 >> (sizeof(B) * 8) + for j in 1:ndivbits + @inbounds if ei >= divmap[ctr] + res |= o << (ctr - 1) + end + ctr += o + end + end + + res +end diff --git a/src/monomials/packedutils.jl b/src/monomials/packedutils.jl index 7657493a..28c6d243 100644 --- a/src/monomials/packedutils.jl +++ b/src/monomials/packedutils.jl @@ -3,10 +3,8 @@ ### # Utilities for packed integers -# Functions here work on integers packed into a single integer of a wider type. - # How many integers of type B can be stored in an integer of type T -function packed_elperchunk(T, B) +function packed_elperchunk(::Type{T}, ::Type{B}) where {T, B} epc = div(sizeof(T), sizeof(B)) @invariant epc * sizeof(B) == sizeof(T) epc @@ -14,8 +12,9 @@ end # the size of total degree. # By default, does not differ from the size of other elements -packed_degsize(T, B, n) = sizeof(B) -packed_nchunks(T, B, n) = div((n - 1) * sizeof(B) + packed_degsize(T, B, n), sizeof(T)) + 1 +packed_degsize(::Type{T}, ::Type{B}, n) where {T, B} = sizeof(B) +packed_nchunks(::Type{T}, ::Type{B}, n) where {T, B} = + div((n - 1) * sizeof(B) + packed_degsize(T, B, n), sizeof(T)) + 1 @inline @generated function packed_unpack!( b::AbstractVector{MH}, diff --git a/src/utils/examples.jl b/src/utils/examples.jl index 2b7582d0..8d64f43c 100644 --- a/src/utils/examples.jl +++ b/src/utils/examples.jl @@ -186,7 +186,7 @@ end function eco10(; np=AbstractAlgebra, k=np.QQ, ordering=:lex) _, (x0, x1, x2, x3, x4, x5, x6, x7, x8, x9) = - np.polynomial_ring(k, ["x$i" for i in 1:11], ordering=ordering) + np.polynomial_ring(k, ["x$i" for i in 1:10], ordering=ordering) [ x0 * x1 * x9 + x1 * x2 * x9 + diff --git a/test/monoms/monom_arithmetic.jl b/test/monoms/monom_arithmetic.jl index ff9c5ab8..9f83dc21 100644 --- a/test/monoms/monom_arithmetic.jl +++ b/test/monoms/monom_arithmetic.jl @@ -17,6 +17,7 @@ implementations_to_test = [ Groebner.PackedTuple1{T, UInt8} where {T}, Groebner.PackedTuple2{T, UInt8} where {T}, Groebner.PackedTuple3{T, UInt8} where {T}, + Groebner.PackedTuple4{T, UInt8} where {T}, Groebner.SparseExponentVector{T} where {T} ] @@ -97,6 +98,19 @@ implementations_to_test = [ @test monom_is_equal(tmp, e) @test !monom_is_gcd_const(a, b) @test monom_is_equal(tmp, e) + + n = 30 + if n > Groebner.monom_max_vars(MonomType) + continue + end + a = monom_construct_from_vector(MonomType, 2 .* ones(Int, n)) + b = monom_construct_from_vector(MonomType, 1 .* ones(Int, n)) + tmp = monom_copy(a) + @test monom_is_divisible(a, b) + flag, tmp = monom_is_divisible!(tmp, a, b) + @test flag + @test monom_is_equal(tmp, b) + @test !monom_is_gcd_const(a, b) end # test that different implementations agree diff --git a/test/monoms/monom_orders.jl b/test/monoms/monom_orders.jl index 86ff1169..f79ea4b2 100644 --- a/test/monoms/monom_orders.jl +++ b/test/monoms/monom_orders.jl @@ -9,6 +9,7 @@ implementations_to_test = [ Groebner.PackedTuple1{T, UInt8} where {T}, Groebner.PackedTuple2{T, UInt8} where {T}, Groebner.PackedTuple3{T, UInt8} where {T}, + Groebner.PackedTuple4{T, UInt8} where {T}, Groebner.SparseExponentVector{T} where {T} ] @@ -186,7 +187,9 @@ end @testset "monoms, variable permutation" begin for T in (UInt64, UInt32, UInt16) for EV in implementations_to_test - if EV{T} <: Groebner.PackedTuple2 || EV{T} <: Groebner.PackedTuple3 + if EV{T} <: Groebner.PackedTuple2 || + EV{T} <: Groebner.PackedTuple3 || + EV{T} <: Groebner.PackedTuple4 continue end diff --git a/test/monoms/packedtuples.jl b/test/monoms/packedtuples.jl index d6c6b508..0c1bf3fe 100644 --- a/test/monoms/packedtuples.jl +++ b/test/monoms/packedtuples.jl @@ -111,8 +111,8 @@ end tmp = similar(x) @test Groebner.monom_to_vector!(tmp, ev) == [1, 2, 3, 0, 4] - @test Groebner.monom_max_vars(Groebner.PackedTuple2{UInt64, UInt8}) == 15 - @test Groebner.monom_max_vars(Groebner.PackedTuple2{UInt64, UInt16}) == 7 + @test Groebner.monom_max_vars(Groebner.PackedTuple3{UInt64, UInt8}) == 23 + @test Groebner.monom_max_vars(Groebner.PackedTuple3{UInt64, UInt16}) == 11 ev = Groebner.monom_construct_from_vector(EV{UInt64, UInt16}, x) @test (ev.a1, ev.a2) == (0x000a000000000004, 0x0000000300020001) @@ -146,3 +146,50 @@ end end end end + +@testset "packed exponent tuple-4" begin + EV{T, U} = Groebner.PackedTuple4{T, U} where {T, U} + + x = [1, 2, 3, 0, 4] + ev = Groebner.monom_construct_from_vector(EV{UInt64, UInt8}, x) + @test ev.a1 == 0x0a00000400030201 + @test ev.a2 == 0x0000000000000000 + @test ev.a3 == 0x0000000000000000 + @test ev.a4 == 0x0000000000000000 + @test typeof(Groebner.monom_totaldeg(ev)) === UInt64 + @test Groebner.monom_totaldeg(ev) == UInt64(10) + @test Groebner.monom_entrytype(ev) === Groebner.MonomHash + tmp = similar(x) + @test Groebner.monom_to_vector!(tmp, ev) == [1, 2, 3, 0, 4] + + @test Groebner.monom_max_vars(Groebner.PackedTuple4{UInt64, UInt8}) == 31 + @test typeof(Groebner.monom_totaldeg(ev)) === UInt64 + @test Groebner.monom_totaldeg(ev) == UInt16(10) + + y = [10, 20, 30, 40, 50] + @test_throws Groebner.MonomialDegreeOverflow Groebner.monom_construct_from_vector( + EV{UInt64, UInt8}, + y + ) + + y = [1, 2, 3, 4, 5] + ev = Groebner.monom_construct_from_vector(EV{UInt64, UInt32}, y) + @test ev.a1 == 0x0000000f00000005 + @test ev.a2 == 0x0000000400000003 + @test ev.a3 == 0x0000000200000001 + + for B in (UInt8, UInt16) + for T in (UInt64, UInt32) + for k in 1:10 + maxvars = 4 * div(sizeof(T), sizeof(B)) - 1 + n = rand(1:maxvars) + @assert maxvars > 0 + x = rand(0:div(typemax(B), 4 * n), n) + tmp = similar(x) + ev = Groebner.monom_construct_from_vector(EV{T, B}, x) + @test x == Groebner.monom_to_vector!(tmp, ev) + @test x == tmp + end + end + end +end