diff --git a/src/Modules/Modules.jl b/src/Modules/Modules.jl index 36d2eef8f135..2250bd9944a2 100644 --- a/src/Modules/Modules.jl +++ b/src/Modules/Modules.jl @@ -1,4 +1,5 @@ include("ModuleTypes.jl") +include("hilbert.jl") include("UngradedModules.jl") include("homological-algebra.jl") include("FreeModElem-orderings.jl") diff --git a/src/Modules/ModulesGraded.jl b/src/Modules/ModulesGraded.jl index 224513ae4906..313f0116bed2 100644 --- a/src/Modules/ModulesGraded.jl +++ b/src/Modules/ModulesGraded.jl @@ -572,7 +572,7 @@ function graded_map(F::FreeMod{T}, A::MatrixElem{T}) where {T <: RingElement} return phi end -function graded_map(F::FreeMod{T}, V::Vector{<:FreeModElem{T}}) where {T <: RingElement} +function graded_map(F::FreeMod{T}, V::Vector{<:AbstractFreeModElem{T}}) where {T <: RingElement} R = base_ring(F) G = grading_group(R) nrows = length(V) diff --git a/src/Modules/hilbert.jl b/src/Modules/hilbert.jl new file mode 100644 index 000000000000..155e7279c06c --- /dev/null +++ b/src/Modules/hilbert.jl @@ -0,0 +1,187 @@ +# This cannot be included in src/Rings/hilbert.jl because the include +# order defined in src/Oscar.jl includes ring stuff before module stuff +# (and it probably would not work if we inverted the inclusion order). + +# ================================================================== +# Hilbert series numerator for modules + +function _power_product(T, expv) + return prod([T[k]^expv[k] for k in 1:length(expv)]); +end + + +function HSNum_module(SubM::SubquoModule{T}, HSRing::Ring, backend::Symbol=:Abbott) where T <: MPolyRingElem + C,phi = present_as_cokernel(SubM, :with_morphism); # phi is the morphism + LM = leading_module(C.quo); + F = ambient_free_module(C.quo); + r = rank(F); + P = base_ring(C.quo); + # short-cut for module R^0 (to avoid problems with empty sum below) + if iszero(r) + return multi_hilbert_series(quo(P,ideal(P,[1]))[1]; parent=HSRing, backend=backend)[1][1] + end + GensLM = gens(LM); + L = [[] for _ in 1:r]; # L[k] will be list of monomial gens for k-th cooord + # Nested loop below extracts the coordinate monomial ideals -- is there a better way? + for g in GensLM + SR = coordinates(ambient_representative(g)); # should have length = 1 + for j in 1:r + if SR[j] != 0 + push!(L[j], SR[j]); + end + end + end + IdealList = [ideal(P,G) for G in L]; + HSeriesList = [multi_hilbert_series(quo(P,I)[1]; parent=HSRing, backend=backend)[1][1] for I in IdealList]; + shifts = [degree(phi(g)) for g in gens(F)]; + @vprintln :hilbert 1 "HSNum_module: shifts are $(shifts)"; + shift_expv = [gen_repr(d) for d in shifts]; + @vprintln :hilbert 1 "HSNum_module: shift_expv are $(shift_expv)"; +### HSeriesRing = parent(HSeriesList[1]); +### @vprintln :hilbert 1 "HSNum_module: HSeriesRing = $(HSeriesRing)"; + t = gens(HSRing); + ScaleFactor = [_power_product(t,e) for e in shift_expv]; + result = sum([ScaleFactor[k]*HSeriesList[k] for k in 1:r]); + return result; +end + +# No longer needed +# function HSNum_module(F::FreeMod{T}; parent::Union{Nothing,Ring} = nothing) where T <: MPolyRingElem +# # ASSUME F is graded free module +# return HSNum_module(sub(F,gens(F))[1]; parent=parent) +# end + +@doc raw""" + multi_hilbert_series(M::SubquoModule; parent::Union{Nothing,Ring} = nothing) + +Compute a pair of pairs `(N ,D), (H ,iso)` where `N` and `D` are the non-reduced numerator and denominator of the Hilbert +series of the subquotient `M`, and `H` is the SNF of the grading group together with the identifying isomorphism `iso`. +If the kwarg `parent` is supplied `N` and `D` are computed in the ring `parent`. + +!!! note + CURRENT LIMITATION: the grading must be over ZZ^m (more general gradings are not yet supported) + +!!! note + Applied to a homogeneous subquotient `M`, the function first computes a Groebner basis to + obtain the leading term module; the rest of the computation uses this latter module + (sliced into ideals, one for each ambient free module component). + +# Examples +```jldoctest +julia> R, _ = polynomial_ring(QQ, ["x", "y", "z"]); + +julia> Rg, (x, y, z) = grade(R, [4,3,2]); + +julia> F = graded_free_module(Rg, 1); + +julia> A = Rg[x; y]; + +julia> B = Rg[x^2; y^3; z^4]; + +julia> M = SubquoModule(F, A, B); + +julia> (num,den),_ = multi_hilbert_series(M); + +julia> num +-t^25 + 2*t^17 + t^16 + t^15 - t^12 - t^11 - t^9 - t^8 - t^7 + t^4 + t^3 + +julia> den +Factored element with data +Dict{AbstractAlgebra.Generic.LaurentMPolyWrap{ZZRingElem, ZZMPolyRingElem, AbstractAlgebra.Generic.LaurentMPolyWrapRing{ZZRingElem, ZZMPolyRing}}, ZZRingElem}(-t^4 + 1 => 1, -t^3 + 1 => 1, -t^2 + 1 => 1) + +``` +""" +function multi_hilbert_series(SubM::SubquoModule{T}; parent::Union{Nothing, Ring} = nothing, backend::Symbol = :Abbott) where T <: MPolyRingElem + R = base_ring(SubM) + @req coefficient_ring(R) isa AbstractAlgebra.Field "The coefficient ring must be a field" + + # Wrap the case where G is abstractly isomorphic to ℤᵐ, but not realized as a + # free Abelian group. + # + # We use the Smith normal form to get there, recreate the graded ring with the + # free grading group, do the computation there and return the isomorphism for + # the grading. + G = grading_group(R) + if !is_zm_graded(R) + error("non-ZZ^m-grading not yet implemented for modules (see issue #2657)") ### Needs improvement to change_base_ring (see below) + H, iso = snf(G) + V = [preimage(iso, x) for x in gens(G)] + isoinv = hom(G, H, V) + W = [isoinv(R.d[i]) for i = 1:ngens(R)] + S, _ = graded_polynomial_ring(coefficient_ring(R), symbols(R), W) + map_into_S = hom(R, S, gens(S)) + SubM2,_ = change_base_ring(map_into_S,SubM) # !!! BUG this seems to forget that things are graded BUG (issue #2657) !!! + (numer, denom), _ = hilbert_series(SubM2; parent=parent, backend=backend) + return (numer, denom), (H, iso) + end + + # Now we may assume that the grading group is free Abelian. + m = ngens(G) + n = ngens(R) + HSRing = _hilbert_series_ring(parent, m) + # Get the weights as Int values: W[k] contains the weight(s) of x[k] + W = [[ Int(R.d[i][j]) for j in 1:m] for i in 1:n] + denom = _hilbert_series_denominator(HSRing, W) + numer = HSNum_module(SubM, HSRing, backend) + return (numer, denom), (G, identity_map(G)) +end + +function multi_hilbert_series(F::FreeMod{T}; parent::Union{Nothing,Ring} = nothing, backend::Symbol = :Abbott) where T <: MPolyRingElem + return multi_hilbert_series(sub(F,gens(F))[1]; parent=parent, backend=backend) +end + + +@doc raw""" + hilbert_series(M::SubquoModule; parent::Union{Nothing,Ring} = nothing) + +Compute a pair `(N,D)` where `N` and `D` are the non-reduced numerator and denominator of the Hilbert +series of the subquotient `M`. If the kwarg `parent` is supplied `N` and `D` are computed in the ring `parent`. + +!!! note + Applied to a homogeneous subquotient `M`, the function first computes a Groebner basis to + obtain the leading term module; the rest of the computation uses this latter module + (sliced into ideals, one for each ambient free module component). + +# Examples +```jldoctest +julia> R, _ = polynomial_ring(QQ, ["x", "y", "z"]); + +julia> Z = abelian_group(0); + +julia> Rg, (x, y, z) = grade(R, [Z[1],Z[1],Z[1]]); + +julia> F = graded_free_module(Rg, 1); + +julia> A = Rg[x; y]; + +julia> B = Rg[x^2; y^3; z^4]; + +julia> M = SubquoModule(F, A, B); + +julia> num,den = hilbert_series(M); + +julia> num +-t^9 + t^7 + 2*t^6 - t^5 - t^3 - 2*t^2 + 2*t + +julia> den +Factored element with data +Dict{AbstractAlgebra.Generic.LaurentPolyWrap{ZZRingElem, ZZPolyRingElem, AbstractAlgebra.Generic.LaurentPolyWrapRing{ZZRingElem, ZZPolyRing}}, ZZRingElem}(-t + 1 => 3) + +``` +""" +function hilbert_series(SubM::SubquoModule{T}; parent::Union{Nothing,Ring} = nothing, backend::Symbol = :Abbott) where T <: MPolyRingElem + @req is_z_graded(base_ring(SubM)) "ring must be ZZ-graded; use `multi_hilbert_series` otherwise" + if parent === nothing + parent, _ = LaurentPolynomialRing(ZZ, :t) + end + HS, _ = multi_hilbert_series(SubM; parent=parent, backend=backend) + return HS +end + +function hilbert_series(F::FreeMod{T}; parent::Union{Nothing,Ring} = nothing, backend::Symbol = :Abbott) where T <: MPolyRingElem + @req is_z_graded(base_ring(F)) "ring must be ZZ-graded; use `multi_hilbert_series` otherwise" + if parent === nothing + parent, _ = LaurentPolynomialRing(ZZ, :t) + end + return hilbert_series(sub(F,gens(F))[1]; parent=parent, backend=backend) +end diff --git a/src/Oscar.jl b/src/Oscar.jl index f7cb2aeb1043..fcbc5b49b24b 100644 --- a/src/Oscar.jl +++ b/src/Oscar.jl @@ -115,6 +115,8 @@ function __init__() add_verbose_scope(:Blowup) add_assert_scope(:Blowup) + add_verbose_scope(:hilbert) + add_assert_scope(:hilbert) add_verbose_scope(:GlobalTateModel) add_verbose_scope(:WeierstrassModel) diff --git a/src/Rings/Rings.jl b/src/Rings/Rings.jl index 441288d21936..021719c2e2c6 100644 --- a/src/Rings/Rings.jl +++ b/src/Rings/Rings.jl @@ -3,6 +3,7 @@ include("orderings.jl") include("mpoly.jl") include("mpoly_types.jl") include("mpoly-graded.jl") +include("hilbert_zach.jl") include("mpoly-ideals.jl") include("groebner.jl") include("solving.jl") @@ -33,4 +34,5 @@ include("AlgClosureFp.jl") include("PBWAlgebra.jl") include("PBWAlgebraQuo.jl") include("FreeAssAlgIdeal.jl") +include("hilbert.jl") diff --git a/src/Rings/hilbert.jl b/src/Rings/hilbert.jl new file mode 100644 index 000000000000..f4bd52e4acca --- /dev/null +++ b/src/Rings/hilbert.jl @@ -0,0 +1,1088 @@ +################################################################## +# Auxiliary functions +# Perhaps document them and make them publicly available? + +# random_subset: +# A function similar to this is in StatsBase.jl but I wish to avoid +# having an external dependency (for just 1 simple fn!) + +# Generate random m-subset of {1,2,3,...,n} +# Result is a Int64[]; entries are NOT SORTED!! +function random_subset(n::Int64, m::Int64) + # assume n >= 1, m >= 0 and m <= n + if m == 0 + return Int64[] + end + L = collect(1:n) + if m == n + return L + end + for j in 1:m + k = rand(j:n) + L[j],L[k] = L[k],L[j] # just a SWAP + end + L = first(L,m) + #?? sort!(L) ?? + return L +end + +# There must be a better way...! +# Split a "list" into 2 parts determined by a predicate. +# Returns 2-tuple: list-of-sat-elems, list-of-unsat-elems +function filter2(pred::Function, L::Vector) + sat = [] + unsat = [] + for x in L + if pred(x) + push!(sat,x) + else + push!(unsat,x) + end + end + return sat,unsat +end + + +################################################################## +# Code for representing & manipulating PPs (power products, aka. monomials) +# All this code is "local" to this file, & not exported! + + +# Each PP is represented as Vector{PP_exponent} + +PP_exponent = Int64; # UInt ??? Strange: Int32 was slower on my machine ?!? + + +#= mutable =# struct PP + expv::Vector{PP_exponent} +end + +function Base.copy(t::PP) + return PP(copy(t.expv)) +end + +function length(t::PP) + return length(t.expv) +end + +function getindex(t::PP, i::Int) + return t.expv[i] +end + +function setindex!(t::PP, i::Int, d::Int) + return t.expv[i] = d +end + +# Should be is_one, but julia complained :-( +function isone(t::PP) + return all(t.expv .== 0) +end + +function degree(t::PP) + return sum(t.expv) +end + + +# RETURN VALUE??? Perhaps index or 0? (or -1?) +function is_simple_power_pp(t::PP) + CountNZ = 0 + for i in 1:length(t) + @inbounds if (t[i] == 0) continue end + if (CountNZ > 0) return false end + CountNZ = i + end + if (CountNZ != 0) return true end # MAYBE RETURN index & exp??? + return false # because t == 1 +end + + +function is_divisible(t::PP, s::PP) # is t divisible by s + nvars = length(t) # assume equal to length(s) + for i in 1:nvars + @inbounds if t[i] < s[i] + return false + end + end + return true +end + + +# modifies first arg +function mult_by_var!(t::PP, j::Int64) + @inbounds t.expv[j] += 1 +end + +function mult(t1::PP, t2::PP) + # ASSUMES: length(t1) == length(t2) + return PP(t1.expv + t2.expv) +end + + +function divide(t1::PP, t2::PP) + # ASSUMES: length(t1) == length(t2), also that t1 is mult of t2 + return PP(t1.expv - t2.expv) +end + +function is_coprime(t1::PP, t2::PP) + # ASSUMES: length(t1) == length(t2) + nvars = length(t1) + for i in 1:nvars + @inbounds if t1[i] != 0 && t2[i] != 0 + return false + end + end + return true +end + +function lcm(t1::PP, t2::PP) + # ASSUMES: length(t1.expv) == length(t2.expv) + nvars = length(t1) + expv = [0 for _ in 1:nvars] + for i in 1:nvars + @inbounds expv[i] = max(t1[i], t2[i]) + end + return PP(expv) +end + +function gcd(t1::PP, t2::PP) + # ASSUMES: length(t1.expv) == length(t2.expv) + nvars = length(t1) + expv = [0 for _ in 1:nvars] + for i in 1:nnvars + @inbounds expv[i] = min(t1[i], t2[i]) + end + return PP(expv) +end + +function gcd3(t1::PP, t2::PP, t3::PP) + # ASSUMES: length(t1) == length(t2) == length(t3) + nvars = length(t1) + expv = [0 for _ in 1:nvars] + for i in 1:nvars + @inbounds expv[i] = min(t1[i], t2[i], t3[i]) + end + return PP(expv) +end + + +# Computes t1/gcd(t1,t2) +function colon(t1::PP, t2::PP) + # ASSUMES: length(t1) == length(t2) + nvars = length(t1) + expv = [0 for _ in 1:nvars] + for i in 1:nvars + @inbounds expv[i] = max(0, t1[i]-t2[i]) + end + return PP(expv) +end + + +# Computes t1/gcd(t1,t2^infty) +function saturatePP(t1::PP, t2::PP) + # ASSUMES: length(t1) == length(t2) + nvars = length(t1) + expv = [0 for _ in 1:nvars] + for i in 1:nvars + @inbounds if t2[i] == 0 + @inbounds expv[i] = t1[i] + end + end + return PP(expv) +end + +# Computes radical = product of vars which divide t +function radical(t::PP) + nvars = length(t) + expv = [0 for _ in 1:nvars] + for i in 1:nvars + @inbounds if t[i] > 0 + expv[i] = 1 + end + end + return PP(expv) +end + + +# Test if t1 < t1 in DegRevLex ordering; +# NB if t1 == t2 then returns false. +function deg_rev_lex_less(t1::PP, t2::PP) + d1 = degree(t1) + d2 = degree(t2) + if d1 != d2 return (d1 < d2) end + nvars = length(t1) + for i in nvars:-1:1 + if t1[i] != t2[i] + return (t1[i] > t2[i]) + end + end + return false +end + + +# Test whether PP involves at least 1 indet with index in IndexList +function involves(t::PP, IndexList::Vector{Int64}) + # ASSUMES: indexes in IndexList are all in range + for i in IndexList + @inbounds if t[i] != 0 + return true + end + end + return false +end + + +function Base.show(io::IO, t::PP) + if all(t.expv .== 0) # isone(PP) # why doesn't this work ??? + print(io, "1") + return + end + str = "" + nvars = length(t) + for i in 1:nvars + if (t[i] != 0) + if (str != "") + str = str * "*" + end + str = str * "x[$(i)]" + if (t[i] > 1) + str = str * "^$(t[i])" + end + end + end + print(io, str) +end + + +####################################################### +# Interreduction of a list of PPs + +# interreduce list of PPs; equiv find min set of gens for monoideal gen by L +function interreduce(L::Vector{PP}) + sort!(L, by=degree) + MinGens = PP[] + for t in L + discard = false + for s in MinGens + if is_divisible(t,s) + discard = true + break + end + end + if !discard + push!(MinGens, t) + end + end + return MinGens +end + + +# Is t a multiple of at least one element of L? +# Add degree truncation??? +function not_mult_of_any(L::Vector{PP}, t::PP) + for s in L + if is_divisible(t,s) + return false + end + end + return true +end + + + +# "project" PP onto sub-monoid of PPs gen by indets in indexes +function project_indets(t::PP, indexes::Vector{Int}) + return PP([t[k] for k in indexes]) +end + + +# NOT SURE THIS IS USEFUL: how many indets are really needed in the list L? +function true_num_vars(L::Vector{PP}) + if isempty(L) + return 0 + end + t = radical(L[1]) + for j in 2:length(L) + t = lcm(t, radical(L[j])) + end + return degree(t) +end + + +################################################################### +# This function is also private/local to this file. + +# Think of a graph where each vertex is labelled by an indeterminate. +# Each PP in the list L is interpreted as saying that all indeterminates +# involved in that PP are "connected". The task is to find (minimal) +# connected components of the entire graph. + +# Input: non-empty list of PPs +# Output: list of lists of var indexes, each sublist is a connected component +function connected_components(L::Vector{PP}) + ConnCompt::Vector{Vector{Int64}} = [] + nvars = length(L[1]) + IgnoreVar = [ false for _ in 1:nvars] + VarAppears = copy(IgnoreVar) + for t in L + for j in 1:nvars + @inbounds if t[j] != 0 + @inbounds VarAppears[j] = true + end + end + end + CountIgnore = 0 + for j in 1:nvars + @inbounds if !VarAppears[j] + @inbounds IgnoreVar[j] = true + CountIgnore += 1 + end + end + while CountIgnore < nvars + j = findfirst(!, IgnoreVar) # j is index of some var which appears in at least one PP + k = findfirst((t -> t[j] != 0), L) # pick some PP involving j-th var + lcm = L[k] + DoAnotherIteration = true + while DoAnotherIteration + DoAnotherIteration = false + for t in L + if is_coprime(lcm,t) continue end + s = saturatePP(t,lcm) + if isone(s) continue end + lcm = mult(lcm,s) ### lcm *= s + DoAnotherIteration = true + end + end #while + vars = filter((k -> lcm[k] > 0), 1:nvars) + # remove conn compt just found from L??? + #seems to be slower with this line ?!? L = filter((t -> is_coprime(t,lcm)), L) + push!(ConnCompt, vars) + for k in vars + IgnoreVar[k] = true + CountIgnore += 1 + end + end #while + return ConnCompt +end + + +############################################################################# +# Implementation of "Hilbert series numerator": +# main reference Bigatti 1997 (JPAA) "Computation of Hilbert-Poincare series" +#----------------------------------------------------------------------------- + +# Two base cases for HibertFn + +# Data needed: +# convention: indets are called x[1] to x[n] -- start from 1 because julia does +# Weight matrix: col k is weight of x[k] (only at topmost call) +# PP: internal repr is expv Vector{PP_exponent} equiv to Vector{Int} or similar +# Gens: +# SimplePPs: PPs which are "simple" power of an indet +# NonSimplePPs: PPs which are divisible by at least 2 distinct indets +# SimplePPs union NonSimplePPs is interreduced & non-empty!! +# HS "indets": seems useful to have power-prods of them by the corr weights -- parameter T in the fn defns + +# Case gens are simple powers +function HSNum_base_SimplePowers(SimplePPs::Vector{PP}, T::Vector{RingElemType}) where {RingElemType <: RingElem} # T is list of HSNum PPs, one for each grading dim + ans = one(T[1]) + for t in SimplePPs + k = findfirst(entry -> (entry > 0), t.expv) + ans = ans * (1 - T[k]^t[k]) # ???? ans -= ans*T[k]^t[k] + end + return ans +end + + +function HSNum_base_case1(t::PP, SimplePPs::Vector{PP}, T::Vector{RingElemType}) where {RingElemType <: RingElem} # T is list of HSNum PPs, one for each grading dim + # t is not a "simple power", all the others are + @vprintln :hilbert 1 "HSNum_base_case1: t = $(t)"; + @vprintln :hilbert 1 "HSNum_base_case1: SimplePPs = $(SimplePPs)"; + ans = HSNum_base_SimplePowers(SimplePPs, T) + ReducedSimplePPs::Vector{PP} = [] + for j in 1:length(SimplePPs) + @inbounds e = SimplePPs[j].expv + k = findfirst((entry -> (entry > 0)), e) # ispositive + if t[k] == 0 + push!(ReducedSimplePPs,SimplePPs[j]) + continue + end # no need to make a copy + tt = copy(SimplePPs[j]) + tt.expv[k] -= t[k] # guaranteed > 0 + push!(ReducedSimplePPs, tt) + end + nvars = length(t) + @inbounds scale = prod([T[k]^t[k] for k in 1:nvars]) + ans = ans - scale * HSNum_base_SimplePowers(ReducedSimplePPs, T) + @vprintln :hilbert 1 "HSNum_base_case1: returning $(ans)"; + return ans +end # function + + +## CC contains at least 2 connected components (each compt repr as Vector{Int64} of the variable indexes in the compt) +function HSNum_splitting_case(CC::Vector{Vector{Int64}}, SimplePPs::Vector{PP}, NonSimplePPs::Vector{PP}, T::Vector{RET}, PivotStrategy::Symbol) where {RET <: RingElem} + @vprintln :hilbert 1 "Splitting case: CC = $(CC)"; + HSNumList = [] ## list of HSNums + # Now find any simple PPs which are indep of the conn compts found + nvars = length(NonSimplePPs[1]) + FoundVars = PP([0 for _ in 1:nvars]) # will be prod of all vars appearing in some conn compt + for IndexSet in CC + for j in IndexSet + mult_by_var!(FoundVars, j) + end + end + IsolatedSimplePPs = filter((t -> is_coprime(t,FoundVars)), SimplePPs) + for IndexSet in CC + SubsetNonSimplePPs = filter((t -> involves(t,IndexSet)), NonSimplePPs) + SubsetSimplePPs = filter((t -> involves(t,IndexSet)), SimplePPs) + push!(HSNumList, HSNum_loop(SubsetSimplePPs, SubsetNonSimplePPs, T, PivotStrategy)) + # Next 3 lines commented out: seemed to bring no benefit (??but why??) + #? SubsetNonSimplePPs = [project_indets(t, IndexSet) for t in SubsetNonSimplePPs] + #? SubsetSimplePPs = [project_indets(t, IndexSet) for t in SubsetSimplePPs] + #? push!(HSNumList, HSNum_loop(SubsetSimplePPs, SubsetGens, [T[k] for k in IndexSet], PivotStrategy)) + end + HSNum_combined = prod(HSNumList) + if !isempty(IsolatedSimplePPs) + HSNum_combined *= HSNum_base_SimplePowers(IsolatedSimplePPs, T) + ##OLD HSNum_combined *= HSNum_loop(IsolatedSimplePPs, PP[], T, PivotStrategy) + end + return HSNum_combined +end + + +function HSNum_total_splitting_case(VarIndexes::Vector{Int64}, SimplePPs::Vector{PP}, NonSimplePPs::Vector{PP}, T::Vector{RET}, PivotStrategy::Symbol) where {RET <: RingElem} + @vprintln :hilbert 1 "Total splitting case: VarIndexes = $(VarIndexes)"; + HSNumList = [] ## list of HSNums + # Now find any simple PPs which are indep of the conn compts found + nvars = length(NonSimplePPs[1]) + FoundVars = PP([0 for _ in 1:nvars]) # will be prod of all vars appearing in some conn compt + for i in VarIndexes + mult_by_var!(FoundVars, i) + end + IsolatedSimplePPs = filter((t -> is_coprime(t,FoundVars)), SimplePPs) + for t in NonSimplePPs + SubsetSimplePPs = filter((s -> !is_coprime(s,t)), SimplePPs) + push!(HSNumList, HSNum_base_case1(t, SubsetSimplePPs, T)) + end + HSNum_combined = prod(HSNumList) + if !isempty(IsolatedSimplePPs) + HSNum_combined *= HSNum_base_SimplePowers(IsolatedSimplePPs, T) + ##OLD HSNum_combined *= HSNum_loop(IsolatedSimplePPs, PP[], T, PivotStrategy) + end + return HSNum_combined +end + + +#------------------------------------------------------------------ +# Pivot selection strategies: (see Bigatti 1997 paper) + +# term-order corr to matrix with -1 on anti-diag: [[0,0,-1], [0,-1,0],...] +function is_rev_lex_smaller(t1::PP, t2::PP) + n = length(t1) + for j in n:-1:1 + if t1[j] != t2[j] + return (t1[j] > t2[j]) + end + end + return false # t1 and t2 were equal (should not happen in this code) +end + +function rev_lex_min(L::Vector{PP}) + # assume length(L) > 0 + if isempty(L) return L[1] end + IndexMin = 1 + for j in 2:length(L) + if !is_rev_lex_smaller(L[j], L[IndexMin]) + IndexMin = j + end + end + return L[IndexMin] +end + +function rev_lex_max(L::Vector{PP}) + # assume length(L) > 0 + if length(L) == 1 return L[1] end + IndexMax = 1 + for j in 2:length(L) + if !is_rev_lex_smaller(L[IndexMax], L[j]) + IndexMax = j + end + end + return L[IndexMax] +end + +function HSNum_bayer_stillman(SimplePPs::Vector{PP}, NonSimplePPs::Vector{PP}, T::Vector{RET}) where {RET <: RingElem} + ##println("HSNum_BS: Simple: $(SimplePPs)") + ##println("HSNum_BS: NonSimple: $(NonSimplePPs)") + # Maybe sort the gens??? + if isempty(NonSimplePPs) + return HSNum_base_SimplePowers(SimplePPs, T) + end + if length(NonSimplePPs) == 1 + return HSNum_base_case1(NonSimplePPs[1], SimplePPs, T) + end + # NonSimplePPs contains at least 2 elements + # BSPivot = last(NonSimplePPs) # pick one somehow -- this is just simple to impl + #?? BSPivot = rev_lex_min(NonSimplePPs) # VERY SLOW on Hilbert-test-rnd6.jl + BSPivot = rev_lex_max(NonSimplePPs) + #VERY SLOW?!? BSPivot = rev_lex_max(vcat(SimplePPs,NonSimplePPs)) # VERY SLOW on Hilbert-test-rnd6.jl + @vprintln :hilbert 2 "BSPivot = $(BSPivot)"; + NonSPP = filter((t -> (t != BSPivot)), NonSimplePPs) + SPP = SimplePPs#filter((t -> (t != BSPivot)), SimplePPs) + part1 = HSNum_loop(SPP, NonSPP, T, :bayer_stillman) + ReducedPPs = interreduce([colon(t,BSPivot) for t in vcat(SPP,NonSPP)]) + NewSimplePPs, NewNonSimplePPs = SeparateSimplePPs(ReducedPPs) + part2 = HSNum_loop(NewSimplePPs, NewNonSimplePPs, T, :bayer_stillman) + nvars = length(BSPivot) + return part1 - prod([T[k]^BSPivot[k] for k in 1:nvars])*part2 +end + +#-------------------------------------------- +# Pivot selection strategies + +# [[AUXILIARY FN]] +# Return 2-tuple: +# (part 1) list of indexes of the indets which appear in +# the greatest number of PPs in gens. +# (part 2) corr freq +function HSNum_most_freq_indets(gens::Vector{PP}) + # ASSUMES: gens is non-empty + nvars = length(gens[1]) + freq = [0 for i in 1:nvars] # Vector{Int} or Vector{UInt} ??? + for t in gens + for i in 1:nvars + @inbounds if (t[i] != 0) + @inbounds freq[i] += 1 + end + end + end + MaxFreq = maximum(freq) + MostFreq = findall((x -> x==MaxFreq), freq) + return MostFreq, MaxFreq +end + + +# Returns index of the indet +function HSNum_most_freq_indet1(gens::Vector{PP}) + # ASSUMES: gens is non-empty + MostFreq,_ = HSNum_most_freq_indets(gens) + return MostFreq[1] +end + +# Returns index of the indet +function HSNum_most_freq_indet_rnd(gens::Vector{PP}) + # ASSUMES: gens is non-empty + MostFreq,_ = HSNum_most_freq_indets(gens) + return rand(MostFreq) +end + + + +# THIS PIVOT STRATEGY WAS AWFULLY SLOW! +# function HSNum_choose_pivot_indet(MostFreq::Vector{Int64}, gens::Vector{PP}) +# PivotIndet = rand(MostFreq) ##HSNum_most_freq_indet_rnd(gens) # or HSNum_most_freq_indet1 +# nvars = length(gens[1]) +# PivotExpv = [0 for _ in 1:nvars] +# PivotExpv[PivotIndet] = 1 +# PivotPP = PP(PivotExpv) +# end + +function HSNum_choose_pivot_simple_power_median(MostFreq::Vector{Int64}, gens::Vector{PP}) + # variant of simple-power-pivot from Bigatti JPAA, 1997 + PivotIndet = rand(MostFreq) + exps = [t[PivotIndet] for t in gens] + exps = filter((e -> e>0), exps) + sort!(exps) + exp = exps[div(1+length(exps),2)] # "median" + nvars = length(gens[1]) + PivotExpv = [0 for _ in 1:nvars] + PivotExpv[PivotIndet] = exp + PivotPP = PP(PivotExpv) +end + +function HSNum_choose_pivot_simple_power_max(MostFreq::Vector{Int64}, gens::Vector{PP}) + # variant of simple-power-pivot from Bigatti JPAA, 1997 + PivotIndet = rand(MostFreq) + exps = [t[PivotIndet] for t in gens] + exp = max(exps...) + nvars = length(gens[1]) + PivotExpv = [0 for _ in 1:nvars] + PivotExpv[PivotIndet] = exp + PivotPP = PP(PivotExpv) +end + +function HSNum_choose_pivot_gcd2simple(MostFreq::Vector{Int64}, gens::Vector{PP}) + # simple-power-pivot from Bigatti JPAA, 1997 + PivotIndet = rand(MostFreq) + cand = filter((t -> t[PivotIndet]>0), gens) + if length(cand) == 1 + error("Failed to detect total splitting case") ### !!SHOULD NEVER GET HERE!! + end + nvars = length(gens[1]) + expv = [0 for _ in 1:nvars] + expv[PivotIndet] = min(cand[1].expv[PivotIndet], cand[2].expv[PivotIndet]) + return PP(expv) +end + +function HSNum_choose_pivot_gcd2max(MostFreq::Vector{Int64}, gens::Vector{PP}) + PivotIndet = rand(MostFreq) + cand = filter((t -> t[PivotIndet]>0), gens) + if length(cand) == 1 + error("Failed to detect total splitting case") ### !!SHOULD NEVER GET HERE!! + end + pick2 = [cand[k] for k in random_subset(length(cand),2)] + t = gcd(pick2...) + nvars = length(t) + d = maximum(t.expv) + for i in 1:nvars + if t[i] < d + t.expv[i]=0 + end + end + return t +end + + +# # May produce a non-simple pivot!!! +# function HSNum_choose_pivot_gcd3(MostFreq::Vector{Int64}, gens::Vector{PP}) +# PivotIndet = rand(MostFreq) +# cand = filter((t -> t[PivotIndet]>0), gens) +# if length(cand) == 1 +# error("Failed to detect total splitting case") ### !!SHOULD NEVER GET HERE!! +# end +# if length(cand) == 2 +# return gcd(cand[1], cand[2]) +# end +# pick3 = [cand[k] for k in random_subset(length(cand),3)] +# return gcd3(pick3[1], pick3[2], pick3[3]) +# end + + +function HSNum_choose_pivot_gcd3simple(MostFreq::Vector{Int64}, gens::Vector{PP}) + PivotIndet = rand(MostFreq) + cand = filter((t -> t[PivotIndet]>0), gens) + if length(cand) == 1 + error("Failed to detect total splitting case") ### !!SHOULD NEVER GET HERE!! + end + if length(cand) == 2 + t = gcd(cand[1], cand[2]) + else + pick3 = [cand[k] for k in random_subset(length(cand),3)] + t = gcd3(pick3...) ## t = gcd3(pick3[1], pick3[2], pick3[3]) + end + # t is gcd of 3 rnd PPs (or of 2 if there are only 2) + # Now set all exps to 0 except the first max one. + j = 1 + d = t[1] + for i in 2:length(t) + if t[i] <= d + t.expv[i] = 0 + continue + end + t.expv[j] = 0 + j = i + d = t[i] + end + return t +end + + +function HSNum_choose_pivot_gcd3max(MostFreq::Vector{Int64}, gens::Vector{PP}) + PivotIndet = rand(MostFreq) + cand = filter((t -> t[PivotIndet]>0), gens) + if length(cand) == 1 + error("Failed to detect total splitting case") ### !!SHOULD NEVER GET HERE!! + end + if length(cand) == 2 + t = gcd(cand[1], cand[2]) + else + pick3 = [cand[k] for k in random_subset(length(cand),3)] + t = gcd3(pick3...) ## t = gcd3(pick3[1], pick3[2], pick3[3]) + end + d = maximum(t.expv) + # Now set to 0 all exps which are less than the max + for i in 1:length(t) + if t[i] < d + t.expv[i]=0 + end + end + return t +end + + +# # May produce a non-simple pivot!!! +# function HSNum_choose_pivot_gcd4(MostFreq::Vector{Int64}, gens::Vector{PP}) +# PivotIndet = rand(MostFreq) +# cand = filter((t -> t[PivotIndet]>0), gens) +# if length(cand) == 1 +# error("Failed to detect total splitting case") ### !!SHOULD NEVER GET HERE!! +# end +# if length(cand) == 2 +# return gcd(cand[1], cand[2]) +# end +# if length(cand) ==3 +# return gcd3(cand[1], cand[2], cand[3]) +# end +# pick4 = [cand[k] for k in random_subset(length(cand),4)] +# return gcd(gcd(pick4[1], pick4[2]), gcd(pick4[3],pick4[4])) +# end + + + +# Assume SimplePPs+NonSimplePPs are interreduced +function HSNum_loop(SimplePPs::Vector{PP}, NonSimplePPs::Vector{PP}, T::Vector{RET}, PivotStrategy::Symbol) where {RET <: RingElem} + @vprintln :hilbert 1 "HSNum_loop: SimplePPs=$(SimplePPs)" + @vprintln :hilbert 1 "HSNum_loop: NonSimplePPs=$(NonSimplePPs)" +# @vprintln :hilbert 1 "LOOP: first <=5 NonSimplePPs=$(first(NonSimplePPs,5))" + # Check if we have base case 0 + if isempty(NonSimplePPs) + @vprintln :hilbert 1 "HSNum_loop: --> delegate base case 0" + return HSNum_base_SimplePowers(SimplePPs, T) + end + # Check if we have base case 1 + if length(NonSimplePPs) == 1 + @vprintln :hilbert 1 "HSNum_loop: --> delegate base case 1" + return HSNum_base_case1(NonSimplePPs[1], SimplePPs, T) + end + # ---------------------- + MostFreq,freq = HSNum_most_freq_indets(NonSimplePPs) + if freq == 1 + return HSNum_total_splitting_case(MostFreq, SimplePPs, NonSimplePPs, T, PivotStrategy) + end + + if PivotStrategy == :bayer_stillman + return HSNum_bayer_stillman(SimplePPs, NonSimplePPs, T) + end + # Check for "splitting case" + if length(NonSimplePPs) <= length(NonSimplePPs[1])#=nvars=# + CC = connected_components(NonSimplePPs) + else CC = [] + end + if length(CC) > 1 + @vprintln :hilbert 1 "HSNum_loop: --> delegate Splitting case" + return HSNum_splitting_case(CC, SimplePPs, NonSimplePPs, T, PivotStrategy) + end + # ---------------------- + # Pivot case: first do the ideal sum, then do ideal quotient + # (the ideas are relatively simple, the code is long, tedious, and a bit fiddly) + if PivotStrategy == :simple_power_median || PivotStrategy == :auto + PivotPP = HSNum_choose_pivot_simple_power_median(MostFreq, NonSimplePPs) + end + if PivotStrategy == :simple_power_max + PivotPP = HSNum_choose_pivot_simple_power_max(MostFreq, NonSimplePPs) + end + if PivotStrategy == :gcd2simple + PivotPP = HSNum_choose_pivot_gcd2simple(MostFreq, NonSimplePPs) + end + if PivotStrategy == :gcd2max + PivotPP = HSNum_choose_pivot_gcd2max(MostFreq, NonSimplePPs) + end + # if PivotStrategy == :gcd3 + # PivotPP = HSNum_choose_pivot_gcd3(MostFreq, NonSimplePPs) + # end + if PivotStrategy == :gcd3simple + PivotPP = HSNum_choose_pivot_gcd3simple(MostFreq, NonSimplePPs) + end + if PivotStrategy == :gcd3max + PivotPP = HSNum_choose_pivot_gcd3max(MostFreq, NonSimplePPs) + end + # if PivotStrategy == :gcd4 + # PivotPP = HSNum_choose_pivot_gcd4(MostFreq, NonSimplePPs) + # end + # end + @vprintln :hilbert 1 "HSNum_loop: pivot = $(PivotPP)" + PivotIsSimple = is_simple_power_pp(PivotPP) + PivotIndex = findfirst((e -> e>0), PivotPP.expv) # used only if PivotIsSimple == true + USE_SAFE_VERSION_SUM=false + if USE_SAFE_VERSION_SUM + # Safe but slow version: just add new gen, then interreduce + RecurseSum = vcat(SimplePPs, NonSimplePPs) + push!(RecurseSum, PivotPP) + RecurseSum = interreduce(RecurseSum) + RecurseSum_SimplePPs_OLD, RecurseSum_NonSimplePPs_OLD = SeparateSimplePPs(RecurseSum) + RecurseSum_NonSimplePPs = RecurseSum_NonSimplePPs_OLD + RecurseSum_SimplePPs = RecurseSum_SimplePPs_OLD + else # use "clever" version: + # We know that SimplePPs + NonSimplePPs are already interreduced + if PivotIsSimple + RecurseSum_SimplePPs_NEW = copy(SimplePPs) + k = findfirst((t -> t[PivotIndex] > 0), SimplePPs) + if k === nothing + push!(RecurseSum_SimplePPs_NEW, PivotPP) + else + RecurseSum_SimplePPs_NEW[k] = PivotPP + end + RecurseSum_NonSimplePPs_NEW = filter((t -> t[PivotIndex] < PivotPP[PivotIndex]), NonSimplePPs) + else # PivotPP is not simple -- so this is the "general case" + RecurseSum_SimplePPs_NEW = copy(SimplePPs) # need to copy? + RecurseSum_NonSimplePPs_NEW = filter((t -> !is_divisible(t,PivotPP)), NonSimplePPs) + push!(RecurseSum_NonSimplePPs_NEW, PivotPP) + end + RecurseSum_NonSimplePPs = RecurseSum_NonSimplePPs_NEW + RecurseSum_SimplePPs = RecurseSum_SimplePPs_NEW + end # USE_SAFE_VERSION_SUM + + # Now do the quotient... + # Now get SimplePPs & NonSimplePPs for the quotient while limiting amount of interreduction + USE_SAFE_VERSION_QUOT = false + if USE_SAFE_VERSION_QUOT + #=SAFE VERSION: simpler but probably slower=# + RecurseQuot = [colon(t,PivotPP) for t in vcat(SimplePPs,NonSimplePPs)] + RecurseQuot = interreduce(RecurseQuot) + RecurseQuot_SimplePPs, RecurseQuot_NonSimplePPs = SeparateSimplePPs(RecurseQuot) + else # Use "smart version" + if !PivotIsSimple # GENERAL CASE (e.g. if not PivotIsSimple) + # Clever approach (for non-simple pivots) taken from Bigatti 1997 paper (p.247 just after Prop 1 to end of sect 5) + # ???Maybe use filter2 (see start of this file) instead of loop below??? + BM = PP[] + NotBM = PP[] + PivotPlus = mult(PivotPP, radical(PivotPP)) + for t in NonSimplePPs + if is_divisible(t,PivotPlus) + push!(BM,t) + else + push!(NotBM,t) + end + end + # BM is short for "big multiple" (see Bigatti 97, p.247 between Prop 1 and Prop 2) + BM = [divide(t,PivotPP) for t in BM] # divide is same as colon here + NotBM = vcat(NotBM, SimplePPs) + NotBM_mixed = PP[] + NotBM_coprime = PP[] + for t in NotBM + if is_coprime(t,PivotPP) + push!(NotBM_coprime,t) + else + push!(NotBM_mixed, colon(t,PivotPP)) + end + end + # At this poiint we have 3 disjoint lists of PPs: + # BM (big multiples), + # NotBM_coprime, + # NotBM_mixed + # In all cases the PPs have been colon-ed by PivotPP + NotBM_mixed = interreduce(NotBM_mixed) # cannot easily be "clever" here + filter!((t -> not_mult_of_any(NotBM_mixed,t)), NotBM_coprime) + RecurseQuot = vcat(NotBM_coprime, NotBM_mixed) # already interreduced + RQ_SimplePPs, RQ_NonSimplePPs = SeparateSimplePPs(RecurseQuot) + RecurseQuot_SimplePPs = RQ_SimplePPs + RecurseQuot_NonSimplePPs = vcat(BM, RQ_NonSimplePPs) + else # Clever approach when PivotIsSimple + # The idea behind this code is fairly simple; sadly the code itself is not :-( + RecurseQuot_SimplePPs = copy(SimplePPs) + k = findfirst((t -> t[PivotIndex] > 0), RecurseQuot_SimplePPs) + if !(k === nothing) + RecurseQuot_SimplePPs[k] = copy(RecurseQuot_SimplePPs[k]) + RecurseQuot_SimplePPs[k].expv[PivotIndex] -= PivotPP[PivotIndex] + end + DegPivot = degree(PivotPP) + NonSimpleTbl = [PP[] for _ in 0:DegPivot] ## WARNING: indexes are offset by 1 -- thanks Julia! + NonSimple1 = PP[] # will contain all PPs divisible by PivotPP^(1+epsilon) + for t in NonSimplePPs + degt = t[PivotIndex] + if degt > DegPivot + push!(NonSimple1, divide(t, PivotPP)) + else + push!(NonSimpleTbl[degt+1], colon(t,PivotPP)) + end + end + NonSimple2 = NonSimpleTbl[DegPivot+1] + for i in DegPivot:-1:1 + NewPPs = filter((t -> not_mult_of_any(NonSimple2,t)), NonSimpleTbl[i]) + NonSimple2 = vcat(NonSimple2, NewPPs) + end + NewSimplePPs = filter(is_simple_power_pp, NonSimple2) ## Use instead + NonSimple2 = filter(!is_simple_power_pp, NonSimple2) ## SeparateSimplePPs??? + if !isempty(NewSimplePPs) + RecurseQuot_SimplePPs = interreduce(vcat(RecurseQuot_SimplePPs, NewSimplePPs)) + end + RecurseQuot_NonSimplePPs = vcat(NonSimple1, NonSimple2) + end #PivotIsSimple + end # USE_SAFE_VERSION_QUOT + # Now put the two pieces together: + nvars = length(PivotPP) + scale = prod([T[k]^PivotPP[k] for k in 1:nvars]) + @vprintln :hilbert 2 "HSNum_loop: SUM recursion: simple $(RecurseSum_SimplePPs)"; + @vprintln :hilbert 2 "HSNum_loop: SUM recursion: nonsimple $(RecurseSum_NonSimplePPs)"; + @vprintln :hilbert 2 "HSNum_loop: QUOT recursion: simple $(RecurseQuot_SimplePPs)"; + @vprintln :hilbert 2 "HSNum_loop: QUOT recursion: nonsimple $(RecurseQuot_NonSimplePPs)"; + + HSNum_sum = HSNum_loop(RecurseSum_SimplePPs, RecurseSum_NonSimplePPs, T, PivotStrategy) + HSNum_quot = HSNum_loop(RecurseQuot_SimplePPs, RecurseQuot_NonSimplePPs, T, PivotStrategy) + @vprintln :hilbert 1 "HSNum_loop: END OF CALL"; + return HSNum_sum + scale*HSNum_quot +end + + +function separate_simple_pps(gens::Vector{PP}) + SimplePPs::Vector{PP} = [] + NonSimplePPs::Vector{PP} = [] + for g in gens + if is_simple_power_pp(g) + push!(SimplePPs, g) + else + push!(NonSimplePPs, g) + end + end + return SimplePPs, NonSimplePPs +end + + +# Check args: either throws or returns nothing. +function HSNum_check_args(gens::Vector{PP}, W::Vector{Vector{Int}}) + if isempty(W) + throw("HSNum: weight matrix must have at least 1 row") + end + nvars = length(W[1]) + if !all((row -> length(row)==nvars), W) + throw("HSNum: weight matrix must have 1 column for each variable") + end + if !all((t -> length(t)==nvars), gens) # OK also if isempty(gens) + throw("HSNum: generators must all have same size exponent vectors") + end + # Args are OK, so simply return (without throwing) +end + + + +function HSNum_abbott_PPs(PP_gens::Vector{PP}, W::Vector{Vector{Int}}, PivotStrategy::Symbol, HSRing::Ring) + # ASSUME W is "rectangular" + @vprintln :hilbert 1 "HSNum: PP_gens = $(PP_gens)"; + @vprintln :hilbert 1 "HSNum: weight matrix W = $(W)"; + HSNum_check_args(PP_gens, W) #throws or does nothing + # Grading is over ZZ^m + m = length(W) + ncols = length(W[1]) + nvars = ncols + t = gens(HSRing) + @assert length(t) >= m "supplied Hilbert series ring contains too few variables" + T = [one(HSRing) for k in 1:nvars] + for k in 1:nvars + s = one(HSRing) + for j in 1:m + s *= t[j]^W[j][k] + end + T[k] = s + end + # Now have T[1] = t1^W[1,1] * t2^W[2,1] * ..., etc + SimplePPs,NonSimplePPs = separate_simple_pps(interreduce(PP_gens)) + sort!(NonSimplePPs, lt=deg_rev_lex_less) # recommended by Bayer+Stillman (for their criterion) + return HSNum_loop(SimplePPs, NonSimplePPs, T, PivotStrategy) +end + +#----------------------------------------------------------------------------- +# This fn copied from GradedModule.jl (in dir OSCAR/HILBERT/) +function gen_repr(d) + grading_dim = length(gens(parent(d))) + return [d[k] for k in 1:grading_dim] +end + +@doc raw""" + HSNum_abbott(A::MPolyQuoRing, HSRing::Ring; pivot_strategy::Symbol = :auto) + +Don't use this internal function: use `hilbert_series` or `multi_hilbert_series` instead. +Compute numerator of Hilbert series of the quotient `A`. +Result is a pair: `N, D` being the numerator `N` (as a laurent polynomial) and the denominator `D` as +a factor list of laurent polynomials. + +!!! note + Applied to an ideal `I`, the function first homogenizes the generators of `I` in the extended ring. + It then creates the ideal generated by these homogenizations, and saturates this ideal + with respect to the ideal which is generated by the homogenizing variables. + +# Examples +```jldoctest +julia> R, (x,y,z) = graded_polynomial_ring(QQ, ["x", "y","z"]) +(Graded multivariate polynomial ring in 3 variables over QQ, MPolyDecRingElem{QQFieldElem, QQMPolyRingElem}[x, y, z]) + +julia> I = ideal(R, [x^3+y^2*z, y^3+x*z^2, z^3+x^2*y]); + +julia> RmodI,_ = quo(R,I); + +julia> HSRing1,_ = polynomial_ring(ZZ, "t"); + +julia> Oscar.HSNum_abbott(RmodI, HSRing1) +-t^9 + 3*t^6 - 3*t^3 + 1 + +julia> R, (x,y,z) = graded_polynomial_ring(QQ, ["x", "y","z"], [1 2 3; 3 2 1]) +(Graded multivariate polynomial ring in 3 variables over QQ, MPolyDecRingElem{QQFieldElem, QQMPolyRingElem}[x, y, z]) + +julia> I = ideal(R, [x*z+y^2, y^6+x^3*z^3, z^6, x^6]); + +julia> RmodI,_ = quo(R,I); + +julia> HSRing2,_ = polynomial_ring(ZZ, ["t[1]", "t[2]"]); + +julia> Oscar.HSNum_abbott(RmodI, HSRing2) +-t[1]^28*t[2]^28 + t[1]^24*t[2]^24 + t[1]^22*t[2]^10 - t[1]^18*t[2]^6 + t[1]^10*t[2]^22 - t[1]^6*t[2]^18 - t[1]^4*t[2]^4 + 1 + +``` +""" +function HSNum_abbott(PmodI::MPolyQuoRing, HSRing::Ring; pivot_strategy::Symbol = :auto) + I = modulus(PmodI) + P = base_ring(I) + nvars = length(gens(P)) + grading_dim = length(gens(Oscar.parent(degree(gen(P,1))))) # better way??? + weights = [degree(var) for var in gens(P)] + W = [[0 for _ in 1:nvars] for _ in 1:grading_dim] + for i in 1:nvars + expv = [Int64(exp) for exp in gen_repr(degree(gen(P,i)))] + for j in 1:grading_dim + W[j][i] = expv[j] + end + end + LTs = gens(leading_ideal(I)) + PPs = [PP(degrees(t)) for t in LTs] + return HSNum_abbott_PPs(PPs, W, pivot_strategy, HSRing) +end + + + +# W[k] is ZZ^m-grading of x[k] +# Expect that HSRing is PolynomialRing or LaurentPolynomialRing +function _hilbert_series_denominator(HSRing::Ring, W::Vector{Vector{Int}}) + n = length(W) + m = length(W[1]) + @req all(r -> length(r)==m, W) "Grading VecVec must be rectangular" + @req length(gens(HSRing)) >= m "Hilbert series ring has too few variables" + t = gens(HSRing) + fac_dict = Dict{elem_type(HSRing), Integer}() + for i = 1:n + # adjoin factor ( 1 - prod(t_j^W[i,j]) ) + new_fac = 1 - prod([t[k]^W[i][k] for k in 1:m]); + # ???BUG??? DOES NOT WORK if HSRing is Univariate poly ring over ZZ + # B = MPolyBuildCtx(HSRing) + # push_term!(B, 1, W[i]) + # new_fac = 1-finish(B) + if haskey(fac_dict, new_fac) + fac_dict[new_fac] += 1 # coalesce identical factors + else + fac_dict[new_fac] = 1 + end + end + return FacElem(HSRing, fac_dict) +end + + +function _hilbert_series_ring(parent::Union{Nothing, Ring}, m::Int) + if parent !== nothing + # Caller supplied a Ring; check that it is reasonable. + # The following line might complain in unforeseen ways in case + # the argument for `parent` was too unreasonable. However, we would + # like to keep the possibilities for that input rather broad. + @req ngens(parent) >= m "parent ring of the output does not contain sufficiently many variables" + return parent + end + # Caller did not supply a Ring, so we create one. + if m == 1 + VAR = [:t] + else + VAR = [_make_variable("t", i) for i = 1:m] + end + HSRing, _ = LaurentPolynomialRing(ZZ, VAR) + return HSRing +end + +# !!! Note +# We create a Laurent poly ring even if all weights are non-negative +# because negative exponents may be needed when dealing with modules +# with shifts. diff --git a/src/Rings/hilbert_zach.jl b/src/Rings/hilbert_zach.jl new file mode 100644 index 000000000000..fe953076ba79 --- /dev/null +++ b/src/Rings/hilbert_zach.jl @@ -0,0 +1,442 @@ +######################################################################## +# Backend functionality for computation of multivariate Hilbert series +# due to Matthias Zach +######################################################################## + + +####################################################################### +# 06.10.2022 +# The following internal routine can probably still be tuned. +# +# Compared to the implementation in Singular, it is still about 10 +# times slower. We spotted the following possible deficits: +# +# 1) The returned polynomials are not yet built with MPolyBuildCtx. +# We tried that, but as for now, the code for the build context +# is even slower than the direct implementation that is in place +# now. Once MPolyBuildCtx is tuned, we should try that again; +# the code snippets are still there. In total, building the return +# values takes up more than one third of the computation time +# at this moment. +# +# 2) Numerous allocations for integer vectors. The code below +# performs lots of iterative allocations for lists of integer +# vectors. If we constructed a container data structure to +# maintain this list internally and do the allocations at once +# (for instance in a big matrix), this could significantly +# speed up the code. However, that is too much work for +# the time being, as long as a high-performing Hilbert series +# computation is not of technical importance. +# +# 3) Singular uses bitmasking for exponent vectors to decide on +# divisibility more quickly. This is particularly important in +# the method _divide_by. For singular, this leads to a speedup +# of factor 5. Here, only less than one third of the time is +# spent in `_divide_by`, so it does not seem overly important +# at this point, but might become relevant in the future. +# A particular modification of the singular version of bitmasking +# is to compute means for the exponents occurring for each variable +# and set the bits depending on whether a given exponent is greater +# or less than that mean value. + +function _hilbert_numerator_from_leading_exponents( + a::Vector{Vector{Int}}, + weight_matrix::Matrix{Int}, + return_ring::Ring, + #algorithm=:generator + #algorithm=:custom + #algorithm=:gcd + #algorithm=:indeterminate + #algorithm=:cocoa + algorithm::Symbol # =:BayerStillmanA, # This is by far the fastest strategy. Should be used. + # short exponent vectors where the k-th bit indicates that the k-th + # exponent is non-zero. + ) + t = gens(return_ring) + ret = _hilbert_numerator_trivial_cases(a, weight_matrix, return_ring, t) + ret !== nothing && return ret + + if algorithm == :BayerStillmanA + return _hilbert_numerator_bayer_stillman(a, weight_matrix, return_ring, t) + elseif algorithm == :custom + return _hilbert_numerator_custom(a, weight_matrix, return_ring, t) + elseif algorithm == :gcd # see Remark 5.3.11 + return _hilbert_numerator_gcd(a, weight_matrix, return_ring, t) + elseif algorithm == :generator # just choosing on random generator, cf. Remark 5.3.8 + return _hilbert_numerator_generator(a, weight_matrix, return_ring, t) + elseif algorithm == :indeterminate # see Remark 5.3.8 + return _hilbert_numerator_indeterminate(a, weight_matrix, return_ring, t) + elseif algorithm == :cocoa # see Remark 5.3.14 + return _hilbert_numerator_cocoa(a, weight_matrix, return_ring, t) + end + error("invalid algorithm") +end + + +######################################################################## +# Implementations of the different algorithms below +# +# Compute the Hilbert series from the monomial leading ideal using +# different bisection strategies along the lines of +# [Kreuzer, Robbiano: Computational Commutative Algebra 2, Springer] +# Section 5.3. +######################################################################## + +function _hilbert_numerator_trivial_cases( + a::Vector{Vector{Int}}, weight_matrix::Matrix{Int}, + S::Ring, t::Vector, oneS = one(S) + ) + length(a) == 0 && return oneS + + # See Proposition 5.3.6 + if _are_pairwise_coprime(a) + return prod(oneS - _expvec_to_poly(S, t, weight_matrix, e) for e in a) + end + + return nothing +end + +function _hilbert_numerator_cocoa( + a::Vector{Vector{Int}}, weight_matrix::Matrix{Int}, + S::Ring, t::Vector + ) + ret = _hilbert_numerator_trivial_cases(a, weight_matrix, S, t) + ret !== nothing && return ret + + n = length(a) + m = length(a[1]) + counters = [0 for i in 1:m] + for e in a + counters += [iszero(k) ? 0 : 1 for k in e] + end + j = _find_maximum(counters) + + p = a[rand(1:n)] + while p[j] == 0 + p = a[rand(1:n)] + end + + q = a[rand(1:n)] + while q[j] == 0 || p == q + q = a[rand(1:n)] + end + + pivot = [0 for i in 1:m] + pivot[j] = minimum([p[j], q[j]]) + + + ### Assembly of the quotient ideal with less generators + rhs = [e for e in a if !_divides(e, pivot)] + push!(rhs, pivot) + + ### Assembly of the division ideal with less total degree + lhs = _divide_by_monomial_power(a, j, pivot[j]) + + f = one(S) + for i in 1:nvars(S) + z = t[i] + f *= z^(sum([pivot[j]*weight_matrix[i, j] for j in 1:length(pivot)])) + end + + return _hilbert_numerator_cocoa(rhs, weight_matrix, S, t) + f*_hilbert_numerator_cocoa(lhs, weight_matrix, S, t) +end + +function _hilbert_numerator_indeterminate( + a::Vector{Vector{Int}}, weight_matrix::Matrix{Int}, + S::Ring, t::Vector + ) + ret = _hilbert_numerator_trivial_cases(a, weight_matrix, S, t) + ret !== nothing && return ret + + e = first(a) + found_at = findfirst(!iszero, e)::Int64 + pivot = zero(e) + pivot[found_at] = 1 + + ### Assembly of the quotient ideal with less generators + rhs = [e for e in a if e[found_at] == 0] + push!(rhs, pivot) + + ### Assembly of the division ideal with less total degree + lhs = _divide_by(a, pivot) + + f = one(S) + for i in 1:nvars(S) + z = t[i] + f *= z^(sum([pivot[j]*weight_matrix[i, j] for j in 1:length(pivot)])) + end + + return _hilbert_numerator_indeterminate(rhs, weight_matrix, S, t) + f*_hilbert_numerator_indeterminate(lhs, weight_matrix, S, t) +end + +function _hilbert_numerator_generator( + a::Vector{Vector{Int}}, weight_matrix::Matrix{Int}, + S::Ring, t::Vector + ) + ret = _hilbert_numerator_trivial_cases(a, weight_matrix, S, t) + ret !== nothing && return ret + + b = copy(a) + pivot = pop!(b) + + f = one(S) + for i in 1:nvars(S) + z = t[i] + f *= z^(sum([pivot[j]*weight_matrix[i, j] for j in 1:length(pivot)])) + end + + c = _divide_by(b, pivot) + p1 = _hilbert_numerator_generator(b, weight_matrix, S, t) + p2 = _hilbert_numerator_generator(c, weight_matrix, S, t) + + return p1 - f * p2 +end + +function _hilbert_numerator_gcd( + a::Vector{Vector{Int}}, weight_matrix::Matrix{Int}, + S::Ring, t::Vector + ) + ret = _hilbert_numerator_trivial_cases(a, weight_matrix, S, t) + ret !== nothing && return ret + + n = length(a) + counters = [0 for i in 1:length(a[1])] + for e in a + counters += [iszero(k) ? 0 : 1 for k in e] + end + j = _find_maximum(counters) + + p = a[rand(1:n)] + while p[j] == 0 + p = a[rand(1:n)] + end + + q = a[rand(1:n)] + while q[j] == 0 || p == q + q = a[rand(1:n)] + end + + pivot = _gcd(p, q) + + ### Assembly of the quotient ideal with less generators + rhs = [e for e in a if !_divides(e, pivot)] + push!(rhs, pivot) + + ### Assembly of the division ideal with less total degree + lhs = _divide_by(a, pivot) + + f = one(S) + for i in 1:nvars(S) + z = t[i] + f *= z^(sum([pivot[j]*weight_matrix[i, j] for j in 1:length(pivot)])) + end + + return _hilbert_numerator_gcd(rhs, weight_matrix, S, t) + f*_hilbert_numerator_gcd(lhs, weight_matrix, S, t) +end + +function _hilbert_numerator_custom( + a::Vector{Vector{Int}}, weight_matrix::Matrix{Int}, + S::Ring, t::Vector + ) + ret = _hilbert_numerator_trivial_cases(a, weight_matrix, S, t) + ret !== nothing && return ret + + p = Vector{Int}() + q = Vector{Int}() + max_deg = 0 + for i in 1:length(a) + b = a[i] + for j in i+1:length(a) + c = a[j] + r = _gcd(b, c) + if sum(r) > max_deg + max_deg = sum(r) + p = b + q = c + end + end + end + + ### Assembly of the quotient ideal with less generators + pivot = _gcd(p, q) + rhs = [e for e in a if !_divides(e, pivot)] + push!(rhs, pivot) + + ### Assembly of the division ideal with less total degree + lhs = _divide_by(a, pivot) + + f = one(S) + for i in 1:nvars(S) + z = t[i] + f *= z^(sum([pivot[j]*weight_matrix[i, j] for j in 1:length(pivot)])) + end + + return _hilbert_numerator_custom(rhs, weight_matrix, S, t) + f*_hilbert_numerator_custom(lhs, weight_matrix, S, t) +end + +function _hilbert_numerator_bayer_stillman( + a::Vector{Vector{Int}}, weight_matrix::Matrix{Int}, + S::Ring, t::Vector, oneS = one(S) + ) + ########################################################################### + # For this strategy see + # + # Bayer, Stillman: Computation of Hilber Series + # J. Symbolic Computation (1992) No. 14, pp. 31--50 + # + # Algorithm 2.6, page 35 + ########################################################################### + ret = _hilbert_numerator_trivial_cases(a, weight_matrix, S, t, oneS) + ret !== nothing && return ret + + # make sure we have lexicographically ordered monomials + sort!(a, alg=QuickSort) + + # initialize the result + h = oneS - _expvec_to_poly(S, t, weight_matrix, a[1]) + linear_mons = Vector{Int}() + for i in 2:length(a) + J = _divide_by(a[1:i-1], a[i]) + empty!(linear_mons) + J1 = filter!(J) do m + k = findfirst(!iszero, m) + k === nothing && return false # filter out zero vector + if m[k] == 1 && findnext(!iszero, m, k + 1) === nothing + push!(linear_mons, k) + return false + end + return true + end + q = _hilbert_numerator_bayer_stillman(J1, weight_matrix, S, t, oneS) + for k in linear_mons + q *= (oneS - prod(t[i]^weight_matrix[i, k] for i in 1:length(t))) + end + + h -= q * _expvec_to_poly(S, t, weight_matrix, a[i]) + end + return h +end + + +######################################################################## +# Auxiliary helper functions below +######################################################################## + +### compute t ^ (weight_matrix * expvec), where t == gens(S) +function _expvec_to_poly(S::Ring, t::Vector, weight_matrix::Matrix{Int}, expvec::Vector{Int}) + o = one(coefficient_ring(S)) # TODO: cache this ?! + return S([o], [weight_matrix * expvec]) +end + +### special case for univariate polynomial ring +function _expvec_to_poly(S::PolyRing, t::Vector, weight_matrix::Matrix{Int}, expvec::Vector{Int}) + @assert length(t) == 1 + @assert size(weight_matrix) == (1, length(expvec)) + + # compute the dot-product of weight_matrix[1,:] and expvec, but faster than dot + # TODO: what about overflows? + s = weight_matrix[1,1] * expvec[1] + for i in 2:length(expvec) + @inbounds s += weight_matrix[1,i] * expvec[i] + end + return t[1]^s +end + +function _find_maximum(a::Vector{Int}) + m = a[1] + j = 1 + for i in 1:length(a) + if a[i] > m + m = a[i] + j = i + end + end + return j +end + +### This implements the speedup from the CoCoA strategy, +# see Exercise 4 in Section 5.3 +function _divide_by_monomial_power(a::Vector{Vector{Int}}, j::Int, k::Int) + divides(a::Vector{Int}, b::Vector{Int}) = all(k->(k>=0), b[1:j-1] - a[1:j-1]) && all(k->(k>=0), b[j+1:end] - a[j+1:end]) + kept = [e for e in a if e[j] == 0] + for i in 1:k + next_slice = [e for e in a if e[j] == i] + still_kept = next_slice + for e in kept + all(v->!divides(v, e), next_slice) && push!(still_kept, e) + end + kept = still_kept + end + still_kept = vcat(still_kept, [e for e in a if e[j] > k]) + + result = Vector{Vector{Int}}() + for e in still_kept + v = copy(e) + v[j] = (e[j] > k ? e[j] - k : 0) + push!(result, v) + end + return result +end + +_divides(a::Vector{Int}, b::Vector{Int}) = all(k->(a[k]>=b[k]), 1:length(a)) + +function _gcd(a::Vector{Int}, b::Vector{Int}) + return [a[k] > b[k] ? b[k] : a[k] for k in 1:length(a)] +end + +function _are_pairwise_coprime(a::Vector{Vector{Int}}) + length(a) <= 1 && return true + n = length(a) + m = length(a[1]) + for i in 1:n-1 + for j in i+1:n + all(k -> a[i][k] == 0 || a[j][k] == 0, 1:m) || return false + end + end + return true +end + +### Assume that a is minimal. We divide by `pivot` and return +# a minimal set of generators for the resulting monomial ideal. +function _divide_by(a::Vector{Vector{Int}}, pivot::Vector{Int}) + # The good ones will contribute to a minimal generating + # set of the lhs ideal. + # + # The bad monomials come from those which hop over the boundaries of + # the monomial diagram by the shift. Their span has a new + # generator which is collected in `bad` a priori. It is checked + # whether they become superfluous and if not, they are added to + # the good ones. + good = sizehint!(Vector{Vector{Int}}(), length(a)) + bad = sizehint!(Vector{Vector{Int}}(), length(a)) + for e in a + if _divides(e, pivot) + push!(good, e - pivot) + else + push!(bad, e) + end + end + + # pre-allocate m so that don't need to allocate it again each loop iteration + m = similar(pivot) + for e in bad + # the next line computers m = [k < 0 ? 0 : k for k in e] + # but without allocations + for i in 1:length(e) + m[i] = max(e[i] - pivot[i], 0) + end + + # check whether the new monomial m is already in the span + # of the good ones. If yes, discard it. If not, discard those + # elements of the good ones that are in the span of m and put + # m in the list of good ones instead. + if all(x->!_divides(m, x), good) + # Remove those 'good' elements which are multiples of m + filter!(x -> !_divides(x, m), good) + push!(good, copy(m)) + end + end + + return good +end + diff --git a/src/Rings/mpoly-affine-algebras.jl b/src/Rings/mpoly-affine-algebras.jl index 567a8987faf5..300c4715c43d 100644 --- a/src/Rings/mpoly-affine-algebras.jl +++ b/src/Rings/mpoly-affine-algebras.jl @@ -179,8 +179,10 @@ end ################################################################################## +# TODO: The function below now also works for rings which are not standard graded +# by virtue of Abbott's implementation. Clean up the docstring accordingly. @doc raw""" - hilbert_series(A::MPolyQuoRing) + hilbert_series(A::MPolyQuoRing; backend::Symbol=:Singular, algorithm::Symbol=:BayerStillmanA) Given a $\mathbb Z$-graded affine algebra $A = R/I$ over a field $K$, where the grading is inherited from a $\mathbb Z$-grading on the polynomial ring $R$ defined by assigning @@ -194,6 +196,11 @@ where $n$ is the number of variables of $R$, and $w_1, \dots, w_n$ are the assig See also `hilbert_series_reduced`. +!!! note + The advanced user can select different backends for the computation (`:Singular` and + `:Abbott` for the moment), as well as different algorithms. The latter might be + ignored for certain backends. + # Examples ```jldoctest julia> R, (w, x, y, z) = graded_polynomial_ring(QQ, ["w", "x", "y", "z"]); @@ -201,31 +208,39 @@ julia> R, (w, x, y, z) = graded_polynomial_ring(QQ, ["w", "x", "y", "z"]); julia> A, _ = quo(R, ideal(R, [w*y-x^2, w*z-x*y, x*z-y^2])); julia> hilbert_series(A) -(2*t^3 - 3*t^2 + 1, t^4 - 4*t^3 + 6*t^2 - 4*t + 1) +(2*t^3 - 3*t^2 + 1, Factored element with data +Dict{ZZPolyRingElem, ZZRingElem}(-t + 1 => 4)) julia> R, (x, y, z) = graded_polynomial_ring(QQ, ["x", "y", "z"], [1, 2, 3]); julia> A, _ = quo(R, ideal(R, [x*y*z])); julia> hilbert_series(A) -(-t^6 + 1, -t^6 + t^5 + t^4 - t^2 - t + 1) +(-t^6 + 1, Factored element with data +Dict{ZZPolyRingElem, ZZRingElem}(-t^2 + 1 => 1, -t + 1 => 1, -t^3 + 1 => 1)) ``` """ -function hilbert_series(A::MPolyQuoRing) - if iszero(A.I) - R = base_ring(A.I) - @req is_z_graded(R) "The base ring must be ZZ-graded" - W = R.d - W = [Int(W[i][1]) for i = 1:ngens(R)] - @req minimum(W) > 0 "The weights must be positive" - Zt, t = ZZ["t"] - den = prod([1-t^Int(w[1]) for w in R.d]) - return (one(parent(t)), den) - end - H = HilbertData(A.I) - return hilbert_series(H) +function hilbert_series(A::MPolyQuoRing; #=backend::Symbol=:Singular, algorithm::Symbol=:BayerStillmanA,=# parent::Union{Nothing,Ring}=nothing) + R = base_ring(A.I) + @req is_z_graded(R) "ring must be graded by the integers" + parent, t = (parent === nothing) ? polynomial_ring(ZZ, "t") : (parent, first(gens(parent))); + W = R.d + W = [Int(W[i][1]) for i = 1:ngens(R)] + @req minimum(W) > 0 "The weights must be positive" + # if iszero(A.I) + # den = prod([1-t^Int(w[1]) for w in R.d]) + # return (one(parent(t)), den) + # end + + (numer, denom), _ = multi_hilbert_series(A; parent=parent) + return numer,denom end +# TODO: The method below is missing. It should be made better and put to the correct place (AA). +ngens(S::AbstractAlgebra.Generic.LaurentMPolyWrapRing) = length(gens(S)) +ngens(S::AbstractAlgebra.Generic.LaurentPolyWrapRing) = 1 +ngens(P::PolyRing) = 1 + @doc raw""" hilbert_series_reduced(A::MPolyQuoRing) @@ -252,7 +267,8 @@ julia> R, (x, y, z) = graded_polynomial_ring(QQ, ["x", "y", "z"], [1, 2, 3]); julia> A, _ = quo(R, ideal(R, [x*y*z])); julia> hilbert_series(A) -(-t^6 + 1, -t^6 + t^5 + t^4 - t^2 - t + 1) +(-t^6 + 1, Factored element with data +Dict{ZZPolyRingElem, ZZRingElem}(-t^2 + 1 => 1, -t + 1 => 1, -t^3 + 1 => 1)) julia> hilbert_series_reduced(A) (t^2 - t + 1, t^2 - 2*t + 1) @@ -291,18 +307,18 @@ julia> hilbert_series_expanded(A, 5) ``` """ function hilbert_series_expanded(A::MPolyQuoRing, d::Int) - if iszero(A.I) - R = base_ring(A.I) - @req is_z_graded(R) "The base ring must be ZZ-graded" - W = R.d - W = [Int(W[i][1]) for i = 1:ngens(R)] - @req minimum(W) > 0 "The weights must be positive" - H = hilbert_series(A) - T, t = power_series_ring(QQ, d+1, "t") - return _rational_function_to_power_series(T, H[1], H[2]) - end - H = HilbertData(A.I) - return hilbert_series_expanded(H, d) + if iszero(modulus(A)) + R = base_ring(A) + @req is_z_graded(R) "The base ring must be ZZ-graded" + W = R.d + W = [Int(W[i][1]) for i = 1:ngens(R)] + @req minimum(W) > 0 "The weights must be positive" + num, denom = hilbert_series(A) + T, t = power_series_ring(QQ, d+1, "t") + return _rational_function_to_power_series(T, num, evaluate(denom)) + end + H = HilbertData(A.I) + return hilbert_series_expanded(H, d) end @doc raw""" @@ -453,12 +469,12 @@ end @doc raw""" - multi_hilbert_series(A::MPolyQuoRing; algorithm::Symbol=:BayerStillmanA) + multi_hilbert_series(A::MPolyQuoRing; algorithm::Symbol=:BayerStillmanA, parent::Union{Nothing,Ring}=nothing) -Return the Hilbert series of the positively graded affine algebra `A`. +Return the Hilbert series of the graded affine algebra `A`. !!! note - The advanced user can select a `algorithm` for the computation; + The advanced user can select an `algorithm` for the computation; see the code for details. # Examples @@ -478,7 +494,8 @@ julia> H[1][1] -t[1]^7*t[2]^-2 + t[1]^6*t[2]^-1 + t[1]^6*t[2]^-2 + t[1]^5*t[2]^-4 - t[1]^4 + t[1]^4*t[2]^-2 - t[1]^4*t[2]^-4 - t[1]^3*t[2]^-1 - t[1]^3*t[2]^-2 + 1 julia> H[1][2] --t[1]^3*t[2]^-1 + t[1]^2 + 2*t[1]^2*t[2]^-1 - 2*t[1] - t[1]*t[2]^-1 + 1 +Factored element with data +Dict{AbstractAlgebra.Generic.LaurentMPolyWrap{ZZRingElem, ZZMPolyRingElem, AbstractAlgebra.Generic.LaurentMPolyWrapRing{ZZRingElem, ZZMPolyRing}}, ZZRingElem}(-t[1] + 1 => 2, -t[1]*t[2]^-1 + 1 => 1) julia> H[2][1] GrpAb: Z^2 @@ -501,78 +518,96 @@ julia> R, (w, x, y, z) = graded_polynomial_ring(QQ, ["w", "x", "y", "z"], W); julia> A, _ = quo(R, ideal(R, [w*y-x^2, w*z-x*y, x*z-y^2])); -julia> H = multi_hilbert_series(A); +julia> (num, den), (H, iso) = multi_hilbert_series(A); -julia> H[1][1] +julia> num 2*t^3 - 3*t^2 + 1 -julia> H[1][2] -t^4 - 4*t^3 + 6*t^2 - 4*t + 1 +julia> den +Factored element with data +Dict{AbstractAlgebra.Generic.LaurentMPolyWrap{ZZRingElem, ZZMPolyRingElem, AbstractAlgebra.Generic.LaurentMPolyWrapRing{ZZRingElem, ZZMPolyRing}}, ZZRingElem}(-t + 1 => 4) -julia> H[2][1] +julia> H GrpAb: Z -julia> H[2][2] +julia> iso Map with following data Domain: ======= -Abelian group with structure: Z +H Codomain: ========= G ``` """ -function multi_hilbert_series(A::MPolyQuoRing; algorithm::Symbol=:BayerStillmanA) - R = base_ring(A) - I = A.I - @req coefficient_ring(R) isa AbstractAlgebra.Field "The coefficient ring must be a field" - @req is_positively_graded(R) "The base ring must be positively graded" - - G = grading_group(R) - if !is_zm_graded(R) - H, iso = snf(G) - V = [preimage(iso, x) for x in gens(G)] - isoinv = hom(G, H, V) - W = R.d - W = [isoinv(W[i]) for i = 1:length(W)] - S, _ = graded_polynomial_ring(coefficient_ring(R), symbols(R), W) - change = hom(R, S, gens(S)) - I = change(A.I) - R = S - else - H, iso = G, identity_map(G) - end - m = ngens(grading_group(R)) - n = ngens(R) - W = R.d - MI = Matrix{Int}(undef, n, m) - for i=1:n - for j=1:m - MI[i, j] = Int(W[i][j]) - end - end - if m == 1 - VAR = [:t] - else - VAR = [_make_variable("t", i) for i = 1:m] - end - S, _ = LaurentPolynomialRing(ZZ, VAR) - q = one(S) - for i = 1:n - e = [Int(MI[i, :][j]) for j = 1:m] - B = MPolyBuildCtx(S) - push_term!(B, 1, e) - q = q*(1-finish(B)) - end - if iszero(I) - p = one(S) - else - LI = leading_ideal(I, ordering=degrevlex(gens(R))) - p = _numerator_monomial_multi_hilbert_series(LI, S, m, algorithm=algorithm) - end - return (p, q), (H, iso) +function multi_hilbert_series( + A::MPolyQuoRing; + algorithm::Symbol=:BayerStillmanA, + backend::Symbol=:Abbott, + parent::Union{Nothing, Ring}=nothing + ) + R = base_ring(A) + I = modulus(A) + @req coefficient_ring(R) isa AbstractAlgebra.Field "The coefficient ring must be a field" + + # Wrap the case where G is abstractly isomorphic to ℤᵐ, but not realized as a + # free Abelian group. + # + # We use the Smith normal form to get there, recreate the graded ring with the + # free grading group, do the computation there and return the isomorphism for + # the grading. + G = grading_group(R) + if !is_zm_graded(R) + H, iso = snf(G) + V = [preimage(iso, x) for x in gens(G)] + isoinv = hom(G, H, V) + W = [isoinv(R.d[i]) for i = 1:length(R.d)] + S, _ = graded_polynomial_ring(coefficient_ring(R), symbols(R), W) + map_into_S = hom(R, S, gens(S)) + J = map_into_S(I) + AA, _ = quo(S, J) + (numer, denom), _ = multi_hilbert_series(AA; algorithm, parent) + return (numer, denom), (H, iso) + end + + # Now we may assume that the grading group is free Abelian. + m = ngens(G) + n = ngens(R) + HSRing = _hilbert_series_ring(parent, m) + + # Get the weights as Int values: W[k] contain the weight(s) of x[k] + W = [[ Int(R.d[i][j]) for j in 1:m] for i in 1:n] + fac_denom = _hilbert_series_denominator(HSRing, W) + + # Old method below without factorization; left for debugging + # q = one(parent) + # for i = 1:n + # e = [Int(MI[i, :][j]) for j = 1:m] + # B = MPolyBuildCtx(parent) + # push_term!(B, 1, e) + # q = q*(1-finish(B)) + # end + # @assert evaluate(fac_denom) == q + + # Shortcut for the trivial case + iszero(I) && return (one(HSRing), fac_denom), (G, identity_map(G)) + + # In general refer to internal methods for monomial ideals + # TODO: Shouldn't the ordering be adapted to the grading in some sense? + numer = one(HSRing) + if backend == :Zach + LI = leading_ideal(I; ordering=degrevlex(gens(R))) # ??? better not to specify the grading ??? + numer = _numerator_monomial_multi_hilbert_series(LI, HSRing, m; algorithm=algorithm) + elseif backend == :Abbott + # TODO: Pass on the `algorithm` keyword argument also here. + numer = HSNum_abbott(A, HSRing) + else + error("backend ($(backend)) not found") + end + return (numer, fac_denom), (G, identity_map(G)) end + ### TODO: original version of multi_hilbert_series based on moving things to the positive orthant #function multi_hilbert_series(A::MPolyQuoRing) @@ -655,6 +690,7 @@ julia> A, _ = quo(R, I); julia> H = multi_hilbert_series_reduced(A); + julia> H[1][1] -t[1]^5*t[2]^-1 + t[1]^3 + t[1]^3*t[2]^-3 + t[1]^2 + t[1]^2*t[2]^-1 + t[1]^2*t[2]^-2 + t[1] + t[1]*t[2]^-1 + 1 @@ -705,7 +741,7 @@ G """ function multi_hilbert_series_reduced(A::MPolyQuoRing; algorithm::Symbol=:BayerStillmanA) (p, q), (H, iso) = multi_hilbert_series(A, algorithm=algorithm) - f = p//q + f = p//evaluate(q) p = numerator(f) q = denominator(f) sig = coeff(q.mpoly, -1 .* q.mindegs) diff --git a/src/Rings/mpoly-graded.jl b/src/Rings/mpoly-graded.jl index bd83f6b97fa4..c9406eb8760a 100644 --- a/src/Rings/mpoly-graded.jl +++ b/src/Rings/mpoly-graded.jl @@ -1543,435 +1543,6 @@ function Base.show(io::IO, h::HilbertData) print(io, "Hilbert Series for $(h.I), data: $(h.data), weights: $(h.weights)") ###new end -######################################################################## -# Compute the Hilbert series from the monomial leading ideal using -# different bisection strategies along the lines of -# [Kreuzer, Robbiano: Computational Commutative Algebra 2, Springer] -# Section 5.3. -######################################################################## - -_divides(a::Vector{Int}, b::Vector{Int}) = all(k->(a[k]>=b[k]), 1:length(a)) - -function _gcd(a::Vector{Int}, b::Vector{Int}) - return [a[k] > b[k] ? b[k] : a[k] for k in 1:length(a)] -end - -function _are_pairwise_coprime(a::Vector{Vector{Int}}) - length(a) <= 1 && return true - n = length(a) - m = length(a[1]) - for i in 1:n-1 - for j in i+1:n - all(k -> a[i][k] == 0 || a[j][k] == 0, 1:m) || return false - end - end - return true -end - -### Assume that a is minimal. We divide by `pivot` and return -# a minimal set of generators for the resulting monomial ideal. -function _divide_by(a::Vector{Vector{Int}}, pivot::Vector{Int}) - # The good ones will contribute to a minimal generating - # set of the lhs ideal. - # - # The bad monomials come from those which hop over the boundaries of - # the monomial diagram by the shift. Their span has a new - # generator which is collected in `bad` a priori. It is checked - # whether they become superfluous and if not, they are added to - # the good ones. - good = sizehint!(Vector{Vector{Int}}(), length(a)) - bad = sizehint!(Vector{Vector{Int}}(), length(a)) - for e in a - if _divides(e, pivot) - push!(good, e - pivot) - else - push!(bad, e) - end - end - - # pre-allocate m so that don't need to allocate it again each loop iteration - m = similar(pivot) - for e in bad - # the next line computers m = [k < 0 ? 0 : k for k in e] - # but without allocations - for i in 1:length(e) - m[i] = max(e[i] - pivot[i], 0) - end - - # check whether the new monomial m is already in the span - # of the good ones. If yes, discard it. If not, discard those - # elements of the good ones that are in the span of m and put - # m in the list of good ones instead. - if all(x->!_divides(m, x), good) - # Remove those 'good' elements which are multiples of m - filter!(x -> !_divides(x, m), good) - push!(good, copy(m)) - end - end - - return good -end - -### This implements the speedup from the CoCoA strategy, -# see Exercise 4 in Section 5.3 -function _divide_by_monomial_power(a::Vector{Vector{Int}}, j::Int, k::Int) - divides(a::Vector{Int}, b::Vector{Int}) = all(k->(k>=0), b[1:j-1] - a[1:j-1]) && all(k->(k>=0), b[j+1:end] - a[j+1:end]) - kept = [e for e in a if e[j] == 0] - for i in 1:k - next_slice = [e for e in a if e[j] == i] - still_kept = next_slice - for e in kept - all(v->!divides(v, e), next_slice) && push!(still_kept, e) - end - kept = still_kept - end - still_kept = vcat(still_kept, [e for e in a if e[j] > k]) - - result = Vector{Vector{Int}}() - for e in still_kept - v = copy(e) - v[j] = (e[j] > k ? e[j] - k : 0) - push!(result, v) - end - return result -end - -function _find_maximum(a::Vector{Int}) - m = a[1] - j = 1 - for i in 1:length(a) - if a[i] > m - m = a[i] - j = i - end - end - return j -end - -####################################################################### -# 06.10.2022 -# The following internal routine can probably still be tuned. -# -# Compared to the implementation in Singular, it is still about 10 -# times slower. We spotted the following possible deficits: -# -# 1) The returned polynomials are not yet built with MPolyBuildCtx. -# We tried that, but as for now, the code for the build context -# is even slower than the direct implementation that is in place -# now. Once MPolyBuildCtx is tuned, we should try that again; -# the code snippets are still there. In total, building the return -# values takes up more than one third of the computation time -# at this moment. -# -# 2) Numerous allocations for integer vectors. The code below -# performs lots of iterative allocations for lists of integer -# vectors. If we constructed a container data structure to -# maintain this list internally and do the allocations at once -# (for instance in a big matrix), this could significantly -# speed up the code. However, that is too much work for -# the time being, as long as a high-performing Hilbert series -# computation is not of technical importance. -# -# 3) Singular uses bitmasking for exponent vectors to decide on -# divisibility more quickly. This is particularly important in -# the method _divide_by. For singular, this leads to a speedup -# of factor 5. Here, only less than one third of the time is -# spent in `_divide_by`, so it does not seem overly important -# at this point, but might become relevant in the future. -# A particular modification of the singular version of bitmasking -# is to compute means for the exponents occurring for each variable -# and set the bits depending on whether a given exponent is greater -# or less than that mean value. - -function _hilbert_numerator_from_leading_exponents( - a::Vector{Vector{Int}}, - weight_matrix::Matrix{Int}, - return_ring::Ring, - #algorithm=:generator - #algorithm=:custom - #algorithm=:gcd - #algorithm=:indeterminate - #algorithm=:cocoa - algorithm::Symbol # =:BayerStillmanA, # This is by far the fastest strategy. Should be used. - # short exponent vectors where the k-th bit indicates that the k-th - # exponent is non-zero. - ) - t = gens(return_ring) - ret = _hilbert_numerator_trivial_cases(a, weight_matrix, return_ring, t) - ret !== nothing && return ret - - if algorithm == :BayerStillmanA - return _hilbert_numerator_bayer_stillman(a, weight_matrix, return_ring, t) - elseif algorithm == :custom - return _hilbert_numerator_custom(a, weight_matrix, return_ring, t) - elseif algorithm == :gcd # see Remark 5.3.11 - return _hilbert_numerator_gcd(a, weight_matrix, return_ring, t) - elseif algorithm == :generator # just choosing on random generator, cf. Remark 5.3.8 - return _hilbert_numerator_generator(a, weight_matrix, return_ring, t) - elseif algorithm == :indeterminate # see Remark 5.3.8 - return _hilbert_numerator_indeterminate(a, weight_matrix, return_ring, t) - elseif algorithm == :cocoa # see Remark 5.3.14 - return _hilbert_numerator_cocoa(a, weight_matrix, return_ring, t) - end - error("invalid algorithm") -end - -# compute t ^ (weight_matrix * expvec), where t == gens(S) -function _expvec_to_poly(S::Ring, t::Vector, weight_matrix::Matrix{Int}, expvec::Vector{Int}) - o = one(coefficient_ring(S)) # TODO: cache this ?! - return S([o], [weight_matrix * expvec]) -end - -# special case for univariate polynomial ring -function _expvec_to_poly(S::PolyRing, t::Vector, weight_matrix::Matrix{Int}, expvec::Vector{Int}) - @assert length(t) == 1 - @assert size(weight_matrix) == (1, length(expvec)) - - # compute the dot-product of weight_matrix[1,:] and expvec, but faster than dot - # TODO: what about overflows? - s = weight_matrix[1,1] * expvec[1] - for i in 2:length(expvec) - @inbounds s += weight_matrix[1,i] * expvec[i] - end - return t[1]^s -end - -function _hilbert_numerator_trivial_cases( - a::Vector{Vector{Int}}, weight_matrix::Matrix{Int}, - S::Ring, t::Vector, oneS = one(S) - ) - length(a) == 0 && return oneS - - # See Proposition 5.3.6 - if _are_pairwise_coprime(a) - return prod(oneS - _expvec_to_poly(S, t, weight_matrix, e) for e in a) - end - - return nothing -end - - -function _hilbert_numerator_cocoa( - a::Vector{Vector{Int}}, weight_matrix::Matrix{Int}, - S::Ring, t::Vector - ) - ret = _hilbert_numerator_trivial_cases(a, weight_matrix, S, t) - ret !== nothing && return ret - - n = length(a) - m = length(a[1]) - counters = [0 for i in 1:m] - for e in a - counters += [iszero(k) ? 0 : 1 for k in e] - end - j = _find_maximum(counters) - - p = a[rand(1:n)] - while p[j] == 0 - p = a[rand(1:n)] - end - - q = a[rand(1:n)] - while q[j] == 0 || p == q - q = a[rand(1:n)] - end - - pivot = [0 for i in 1:m] - pivot[j] = minimum([p[j], q[j]]) - - - ### Assembly of the quotient ideal with less generators - rhs = [e for e in a if !_divides(e, pivot)] - push!(rhs, pivot) - - ### Assembly of the division ideal with less total degree - lhs = _divide_by_monomial_power(a, j, pivot[j]) - - f = one(S) - for i in 1:nvars(S) - z = t[i] - f *= z^(sum([pivot[j]*weight_matrix[i, j] for j in 1:length(pivot)])) - end - - return _hilbert_numerator_cocoa(rhs, weight_matrix, S, t) + f*_hilbert_numerator_cocoa(lhs, weight_matrix, S, t) -end - -function _hilbert_numerator_indeterminate( - a::Vector{Vector{Int}}, weight_matrix::Matrix{Int}, - S::Ring, t::Vector - ) - ret = _hilbert_numerator_trivial_cases(a, weight_matrix, S, t) - ret !== nothing && return ret - - e = first(a) - found_at = findfirst(!iszero, e)::Int64 - pivot = zero(e) - pivot[found_at] = 1 - - ### Assembly of the quotient ideal with less generators - rhs = [e for e in a if e[found_at] == 0] - push!(rhs, pivot) - - ### Assembly of the division ideal with less total degree - lhs = _divide_by(a, pivot) - - f = one(S) - for i in 1:nvars(S) - z = t[i] - f *= z^(sum([pivot[j]*weight_matrix[i, j] for j in 1:length(pivot)])) - end - - return _hilbert_numerator_indeterminate(rhs, weight_matrix, S, t) + f*_hilbert_numerator_indeterminate(lhs, weight_matrix, S, t) -end - -function _hilbert_numerator_generator( - a::Vector{Vector{Int}}, weight_matrix::Matrix{Int}, - S::Ring, t::Vector - ) - ret = _hilbert_numerator_trivial_cases(a, weight_matrix, S, t) - ret !== nothing && return ret - - b = copy(a) - pivot = pop!(b) - - f = one(S) - for i in 1:nvars(S) - z = t[i] - f *= z^(sum([pivot[j]*weight_matrix[i, j] for j in 1:length(pivot)])) - end - - c = _divide_by(b, pivot) - p1 = _hilbert_numerator_generator(b, weight_matrix, S, t) - p2 = _hilbert_numerator_generator(c, weight_matrix, S, t) - - return p1 - f * p2 -end - -function _hilbert_numerator_gcd( - a::Vector{Vector{Int}}, weight_matrix::Matrix{Int}, - S::Ring, t::Vector - ) - ret = _hilbert_numerator_trivial_cases(a, weight_matrix, S, t) - ret !== nothing && return ret - - n = length(a) - counters = [0 for i in 1:length(a[1])] - for e in a - counters += [iszero(k) ? 0 : 1 for k in e] - end - j = _find_maximum(counters) - - p = a[rand(1:n)] - while p[j] == 0 - p = a[rand(1:n)] - end - - q = a[rand(1:n)] - while q[j] == 0 || p == q - q = a[rand(1:n)] - end - - pivot = _gcd(p, q) - - ### Assembly of the quotient ideal with less generators - rhs = [e for e in a if !_divides(e, pivot)] - push!(rhs, pivot) - - ### Assembly of the division ideal with less total degree - lhs = _divide_by(a, pivot) - - f = one(S) - for i in 1:nvars(S) - z = t[i] - f *= z^(sum([pivot[j]*weight_matrix[i, j] for j in 1:length(pivot)])) - end - - return _hilbert_numerator_gcd(rhs, weight_matrix, S, t) + f*_hilbert_numerator_gcd(lhs, weight_matrix, S, t) -end - -function _hilbert_numerator_custom( - a::Vector{Vector{Int}}, weight_matrix::Matrix{Int}, - S::Ring, t::Vector - ) - ret = _hilbert_numerator_trivial_cases(a, weight_matrix, S, t) - ret !== nothing && return ret - - p = Vector{Int}() - q = Vector{Int}() - max_deg = 0 - for i in 1:length(a) - b = a[i] - for j in i+1:length(a) - c = a[j] - r = _gcd(b, c) - if sum(r) > max_deg - max_deg = sum(r) - p = b - q = c - end - end - end - - ### Assembly of the quotient ideal with less generators - pivot = _gcd(p, q) - rhs = [e for e in a if !_divides(e, pivot)] - push!(rhs, pivot) - - ### Assembly of the division ideal with less total degree - lhs = _divide_by(a, pivot) - - f = one(S) - for i in 1:nvars(S) - z = t[i] - f *= z^(sum([pivot[j]*weight_matrix[i, j] for j in 1:length(pivot)])) - end - - return _hilbert_numerator_custom(rhs, weight_matrix, S, t) + f*_hilbert_numerator_custom(lhs, weight_matrix, S, t) -end - -function _hilbert_numerator_bayer_stillman( - a::Vector{Vector{Int}}, weight_matrix::Matrix{Int}, - S::Ring, t::Vector, oneS = one(S) - ) - ########################################################################### - # For this strategy see - # - # Bayer, Stillman: Computation of Hilber Series - # J. Symbolic Computation (1992) No. 14, pp. 31--50 - # - # Algorithm 2.6, page 35 - ########################################################################### - ret = _hilbert_numerator_trivial_cases(a, weight_matrix, S, t, oneS) - ret !== nothing && return ret - - # make sure we have lexicographically ordered monomials - sort!(a; alg=QuickSort) - - # initialize the result - h = oneS - _expvec_to_poly(S, t, weight_matrix, a[1]) - linear_mons = Vector{Int}() - for i in 2:length(a) - J = _divide_by(a[1:i-1], a[i]) - empty!(linear_mons) - J1 = filter!(J) do m - k = findfirst(!iszero, m) - k === nothing && return false # filter out zero vector - if m[k] == 1 && findnext(!iszero, m, k + 1) === nothing - push!(linear_mons, k) - return false - end - return true - end - q = _hilbert_numerator_bayer_stillman(J1, weight_matrix, S, t, oneS) - for k in linear_mons - q *= (oneS - prod(t[i]^weight_matrix[i, k] for i in 1:length(t))) - end - - h -= q * _expvec_to_poly(S, t, weight_matrix, a[i]) - end - return h -end - ############################################################################ ### Homogenization and Dehomogenization ############################################################################ diff --git a/test/Modules/hilbert.jl b/test/Modules/hilbert.jl new file mode 100644 index 000000000000..839744ca12f7 --- /dev/null +++ b/test/Modules/hilbert.jl @@ -0,0 +1,25 @@ +@testset "Hilbert series and free resolution" begin + R, _ = polynomial_ring(QQ, ["x", "y", "z"]); + Z = abelian_group(0); + Rg, (x, y, z) = grade(R, [3*Z[1],Z[1],-Z[1]]); + F = graded_free_module(Rg, 1); + A = Rg[x; y]; + B = Rg[x^2+y^6; y^7; z^4]; + M = SubquoModule(F, A, B); + fr = free_resolution(M) + + # There is no length(fr), so instead we do the following + indexes = range(fr.C); + fr_len = first(indexes) + + num,_ = hilbert_series(fr[fr_len]) + for i in fr_len:-1:1 + phi = map(fr,i) + N = cokernel(phi) + numer,denom = hilbert_series(N) + numer_next,_ = hilbert_series(fr[i-1]) + @test numer == numer_next- num + num = numer + end + +end diff --git a/test/Modules/runtests.jl b/test/Modules/runtests.jl index 3169bee901cf..2584fb970b49 100644 --- a/test/Modules/runtests.jl +++ b/test/Modules/runtests.jl @@ -4,6 +4,7 @@ using Test include("UngradedModules.jl") include("FreeModElem-orderings.jl") include("ModulesGraded.jl") +include("hilbert.jl") include("module-localizations.jl") include("local_rings.jl") include("MPolyQuo.jl") diff --git a/test/Rings/hilbert.jl b/test/Rings/hilbert.jl new file mode 100644 index 000000000000..0dce2c506af1 --- /dev/null +++ b/test/Rings/hilbert.jl @@ -0,0 +1,115 @@ +@testset "Hilbert series" begin + P, x = graded_polynomial_ring(QQ, 5, "x", [1 1 1 1 1; 1 -2 3 -4 5]); + G = + [ + x[1]^30*x[2]^16*x[3]^7*x[4]^72*x[5]^31, + x[1]^11*x[2]^20*x[3]^29*x[4]^93*x[5]^5, + x[1]^5*x[2]^9*x[3]^43*x[4]^7*x[5]^97, + x[1]^57*x[2]^6*x[3]^7*x[4]^56*x[5]^37, + x[1]^22*x[2]^40*x[3]^27*x[4]^67*x[5]^18, + x[1]^5*x[2]^5*x[3]^7*x[4]^78*x[5]^81, + x[1]^85*x[2]^15*x[3]^12*x[4]^41*x[5]^37, + x[1]^13*x[2]^35*x[3]^2*x[4]^64*x[5]^84, + x[1]^31*x[2]^57*x[3]^12*x[4]^93*x[5]^6, + x[1]^37*x[2]^59*x[3]^88*x[4]^16*x[5]^7, + x[1]^47*x[2]^63*x[3]^32*x[4]^28*x[5]^41, + x[1]^18*x[2]^19*x[3]^95*x[4]^7*x[5]^73, + x[1]^37*x[2]^84*x[3]^11*x[4]^48*x[5]^32, + x[1]^85*x[2]*x[3]^73*x[4]^24*x[5]^31, + x[1]^59*x[2]^56*x[3]^41*x[4]^12*x[5]^50, + x[1]^58*x[2]^36*x[3]*x[4]^35*x[5]^89, + x[1]^71*x[2]^12*x[3]^36*x[4]^76*x[5]^25, + x[1]^58*x[2]^32*x[3]^85*x[4]^44*x[5], + x[1]^61*x[2]^81*x[3]^15*x[4]^59*x[5]^8, + x[1]^84*x[2]^49*x[3]^32*x[4]^52*x[5]^7, + x[1]^40*x[2]^3*x[3]^54*x[4]^39*x[5]^89, + x[1]^12*x[2]^52*x[3]^100*x[4]^49*x[5]^14, + x[1]^75*x[2]^78*x[3]^34*x[4]*x[5]^41, + x[1]^46*x[2]^22*x[3]^99*x[4]^49*x[5]^14, + x[1]^95*x[2]^11*x[3]^63*x[4]^52*x[5]^10, + x[1]^32*x[2]^47*x[3]^97*x[4]^32*x[5]^26, + x[1]^52*x[2]^64*x[3]^62*x[4]^13*x[5]^47, + x[1]^70*x[2]^46*x[3]^90*x[4]^13*x[5]^21, + x[1]^3*x[2]^67*x[3]^90*x[4]^45*x[5]^52, + x[1]^17*x[2]^100*x[3]^58*x[4]^62*x[5]^21, + x[1]^45*x[2]^54*x[3]^65*x[4]^64*x[5]^32, + x[1]^24*x[2]^85*x[3]^27*x[4]^49*x[5]^76, + x[1]^28*x[2]^83*x[3]^9*x[4]^97*x[5]^45, + x[1]^40*x[2]^52*x[3]^99*x[4]^27*x[5]^49, + x[1]^88*x[2]^36*x[3]^26*x[4]^30*x[5]^90, + x[1]^36*x[2]^68*x[3]^76*x[4]^81*x[5]^9, + x[1]^48*x[2]^25*x[3]^80*x[4]^40*x[5]^83, + x[1]^93*x[2]^14*x[3]^80*x[4]^89*x[5]^3, + x[1]^5*x[2]^61*x[3]^65*x[4]^65*x[5]^94, + x[1]^48*x[2]^91*x[3]^70*x[4]^66*x[5]^16, + x[1]^39*x[2]^79*x[3]^98*x[4]^2*x[5]^75, + x[1]^4*x[2]^60*x[3]^74*x[4]^56*x[5]^100, + x[1]^59*x[2]^100*x[3]^74*x[4]^51*x[5]^12, + x[1]^90*x[2]^61*x[3]^85*x[4]^43*x[5]^19, + x[1]^44*x[2]^97*x[3]^39*x[4]^27*x[5]^97, + x[1]^91*x[2]^37*x[3]^2*x[4]^97*x[5]^80, + x[1]^67*x[2]^44*x[3]^99*x[4]^5*x[5]^93, + x[1]^95*x[2]^93*x[3]^88*x[4]^30*x[5]^2, + x[1]^91*x[3]^84*x[4]^59*x[5]^76, + x[1]^62*x[2]^3*x[3]^96*x[4]^72*x[5]^84 + ]; + + I = ideal(P,G); + PmodI, _ = quo(P,I); + (num1,_), (_,_) = multi_hilbert_series(PmodI); + S = parent(num1) + (num2,_), (_,_) = multi_hilbert_series(PmodI; parent=S, backend=:Zach); + @test num2 == num1 + + P2, x = graded_polynomial_ring(QQ, 5, "x", [1 1 1 1 1]) + G = P2.(G); + + I = ideal(P2,G); + PmodI, _ = quo(P2,I); + num1, _ = hilbert_series(PmodI); + L, q = LaurentPolynomialRing(ZZ, "q"); + num2, _ = hilbert_series(PmodI; parent=L); + @test num2 == evaluate(num1, q) +end + +@testset "second round of Hilbert series" begin + # Standard graded + R, (x, y, z) = polynomial_ring(QQ, ["x", "y", "z"]); + J = ideal(R, [y-x^2, z-x^3]); + I = homogenization(J, "w"); + A, _ = quo(base_ring(I), I); + numer1, denom1 = hilbert_series(A); + S, t = LaurentPolynomialRing(ZZ, "t"); + numer2, denom2 = hilbert_series(A; parent=S); + Smult, (T,) = polynomial_ring(ZZ, ["t"]); + numer3, denom3 = hilbert_series(A; parent=Smult); + @test numer3 == evaluate(numer1, T) + @test numer2 == evaluate(numer1, t) + + # Negative grading + RR, (X, Y) = graded_polynomial_ring(QQ, ["X", "Y"], [-1, -1]) + JJ = ideal(RR, X^2 - Y^2); + A, _ = quo(base_ring(JJ), JJ); + (numer1, denom1), _ = multi_hilbert_series(A); + S, t = LaurentPolynomialRing(ZZ, "t"); + (numer2, denom2), _ = multi_hilbert_series(A; parent=S); + Smult, (T,) = polynomial_ring(ZZ, ["t"]); + @test_throws DomainError multi_hilbert_series(A; parent=Smult) + + # Graded by commutative group + G = free_abelian_group(2); + G, _ = quo(G, [G[1]-3*G[2]]); + RR, (X, Y) = graded_polynomial_ring(QQ, ["X", "Y"], [G[1], G[2]]); + JJ = ideal(RR, X^2 - Y^6); + A, _ = quo(base_ring(JJ), JJ); + (numer1, denom1), (H, iso) = multi_hilbert_series(A); + @test is_free(H) && isone(rank(H)) + S, t = LaurentPolynomialRing(ZZ, "t"); + (numer2, denom2), (H, iso) = multi_hilbert_series(A; parent=S); + @test is_free(H) && isone(rank(H)) + Smult, (T,) = polynomial_ring(ZZ, ["t"]); + (numer3, denom3), (H, iso) = multi_hilbert_series(A; parent=Smult); + @test is_free(H) && isone(rank(H)) + @test numer1 == evaluate(numer2, T) + @test numer3 == evaluate(numer2, first(gens(parent(numer3)))) +end diff --git a/test/Rings/runtests.jl b/test/Rings/runtests.jl index 1cf01fd9a6f9..be400c544920 100644 --- a/test/Rings/runtests.jl +++ b/test/Rings/runtests.jl @@ -30,3 +30,4 @@ include("PBWAlgebraQuo.jl") include("FreeAssAlgIdeal.jl") include("binomial-ideals.jl") +include("hilbert.jl")