Skip to content

Commit

Permalink
Merge pull request #151 from sumiya11/cache
Browse files Browse the repository at this point in the history 
Faster monomial arithmetic
  • Loading branch information
@sumiya11
sumiya11 authored Sep 1, 2024
2 parents d8a14fa + 1f3a025 commit 7d2c65f
Show file tree
Hide file tree
Showing 9 changed files with 285 additions and 431 deletions.
44 changes: 44 additions & 0 deletions benchmark/CI-scripts/run_benchmarks.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -74,6 +74,7 @@ function compute_gb(system, trials=7; kws...)
end

# Compute Groebner bases over integers modulo a large prime
compute_gb(Groebner.Examples.katsuran(9, internal_ordering=:degrevlex, k=GF(2^31 - 1)))
push!(
suite,
(
Expand Down Expand Up @@ -139,6 +140,48 @@ push!(
)
)
)
push!(
suite,
(
problem_name="groebner,AA,2^30+3,1000xcyclic 6",
type=:time,
result=[
sum(
compute_gb(
Groebner.Examples.cyclicn(
6,
internal_ordering=:degrevlex,
k=GF(2^30 + 3)
),
2000
)
)
]
)
)
push!(
suite,
(
problem_name="groebner,AA,QQ,100xchandra 7",
type=:time,
result=[
sum(
compute_gb(
Groebner.Examples.chandran(7, internal_ordering=:degrevlex, k=QQ),
100
)
)
]
)
)
push!(
suite,
(
problem_name="groebner,AA,QQ,HIV2",
type=:time,
result=compute_gb(Groebner.Examples.HIV2(internal_ordering=:degrevlex), 3)
)
)
push!(
suite,
(
Expand Down Expand Up @@ -172,6 +215,7 @@ function n_variable_set(n; internal_ordering=:degrevlex, k=GF(2^31 - 1))
f
end

compute_gb(n_variable_set(100, internal_ordering=:degrevlex, k=GF(2^31 - 1)))
push!(
suite,
(
Expand Down
2 changes: 1 addition & 1 deletion src/Groebner.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -96,6 +96,7 @@ end
include("utils/logging.jl")
include("utils/invariants.jl")
include("utils/simd.jl")
include("utils/packed.jl")
include("utils/plots.jl")

# Test systems, such as katsura, cyclic, etc
Expand All @@ -108,7 +109,6 @@ include("utils/keywords.jl")

# Monomial implementations
include("monomials/exponent_vector.jl")
include("monomials/packed_utils.jl")
include("monomials/packed_vector.jl")

include("arithmetic/CompositeNumber.jl")
Expand Down
40 changes: 1 addition & 39 deletions src/f4/sort.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -27,12 +27,10 @@ function sort_part!(
nothing
end

# Also sorts any arrays passed in the `abc` optional argument in the same order.
function sort_polys_by_lead_increasing!(
basis::Basis,
hashtable::MonomialHashtable,
changematrix::Bool,
abc...;
changematrix::Bool;
ord::Ord=hashtable.ord
) where {Ord <: AbstractMonomialOrdering}
b_monoms = basis.monoms
Expand All @@ -52,10 +50,6 @@ function sort_polys_by_lead_increasing!(
# (seems to compile to better code)
basis.monoms[1:(basis.n_filled)] = basis.monoms[permutation]
basis.coeffs[1:(basis.n_filled)] = basis.coeffs[permutation]
@inbounds for a in abc
@invariant length(a) >= length(permutation)
a[1:(basis.n_filled)] = a[permutation]
end

if changematrix
@invariant length(basis.changematrix) >= basis.n_filled
Expand Down Expand Up @@ -252,35 +246,3 @@ function sort_input_terms_to_change_ordering!(
end
permutations
end

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

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

sort_part!(monoms, 1, cnt, lt=cmps, alg=_default_sorting_alg())
end

function sort_term_indices_decreasing!(
monoms::Vector{MonomId},
coeffs::Vector{C},
hashtable::MonomialHashtable,
ord::AbstractMonomialOrdering
) where {C <: Coeff}
exps = hashtable.monoms

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

inds = collect(1:length(monoms))

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

monoms[1:end] = monoms[inds]
coeffs[1:end] = coeffs[inds]
end
161 changes: 3 additions & 158 deletions src/input_output/AbstractAlgebra.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -109,7 +109,9 @@ function io_extract_coeffs_ir_qq(ring::PolyRing, polys)
end

function io_extract_monoms_ir(ring::PolyRing, polys)
var_to_index = get_var_to_index(AbstractAlgebra.parent(polys[1]))
ring_aa = AbstractAlgebra.parent(polys[1])
v = AbstractAlgebra.gens(ring_aa)
var_to_index = Dict{elem_type(ring_aa), Int}(v .=> 1:AbstractAlgebra.nvars(ring_aa))
res = Vector{Vector{Vector{UInt64}}}(undef, length(polys))
@inbounds for i in 1:length(polys)
poly = polys[i]
Expand Down Expand Up @@ -148,14 +150,6 @@ function _io_check_input(polynomials::Vector{T}) where {T}
true
end

# Determines the monomial ordering of the output,
# given the original ordering `origord` and the targer ordering `targetord`
ordering_typed2sym(origord, targetord::Lex) = :lex
ordering_typed2sym(origord, targetord::DegLex) = :deglex
ordering_typed2sym(origord, targetord::DegRevLex) = :degrevlex
ordering_typed2sym(origord) = origord
ordering_typed2sym(origord, targetord::AbstractMonomialOrdering) = origord

function ordering_sym2typed(ord::Symbol)
if !(ord in aa_supported_orderings)
__throw_input_not_supported(ord, "Not a supported ordering.")
Expand Down Expand Up @@ -218,158 +212,9 @@ function io_extract_coeffs_raw_X!(trace, coeffs)
true
end

###

# specialization for multivariate polynomials
function io_extract_coeffs_qq(representation, ring::PolyRing, poly)
iszero(poly) && (return Vector{representation.coefftype}())
n = length(poly)
arr = Vector{Rational{BigInt}}(undef, n)
@inbounds for i in 1:n
arr[i] = Rational{BigInt}(AbstractAlgebra.coeff(poly, i))
end
arr
end

function get_var_to_index(
aa_ring::Union{AbstractAlgebra.MPolyRing{T}, AbstractAlgebra.PolyRing{T}}
) where {T}
v = AbstractAlgebra.gens(aa_ring)
Dict{elem_type(aa_ring), Int}(v .=> 1:AbstractAlgebra.nvars(aa_ring))
end

function io_extract_monoms(
representation::PolynomialRepresentation,
ring::PolyRing,
poly::T
) where {T}
exps = Vector{representation.monomtype}(undef, length(poly))
_io_extract_monoms!(representation.monomtype, exps, poly)
exps
end

function _io_extract_monoms!(::Type{MonomType}, exps, poly) where {MonomType}
@inbounds for j in 1:length(exps)
exps[j] =
monom_construct_from_vector(MonomType, AbstractAlgebra.exponent_vector(poly, j))
end
nothing
end

function io_extract_monoms(
representation::PolynomialRepresentation,
ring::PolyRing,
poly::P
) where {P <: AbstractAlgebra.Generic.PolyRingElem}
exps = Vector{representation.monomtype}(undef, 0)
@inbounds while !AbstractAlgebra.iszero(poly)
push!(
exps,
monom_construct_from_vector(
representation.monomtype,
[AbstractAlgebra.degree(poly)]
)
)
poly = AbstractAlgebra.tail(poly)
end
exps
end

function io_extract_monoms(
representation::PolynomialRepresentation,
ring::PolyRing,
orig_polys::Vector{T}
) where {T}
npolys = length(orig_polys)
var_to_index = get_var_to_index(AbstractAlgebra.parent(orig_polys[1]))
exps = Vector{Vector{representation.monomtype}}(undef, npolys)
@inbounds for i in 1:npolys
poly = orig_polys[i]
exps[i] = io_extract_monoms(representation, ring, poly)
end
false, var_to_index, exps
end

function io_extract_monoms(
representation::PolynomialRepresentation,
ring::PolyRing,
orig_polys::Vector{T},
::DegLex
) where {T}
npolys = length(orig_polys)
var_to_index = get_var_to_index(AbstractAlgebra.parent(orig_polys[1]))
exps = Vector{Vector{representation.monomtype}}(undef, npolys)
@inbounds for i in 1:npolys
poly = orig_polys[i]
exps[i] = Vector{representation.monomtype}(undef, length(poly))
for j in 1:length(poly)
exps[i][j] = monom_construct_from_vector(
representation.monomtype,
poly.exps[(end - 1):-1:1, j]
)
end
end
false, var_to_index, exps
end

function io_extract_monoms(
representation::PolynomialRepresentation,
ring::PolyRing,
orig_polys::Vector{T},
::Lex
) where {T}
npolys = length(orig_polys)
var_to_index = get_var_to_index(AbstractAlgebra.parent(orig_polys[1]))
exps = Vector{Vector{representation.monomtype}}(undef, npolys)
@inbounds for i in 1:npolys
poly = orig_polys[i]
exps[i] = Vector{representation.monomtype}(undef, length(poly))
for j in 1:length(poly)
exps[i][j] = monom_construct_from_vector(
representation.monomtype,
poly.exps[end:-1:1, j]
)
end
end
false, var_to_index, exps
end

function io_extract_monoms(
representation::PolynomialRepresentation,
ring::PolyRing,
orig_polys::Vector{T},
::DegRevLex
) where {T}
npolys = length(orig_polys)
var_to_index = get_var_to_index(AbstractAlgebra.parent(orig_polys[1]))
exps = Vector{Vector{representation.monomtype}}(undef, npolys)
@inbounds for i in 1:npolys
poly = orig_polys[i]
exps[i] = Vector{representation.monomtype}(undef, length(poly))
for j in 1:length(poly)
exps[i][j] = monom_construct_from_vector(
representation.monomtype,
poly.exps[1:(end - 1), j]
)
end
end
false, var_to_index, exps
end

###
# Converting from internal representation to AbstractAlgebra.jl

function _io_convert_ir_to_polynomials(
ring::PolyRing,
polynomials,
monoms::Vector{Vector{M}},
coeffs::Vector{Vector{C}},
params
) where {M <: Monom, C <: Coeff}
origring = AbstractAlgebra.parent(first(polynomials))
_io_convert_ir_to_polynomials(origring, monoms, coeffs, params)
end

# Specialization for univariate polynomials
function _io_convert_ir_to_polynomials(
origring::R,
Expand Down
Loading

0 comments on commit 7d2c65f

Please sign in to comment.