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);