"""
R.<x,y> = PolynomialRing(QQ)
f = x^4 + y^4 + x^2 + y^2
G = [x^2 + 1, x*y^2 + y^3]
div(f,G)
GG = gb(G); GG
div(f,GG)
"""

LCM = lambda a, b: a.lcm(b)

def S(f,g):
   lf, lg = LT(f), LT(g)
   R = parent(f)
   M = LCM(lf,lg)
   return (M//lf)*f - (M//lg)*g

LT = lambda f: f.lt()

def divides(M,N):
	q, r = N.quo_rem(M)
	if r == 0:
		return True
	else:
		return False

def div(f,g):
   n = len(g)
   a = [0 for x in range(0,n)]
   r = 0
   p = f
   while p != 0:
      i = 0
      divisionoccured = False
      while i < n and divisionoccured == False:
         if divides( LT(g[i]), LT(p) ):
            a[i] = a[i] + LT(p)//LT(g[i])
            p = p - (LT(p)//LT(g[i]))*g[i]
            divisionoccured = True
         else:
            i = i + 1
      if divisionoccured == False:
         r = r + LT(p)
         p = p - LT(p)
   return a, r

def setequal(A,B):
   if len(A) != len(B):
      return False
   else:
      result = True
      for x in A:
        if x not in B:
           result = False
           break
      return result

def SREM(f,g,H):
   spoly = S(f,g)
   return div(spoly, H)

def gb(F):
   G = [p for p in F]
   nn, n = -1, -2
   while nn != n:
      n = len(G)
      for i in range(0,n):
         for j in range(i+1,n):
            spoly = S(G[i], G[j])
            spoly = div(spoly,G)[1]
            if spoly != 0:
               G += [spoly]
      nn = len(G)
   return G

RL = lambda a, b: map( lambda x: R.random_element(x), [ a for i in range(0,b) ])