class PuiseuxTerm:
   """
   Container for data produced at each stage of the Newton recursion.
   """
   def __init__(self, p, q, r, coeff, remainder):
      self.p, self.q, self.r = p, q, r
      self.coeff, self.remainder = coeff, remainder
      
class PuiseuxTrunk:
   """
   Iterator for Puiseux expansions.
   Initialize with a 2-variable polynomial.
   Each iteration returns a list of PuiseuxTerms.
   """
   def __init__(self, f):
      self.field = f.base_ring()
      self.ring = self.field['x','y']
      self.poly = self.ring(f)
      
   def newton_normals(self):
      """
      Return the reduced integer normal vectors to the lower
      edges of the Newton polygon of a 2-variable polynomial f.
      """
      halfspaces = self.poly.newton_polytope().Hrepresentation()
      normals = [h.vector()[1:] for h in halfspaces]
      # I believe these always point into the polyhedron.
      # Maybe they are always primitive integer vectors --
      # I was not able to find out.
      result = []
      for v in normals:
         if not v > 0:
            continue
         v *= v.denominator()
         v  = v.change_ring(ZZ)
         v /= gcd(v)
         result.append(v)
      return result

   def edge_step(self, p, q):
      f = self.poly
      result = []
      F = f.base_ring()
      s = F['s'].gen()
      t = S['t'].gen()
      e = f(t**p,s*t**q)
      m = e.valuation()
      e = e.coefficients()[0]
      # remove s factors
      e = e.shift(-e.valuation())
      for poly, multiplicity in e.factor():
         K = F.extension(poly, 'a'); a = K.gen()
         K_xy = PolynomialRing(K, ['x','y'])
         x, y = K_xy.gens()
         g = K_xy( f(x**p,x**q*(a+y))/x**m )
         r = min([d[1] for d in g.exponents()])
         result.append( (a, p, q, K_xy(g/y**r)) )
      return result