import Base./, Base.%


function S_poly(A::Multinomial, B::Multinomial)
    x = A[end].indeterminate
    tmp = [max(A[end].exponent[i], B[end].exponent[i]) for i in 1:x.degree]
    return B[end].coefficient*x^(tmp-A[end].exponent)*A -
        A[end].coefficient*x^(tmp-B[end].exponent)*B
end

"""
Does monomial a divide monomial b?
"""
function divides(a::Monomial, b::Monomial)
   !isapprox(a.coefficient, 0) &&
       all(a.exponent .<= b.exponent) &&
       a.indeterminate == b.indeterminate
end

"""
Quotient of b and a. Assumes divisibility
"""
function /(b::Monomial, a::Monomial)
    Monomial(b.indeterminate,
             b.exponent-a.exponent,
             b.coefficient/a.coefficient)
end

""" function %(A::Multinomial, G::Array{Multinomial,1})

Divide A by the array of multinomials, G returning the remainder.
"""
function %(A::Multinomial, G::Array{Multinomial,1})
    h = deepcopy(A)
    i = 1
    while i <= length(G)
        g = G[i]
        # find first j, if any, such that g[end] divides h[end-j]
        j = 0
        while j < length(h)
            if divides(g[end], h[end-j])
                break
            else
                j += 1
            end
        end
        if j < length(h)
            h = h - (h[end-j]/g[end])*g
            i = 1
        else
            i += 1
        end
    end
    return h
end

""" function iszero(A::Multinomial)

Is multinomial A zero?
"""
function iszero(A::Multinomial)
    length(A)==0 || (length(A)==1 && isapprox(A[1].coefficient, 0))
end

"""
function buchberger(F::Array{Multinomial,1})

Buchberger's algorithm for finding a Groebner basis for the ideal generated by the system of multinomials, F.

This is a straightforward, unoptimized implementation of the algorithm as described by Buchberger himself in [Gröbner Bases: A Short Introduction for Systems Theorists](http://link.springer.com/chapter/10.1007/3-540-45654-6_1)[PDF](http://people.reed.edu/~davidp/pcmi/buchberger.pdf) as follows:

Start with G = F .
For any pair of polynomials f1 , f2 in G:
  Compute the S✁polynomial of f1, f2
  and reduce it to a reduced form h w.r.t. G.
  If h = 0, consider the next pair.
  If h is not 0, add h to G and iterate.

Note that G never shrinks. In the worst case it can become intractably large. Moreover, the final G may contain many redundant multinomials. Neither of these problems aare faced here.
"""
function buchberger(F::Array{Multinomial,1};
                    increment::Int=1000,
                    maxiterations::Int=1000,
                    # for re-entry:
                    p1::Int=1,
                    p2::Int=2)
    ## Initialization ##
    # Expandable storage for Groebner basis
    nG = max(length(F), increment)
    iG = length(F)
    G = Array(Multinomial, nG)
    # Initialize G = F
    for i in 1:iG G[i] = deepcopy(F[i]) end
    ## Main loop ##
    iteration = 1
    # Continue while new pairs are available and max iterations
    # are not exceeded.
    while p2 <= iG && iteration <= maxiterations
        f1 = G[p1]
        f2 = G[p2]
        # Compute their S-polynomial
        S = S_poly(f1,f2)
        # Find S mod G
        h = S%G[1:iG]
        if !iszero(h)
            iG += 1
            if iG > length(G)
                G = vcat(G,Array(Monomial, increment))
            end
            G[iG] = h
        end
        # bookkeeping
        iteration += 1
        # next pair: note p2 increments more slowly than p1
        # so that all pairs p1 < p2 <= k are covered before
        # pairs p1 < p2 = k+1 are reached.
        p1, p2 = p1+1 < p2 ? (p1+1,p2) : (1, p2+1)
    end
    # Return G, p1, p2
    return G, iG, p1, p2
end

""" normalize!(G!::Array{Multinomial,1})

Convert all leading coefficients in  Groebner basis, G, to 1 modififying G in place. 
"""
function normalize!(G::Array{Multinomial,1})
    for g in G
        d = g[end].coefficient
        for m in g
            m.coefficient /= d
        end
    end
    nothing
end

function mark_nonredundant(G::Array{Multinomial,1})
    marks = trues(length(G))
    for i in eachindex(G)
        m = G[i][end]
        for j in eachindex(G)
            if j==i continue end
            if divides(G[j][end], m)
                marks[i] = false
                break
            end
        end
    end
    return marks
end

function reduce_Groebner(G::Array{Multinomial,1})
    Gtmp = deepcopy(G)
    normalize!(Gtmp)
    Gtmp = Gtmp[mark_nonredundant(Gtmp)]
    idx = trues(length(Gtmp))
    for i in eachindex(Gtmp)
        idx[i] = false
        Gtmp[i] = Gtmp[i]%Gtmp[idx]
        idx[i] = true
    end
    return Gtmp
end

function demo_mod()
    x = Indeterminate(2)
    P = 1 + x^[1,2]
    println("Example 1\n")
    @show P
    println()
    G = Array(Multinomial,2)
    G[1]= 1 + x^[1,1]
    G[2]= 1 + x^[0,1]
    @show G[1]
    @show G[2]
    println()
    @show P%G
    println("\nExample 2\n")
    P = x^[1,2] + x^[2,1] + x^[0,2]
    @show P
    println()
    G[1] = x^[1,1]-1
    G[2] = x^[0,2]-1
    @show G[1]
    @show G[2]
    println()
    @show P%G
    nothing
end

function demo_alg()
    @show x = Indeterminate(3)
    println()
    @show A = 1 + 2*x^[1,2,3] + 3*x^[3,2,1]
    println()
    @show B = 2 + 5* x^[1,3,2]
    println()
    @show S_poly(A,B)
    F = Array(Multinomial,2)
    F[1] = A; F[2] = B
    println()
    @show F
    println()
    G, iG, p1, p2 = buchberger(F, maxiterations=10)
    @show iG, p1, p2
    println()
    @show G[iG]
    nothing
end

function demo_bb()
    @show x = Indeterminate(2)
    println()
    F = Array(Multinomial,2)   
    F[1] = x^[1,1] - x^[1,0]
    F[2] = x^[2,0] - x^[0,1]
    @show F[1]
    @show F[2]
    println("\nRunning buchberger on F. Result:\n")
    G, iG, p1, p2 = buchberger(F, maxiterations=50)
    @show iG, p1, p2
    println()
    for i in 1:iG @show G[i] end
    println("\nReducing the Groebner basis found above.\n")
    Gred = reduce_Groebner(G[1:iG])
    for i in eachindex(Gred) @show Gred[i] end
    nothing
end

function demo_reduce()
    @show x = Indeterminate(2)
    println()
    G = Array(Multinomial,5)   
    G[1] = x^[3,0] - 2*x^[1,1]
    G[2] = x^[2,1] - 2*x^[0,2] + x^[1,0]
    G[3] = [x^[2,0]]
    G[4] = [x^[1,1]]
    G[5] = x^[0,2] - (1//2)*x^[1,0]
    @show G[1]
    @show G[2]
    @show G[3]
    @show G[4]
    @show G[5]
    println("\nReducing the Groebner basis shown above.\n")
    Gred = reduce_Groebner(G[1:5])
    for i in eachindex(Gred) @show Gred[i] end
    nothing
end

function bug1()
    x=Indeterminate(2)
    @show -x^[0,2]
    println()
    @show x^[0,4]-x^[0,2]
    println()
    @show -(x^[0,4],x^[0,2])
    nothing
end

    
    
    
    


