import AbstractAlgebra import Singular zero_reductions = 0 non_zero_reductions = 0 prod_criteria = 0 # structure for pairset mutable struct Pairset len::Int64 pairs::Array{Tuple{Int64,Int64}, 1} end # criteria for useless pairs function product_criterion(pair::Tuple{Int64,Int64}, G::Array{Singular.spoly,1}) lm1 = Singular.lm(G[pair[1]]) lm2 = Singular.lm(G[pair[2]]) if Singular.lcm(lm1, lm2) == lm1 * lm2 return true end return false end # initialize pairset function init_pairset(idx::Int64) len = 0 for k in 1:idx-1 len += k end p = Array{Tuple{Int64,Int64},1}(undef, len) P = Pairset(len, p) k = 1 for i = 2:idx for j = 1:i-1 P.pairs[k] = (i,j) k = k+1 end end return P end # add new pairs to pairset function new_pairs(P::Pairset, idx::Int64) k = P.len + 1 P.len = P.len + idx - 1 resize!(P.pairs, P.len) for i in 1:idx-1 P.pairs[k] = (idx,i) k = k+1 end end # multiplier for reduction function multiplier_exp(h::Singular.spoly, g::Singular.spoly) x = Singular.gens(parent(h)) eh = Singular.lead_exponent(h) eg = Singular.lead_exponent(g) mul = one(parent(h)) for k in 1:length(eh) diff = eh[k] - eg[k] if diff < 0 return 0 end if diff > 0 mul *= x[k]^diff continue end end return mul end # multipliers for spolynomials function muls_exp(i::Int64, j::Int64, G::Array{Singular.spoly,1}) x = Singular.gens(parent(G[1])) ei = Singular.lead_exponent(G[i]) ej = Singular.lead_exponent(G[j]) mji = one(parent(G[1])) mij = one(parent(G[1])) for k in 1:length(ei) eij = ej[k] - ei[k] if eij > 0 mji *= x[k]^eij continue end if eij < 0 mij *= x[k]^(-eij) end end return mji, mij end function muls(i::Int64, j::Int64, G::Array{Singular.spoly,1}) lmi = Singular.lm(G[i]) lmj = Singular.lm(G[j]) lcm = Singular.lcm(lmi, lmj) mji = Singular.div(lcm, lmi) mij = Singular.div(lcm, lmj) return mji, mij end # generation of spolynomial of G[i] and G[j] function spoly(pair::Tuple{Int64,Int64}, G::Array{Singular.spoly,1}) i = pair[1] j = pair[2] mji, mij = muls(i, j, G) h = 1//Singular.lc(G[i]) * mji * G[i] - 1//Singular.lc(G[j]) * mij * G[j] return h end # compute normal form resp. standard expression of h w.r.t. G function normal_form(h::Singular.spoly, G::Array{Singular.spoly,1}) i = 0 while true if h == 0 global zero_reductions += 1 return 0 end i = 1 while i <= length(G) mul = Singular.div(Singular.lm(h), Singular.lm(G[i])) # mul = multiplier_exp(h, G[i]) if mul != 0 h = h - Singular.lc(h) * 1//Singular.lc(G[i]) * mul * G[i] i = 1 break end i = i + 1 end if i > length(G) global non_zero_reductions += 1 return h end end end function bba(id::Singular.sideal{T}) where T <: AbstractAlgebra.RingElem global zero_reductions = 0 global non_zero_reductions = 0 global prod_criteria = 0 # get number of generators in ideal l = Singular.ngens(I) R = I.base_ring # array for storing gb G = Array{Singular.spoly,1}(undef, l) # put all initial generators in G for i in 1:l G[i] = I[i] end # generate list of tuples of indices for s-polynomials P = init_pairset(l) # as long as the pairset is not empty while P.len != 0 # always take the last element and adjust length of P pair = P.pairs[P.len] P.len = P.len - 1 if product_criterion(pair, G) global prod_criteria += 1 continue end h = spoly(pair, G) h = normal_form(h, G) # if Buchberger's criterion is not fulfilled if h != 0 new_pairs(P, length(G)+1) push!(G, h) end end # return gb with setting isGB flag set GB = Singular.Ideal(R, G...) GB.isGB = true println("#zero reductions ", zero_reductions) println("#non-zero reductions ", non_zero_reductions) println("#product criteria ", prod_criteria) return GB end # R, (x,y,z) = Singular.PolynomialRing(Singular.QQ,["x","y","z"]); # I = Singular.Ideal(R, x^2+y*z, y^2*x-3*x, x*y*z-y^3);