# implementation of the Staggered Linear Basis Algorithm of Gebauer and M\"oller
# "Buchberger's algorithm and Staggered Linear Bases"
# in Proceedings of SYMSAC 1986 (Waterloo/Ontario), pages 218--221
# ACM Press, 1986

# Originally from http://www.math.usm.edu/perry/Research/staggered_linear_basis.py
# small changes for JAS compatibility, changes are labled with JAS
# $Id: staggered_linear_basis.py 4986 2014-11-17 23:03:07Z kredel $

# usage: staglinbasis(F) where F is a list of polynomials
# result: groebner basis of F
# example:
#   sage: R = PolynomialRing(GF(32003),x,5)
#   sage: I = sage.rings.ideal.Cyclic(R,5).homogenize()
#   sage: F = list(I.gens())
#   sage: B = staglinbasis(F)
#   >> reduction to zero: (9, 12, x1*x2*x3^3*x4^2) x2
#   >> reduction to zero: (7, 12, x1*x2*x3^3*x4^2) x2
#   >> reduction to zero: (5, 12, x1^2*x3^3*x4^2) x1
#   ... [more omitted]
#   >> 32 reductions to zero
#   sage: len(B)
#   >> 42

# minimal_pair
# used to select a minimal critical pair from a list;
# form of the list is (i,j,t) where t=lcm(lt(g[i]),lt(g[j]));
# chosen by smallest lcm
# inputs: list of pairs P
# outputs: minimal pair in P
def minimal_pair(P):
  result = P.pop()
  P.add(result)
  for p in P:
    if p[2] < result[2]:
      result = p
  return result

# spol
# compute the s-polynomial corresponding to the critical pair p
# inputs: critical pair p, basis G
# outputs: s-polynomial
def spol(p,G):
  i,j,tij = p
  R = G[i].parent() # JAS
  ti = R.monomial_quotient(tij,G[i].lm())
  tj = R.monomial_quotient(tij,G[j].lm())
  return ti*G[i] - tj*G[j]

# monid_quotient
# computes the quotient I:J of the monomial ideals I and J
# for each tI in I, and for each tJ in J, it computes lcm(tI,tJ)/tJ
# this corresponds to computing the u's such that uv is in I for each v in J
# note that the ideals must be simplified for this to work correctly
# (see below)
# inputs: lists of monomials I and J (generators of ideals)
# outputs: list of monomials that generate I:J
def monid_quotient(I,J,R): # JAS
  result = set()
  for tI in I:
    # R = tI.parent() #JAS
    for tJ in J:
      u = R.monomial_quotient(tI.lcm(tJ),tJ).lm() # JAS
      if not any(R.monomial_divides(t,u) for t in result):
        result.add(u)
  return result

# monid_add
# expands the monomial ideal I by the monomial t
# simplifies the result
# inputs: list of monomials I (generators of ideal), monomial t
# outputs: simplified I+(t)
def monid_add(I,t):
  result = set()
  R = t.parent()
  for u in I:
    if not(R.monomial_divides(t,u)) and not(R.monomial_divides(u,t)):
      result.add(u)
  if not any(R.monomial_divides(u,t) for u in result):
    result.add(t)
  return result

# simplify_monid
# simplifies the list of generators of the monomial ideal I
# by pruning redundant elements
# inputs: list of monomials I (generators of ideal)
# outputs: simplified I
def simplify_monid(I,R): # JAS
  result = set()
  for t in I:
    #R = t.parent() #JAS
    if not any(u != t and R.monomial_divides(u,t) for u in I):
      result.add(t)
  return result

# staglinbasis
# computes the Groebner Basis of a polynomial ideal using the
# Staggered Linear Basis algorithm of Gebauer-Moeller
# inputs: list of polynomials F
# outputs: Groebner basis of F
def staglinbasis(F):
  # tracks the number of zero reductions
  num_zrs = 0
  # get polynomial ring
  R = F[0].parent()
  G = list(F) # JAS
  m = len(G)
  # Z is the set of monomial used to detect principal syzygies
  # in general, Z should be simplified anytime you modify it
  Z = [set([G[j].lm() for j in xrange(i)]) for i in xrange(m)]
  Z = [simplify_monid(ZI,R) for ZI in Z] # JAS
  # initial set of critical pairs
  P = set()
  for i in xrange(m-1):
    for j in xrange(i+1,m):
      tij = G[i].lm().lcm(G[j].lm())
      tj = R.monomial_quotient(tij,G[j].lm()) # JAS
      if not any(R.monomial_divides(t,tj) for t in Z[j]):
        P.add((i,j,tij))
  # main loop
  while len(P) != 0:
    # choose the pair w/smallest lcm
    p = minimal_pair(P)
    P.remove(p)
    # compute S-polynomial, reduce it
    s = spol(p,G)
    i,j,tij = p
    r = s.reduce(G)
    if r == 0:
      print "reduction to zero:", p, R.monomial_quotient(tij,G[j].lm())
      num_zrs += 1
    else:
      # new polynomial; add to basis
      r = r.monic() # JAS #r *= r.lc()**(-1)
      G.append(r)
      # update ideals here, also below
      # need new ideal for G[-1]=r, essentially
      #   Znew = (Z[j] + t[i])):(lcm(t[i],t[j])/t[j]) + (t[1], ..., t[m])
      # where t[i] = lm(g[i])
      Znew = Z[j].copy() # JAS
      Znew.update([G[i].lm()])
      Znew = simplify_monid(Znew,R) #JAS
      tj = R.monomial_quotient(tij,G[j].lm()).lm() # JAS
      Znew = monid_quotient(Znew,[tj],R) # JAS
      Znew.update([g.lm() for g in G[:-1]])
      Znew = simplify_monid(Znew,R)
      Z.append(Znew)
      # compute new critical pairs w/new polynomial
      m = len(G) - 1
      for i in xrange(m):
        tim = G[i].lm().lcm(G[m].lm())
        tm = R.monomial_quotient(tim,G[m].lm())
        # don't add pairs detected by Z[m]
        if not any(R.monomial_divides(t,tm) for t in Z[m]):
          P.add((i,m,tim))
      # update Z[j]
      Z[j].add(R.monomial_quotient(tij,G[j].lm()).lm()) # JAS
      Z[j] = simplify_monid(Z[j],R) # JAS
      # prune pairs detected by new Z[j]
      Q = set()
      for (i,j,tij) in P:
        tj = R.monomial_quotient(tij,G[j].lm())
        if any(R.monomial_divides(t,tj) for t in Z[j]):
          Q.add((i,j,tij))
      P.difference_update(Q)
  print num_zrs, "reductions to zero"
  return G