# implementation of Arri's algorithm to compute a Groebner basis
# "The F5 Criterion revised"
# paper submitted to JSC in 2009, accepted conditionally in 2010 subject to revision,
# resubmitted 2010
# http://arxiv.org/abs/1012.3664

# usage: f5_arri(F) where F is a list of polynomials
# result: groebner basis of F, with signatures
# example:
#   sage: R = PolynomialRing(GF(32003),x,5)
#   sage: I = sage.rings.ideal.Cyclic(R,5).homogenize()
#   sage: F = list(I.gens())
#   sage: B = f5_arri(F)
#   >> There were 0 reductions to zero.
#   sage: len(B)
#   >> 166
#   sage: F.reverse()
#   sage: B = f5_arri(F)
#   >> There were 0 reductions to zero.
#   sage: len(B)
#   >> 39
# as with all incremental algorithms, efficiency depends on the ordering of F

verbose_zero_reduction = 0

verbose_generation = 1
verbose_reduction = 1

verbose_rejection = 1

# polysort
# used to sort polynomials by increasing signature
# inputs: (f,S(f)) and (g,S(g)), labeled polynomials
# outputs: 1 if (f,S(f)) > (g,S(g))
#        -1 if (f,S(f)) < (g,S(g))
#         0 if (f,S(f)) = (g,S(g))
def polysort(f,g):
  fp,Sf = f
  gp,Sg = g
  sigord = sigsort(Sf,Sg)
  if sigord != 0:
    return sigord
  else:
    if fp.lm() > gp.lm():
      return 1
    else:
      return -1
  return 0

# sigsort
# used to sort signatures
# inputs: S(f) and S(g), labeled polynomials
# outputs: 1 if S(f) > S(g)
#        -1 if S(f) < S(g)
#         0 if S(f) = S(g)
def sigsort(s1,s2):
  t1,i1 = s1
  t2,i2 = s2
  if i1 < i2:
    return 1
  if i1 > i2:
    return -1
  if t1 > t2:
    return 1
  if t1 < t2:
    return -1
  return 0

# is_primitive
# tests whether a polynomial is primitive
# inputs: polynomial f, S(f), basis G
# outputs: whether uf*f is primitive
def is_primitive(f, Sf, G):
  R = f.parent()
  for (h,Sh) in G:
    if h != f and Sh[1] == Sf[1] and R.monomial_divides(Sh[0],Sf[0]):
      t = R.monomial_quotient(Sf[0],Sh[0])
      if t*h.lm() == f.lm():
        return False
  return True

# prune_rewritables
# removes rewritable elements from B
# this is NOT the same as Faugere's rewritable algorithm;
# see the paper for details
# inputs: polynomial ring R, wokring basis G, S-polynomials B
# outputs: (probably) smaller B
def prune_rewritables(R,G,B):
  Bnew = set()
  for (f,Sf) in B:
    rewritable = False
    t = f.lm()
    for (g,Sg) in B:
      if Sg[1] == Sf[1] and R.monomial_divides(Sg[0],Sf[0]):
        u = R.monomial_quotient(Sf[0],Sg[0])
        if u*g.lm() < t:
          rewritable = True
          continue
    if not rewritable:
      for (g,Sg) in G:
        if Sg[1] == Sf[1] and R.monomial_divides(Sg[0],Sf[0]):
          u = R.monomial_quotient(Sf[0],Sg[0])
          if u*g.lm() < t:
            rewritable = True
            continue
      if not rewritable:
        Bnew.add((f,Sf))
  return Bnew

# is_normal
# completes test whether a critical pair is a normal pair in Arri's definition
# (beyond inequality of signature)
# inputs: polynomials f and g, S(f), S(g) monomials uf and ug, basis G
# outputs: whether (f,g) is a normal pair
def is_normal(f, Sf, uf, G, B):
  R = f.parent()
  # Faugere's normal?
  if not any(Sh[1]>Sf[1] and R.monomial_divides(h.lm(),uf*Sf[0]) for (h,Sh) in G):
    # primitive?
    return is_primitive(f,Sf,G) and is_primitive(f,Sf,B)
    #return True
  # not even Faugere's normal, let alone primitive
  return False

# Sreduce
# inputs: polynomial f, S(f), basis G, rewritings RW
# outputs: S-irreducible form of f
def Sreduce(f,Sf,G,B):
  R = f.parent()
  reduced = True
  while reduced:
    reduced = False
    for (g,Sg) in G:
      if R.monomial_divides(g.lm(),f.lm()):
        t = R.monomial_quotient(f.lm(),g.lm())
        if sigsort(Sf,(t*Sg[0],Sg[1])) == 1:
          if is_normal(g, Sg, t, G, B) and is_primitive(g, Sg, G) and is_primitive(g, Sg, B):
            reduced = True
            f -= t*g
            if f == 0:
              return f
            f *= f.lc()**(-1)
          else:
            verbose("reduction of "+str(Sf)+" by "+str(Sg)+" rejected: not normal",level=verbose_rejection)
        else:
          verbose("reduction of "+str(Sf)+" by "+str(Sg)+" rejected: too large",level=verbose_rejection)
  return f

# UpdatePairs
# inputs: basis G, known new polynomials B, polynomial f, S(f)
# outputs: new S-polynomials
def UpdatePairs(G,B,f,Sf):
  Bnew = set()
  if is_primitive(f,Sf,G) and is_primitive(f,Sf,B):
    R = f.parent()
    for (g,Sg) in G:
      tf = f.lm()
      tg = g.lm()
      t = tf.lcm(tg)
      uf = R.monomial_quotient(t,tf)
      ug = R.monomial_quotient(t,tg)
      # normal pair?
      siguf = (uf*Sf[0],Sf[1])
      sigug = (ug*Sg[0],Sg[1])
      if siguf != sigug:
        if is_normal(f,Sf,uf,G,B) and is_normal(g,Sg,ug,G,B):
          if sigsort(siguf,sigug) == 1:
            sig = siguf
          else:
            sig = sigug
          verbose("generating "+str(sig)+" from "+str(Sf)+" and "+str(Sg),level=verbose_generation)
          S = uf*f - ug*g
          S *= S.lc()**(-1)
          Bnew.add((S,sig))
        else:
          verbose("pair "+str(Sf)+", "+str(Sg)+" rejected: not normal",level=verbose_rejection)
      else:
        verbose("pair "+str(Sf)+", "+str(Sg)+" rejected: equal signature",level=verbose_rejection)
  return Bnew
  
# prune_not_normal
# prunes polynomials of signature that are not normal
# inputs list of polynomials B, list of polynomials G
# outputs: pruned version of B
def prune_not_normal(B,G):
  Bnew = set()
  while len(B) != 0:
    f,Sf = B.pop()
    if is_normal(f,Sf,Integer(1),G,Bnew.union(B)):
      Bnew.add((f,Sf))
  return list(Bnew)

# f5_arri
# compute groebner basis
# inputs: list of polynomials F
# outputs: groebner basis (we hope) of ideal of F
def f5_arri(F):
  L = set()
  G = set()
  R = F[0].parent()
  B = [(F[i]*F[i].lc()**(-1),(R(1),i)) for i in xrange(len(F))]
  while len(B) != 0:
    Bnew = set()
    for (f,Sf) in B:
      if not any(Sg[1]==Sf[1] and R.monomial_divides(Sg[0],Sf[0]) for Sg in L):
        Bnew.add((f,Sf))
    B = Bnew
    #i = len(B)
    B = prune_rewritables(R,G,B)
    #if len(B) < i:
    #  print "pruned", i - len(B)
    B = prune_not_normal(B,G)
    B.sort(polysort)
    f,Sf = B.pop(0)
    while len(B) != 0 and B[0][1] == Sf:
      B.pop(0)
    f = Sreduce(f,Sf,G,B)
    verbose("reduced "+str(Sf)+" to "+str(f.lm()),level=verbose_reduction)
    if f != 0:
      if is_primitive(f,Sf,G) and is_primitive(f,Sf,B):
        B.extend(UpdatePairs(G,B,f,Sf))
      G.add((f,Sf))
    else:
      L.add(Sf)
      verbose("zero reduction: "+str(Sf),level=verbose_zero_reduction)
  print "There were", len(L), "reductions to zero."
  return G