text: 
  from sympy.core.add import Add
  from sympy.core.mul import Mul
  from sympy.core.power import Pow
  from sympy.core.sympify import sympify
  from sympy.core.basic import Basic, S, C
  from sympy.core.symbol import Symbol, Wild
  from sympy.core.numbers import Integer, igcd, ilcm
  
  from sympy.utilities import all, any
  
  from sympy.polys.monomial import monomial_cmp, monomial_mul, \
      monomial_div, monomial_max, monomial_min, monomial_as_basic
  
  import sympy.polys
  
  import math
  
  class SymbolsError(Exception):
      pass
  
  class PolynomialError(Exception):
      pass
  
  class CoefficientError(PolynomialError):
      pass
  
  class UnivariatePolyError(PolynomialError):
  
      def __init__(self, poly):
          self.message = "%s is an univariate polynomial" % poly
  
      def __str__(self):
          return self.message
  
  class MultivariatePolyError(PolynomialError):
  
      def __init__(self, poly):
          self.message = "%s is a multivariate polynomial" % poly
  
      def __str__(self):
          return self.message
  
  ##
  ## TODO:
  ##
  ##  [!] Copy + Write tests for everything.
  ##  [2] Analyze dense MV multiplication.
  ##  [~] Create MV poly with int coeffs.
  ##  [4] Think on GF polys (for factoring).
  ##  [5] Improve coefficients analysis.
  ##  [6] Implement FactoredPoly.
  ##  [7] Concept of Monomial.
  ##
  
  class Poly(Basic):
      """Represents polynomials with symbolic coefficients.
  
         Polynomials are internally represented as two lists containing
         coefficients and monomials (tuples of exponents) respectively.
         Stored are only terms with non-zero coefficients, so generally
         all polynomials are considered sparse. However algorithms will
         detect dense polynomials and use the appropriate method to
         solve the given problem in the most efficient way.
  
         The most common way to initialize a polynomial instance is to
         provide a valid expression together witch a set of ordered
         symbols and, additionally, monomial ordering:
  
              Poly(expression, x_1, x_2, ..., x_n, order='grlex')
  
         If the given expression is not a polynomial with respect to
         the given variables, an exception is raised. Alternatively
         you can use Basic.as_poly to avoid exception handling.
  
         By default ordering of monomials can be omitted. In this case
         graded lexicographic order will be used. Anyway remember that
         'order' is a keyword argument.
  
         Currently there are supported four standard orderings:
  
             [1] lex       -> lexicographic order
             [2] grlex     -> graded lex order
             [3] grevlex   -> reversed grlex order
             [4] 1-el      -> first elimination order
  
         Polynomial can be also constructed explicitly by specifying
         a collection of coefficients and monomials. This can be done
         in at least three different ways:
  
             [1] [(c_1, M_1), (c_2, M_2), ..., (c_1, M_1)]
  
             [2] (c_1, c2, ..., c_n), (M_1, M_2, ..., M_n)
  
             [3] { M_1 : c_1, M_2 : c_2, ..., M_n : c_n }
  
         Although all three representation look similar, they are
         designed for different tasks and have specific properties:
  
             [1] All coefficients and monomials  are validated before
                 polynomial instance is created. Monomials are sorted
                 with respect to the given order.
  
             [2] No validity checking, no sorting, just a raw copy.
  
             [3] Also no validity checking however monomials are
                 sorted with respect to the given order.
  
         For interactive usage choose [1] as it's safe and relatively
         fast. Use [2] or [3] internally for time critical algorithms,
         when you know that coefficients and monomials will be valid.
  
         Implemented methods:
  
             U --> handles uniformly int | long and Basic coefficients
                   (other methods need to be overridden in IntegerPoly)
  
             P --> a property (with appropriate decorator)
  
             - --> otherwise
  
         [1] Representation conversion:
  
            [1.1] [--] as_basic   --> converts to a valid sympy expression
            [1.2] [U-] as_dict    --> returns { M_1 : c_1, ..., M_1 : c_1 }
            [1.3] [U-] as_uv_dict --> returns dict with integer keys
  
         [2] Polynomial internals:
  
            [2.1] [UP] coeffs  --> list of coefficients
            [2.2] [UP] monoms  --> list of monomials
            [2.3] [UP] symbols --> an ordered list of symbols
            [2.4] [UP] order   --> the ordering of monomials
            [2.5] [UP] stamp   --> frozen set of polynomial's variables
            [2.6] [UP] domain  --> coefficient domain, see [8] for details
            [2.7] [UP] args    --> required for printing and hashing
            [2.8] [UP] flags   --> dict of keyword arguments to Poly
  
         [3] General characteristics:
  
            [3.1] [UP] norm         --> computes oo-norm of coefficients
            [3.2] [UP] length       --> counts the total number of terms
            [3.3] [UP] degree       --> returns degree of the leading monomial
            [3.4] [UP] total_degree --> returns true total degree of a polynomial
  
         [4] Items generators:
  
            [4.1] [U-] iter_coeffs     --> iterate over coefficients
            [4.2] [U-] iter_monoms     --> iterate over monomials
            [4.3] [U-] iter_terms      --> iterate over terms
  
            [4.4] [U-] iter_all_coeffs --> iterate over all coefficients
            [4.5] [U-] iter_all_monoms --> iterate over all monomials
            [4.6] [U-] iter_all_terms  --> iterate over all terms
  
         [5] Leading and trailing items:
  
            [5.1] [UP] lead_coeff or LC --> returns the leading coefficient
            [5.2] [UP] lead_monom or LM --> returns the leading monomial
            [5.3] [UP] lead_term  or LT --> returns the leading term
  
            [5.4] [UP] tail_coeff or TC --> returns the trailing coefficient
            [5.5] [UP] tail_monom or TM --> returns the trailing monomial
            [5.6] [UP] tail_term  or TT --> returns the trailing term
  
         [6] Arithmetic operators:
  
            [6.1]  [U-] __abs__     --> for all i, abs(c_i)
            [6.2]  [U-] __neg__     --> polynomial negation
            [6.3]  [U-] __add__     --> polynomial addition
            [6.4]  [U-] __sub__     --> polynomial subtraction
            [6.5]  [--] __mul__     --> polynomial multiplication
            [6.6]  [U-] __pow__     --> polynomial exponentiation
            [6.7]  [U-] __div__     --> polynomial quotient
            [6.8]  [U-] __mod__     --> polynomial remainder
            [6.9]  [U-] __divmod__  --> both 'div' and 'mod'
            [6.10] [U-] __truediv__ --> the same as 'div'
  
         [7] Polynomial properties:
  
            [7.1]  [UP] is_zero          --> only one term c_0 == 0
            [7.2]  [UP] is_one           --> only one term c_0 == 1
            [7.3]  [-P] is_number        --> only numeric constant term
            [7.4]  [UP] is_constant      --> only arbitrary constant term
            [7.5]  [UP] is_monomial      --> number of terms == 1
            [7.6]  [UP] is_univariate    --> number of variables == 1
            [7.7]  [UP] is_multivariate  --> number of variables > 1
            [7.8]  [UP] is_homogeneous   --> has constant term
            [7.9]  [UP] is_inhomogeneous --> no constant term
            [7.10] [UP] is_sparse        --> filled with less than 90% of terms
            [7.11] [UP] is_dense         --> filled with more than 90% of terms
            [7.12] [UP] is_monic         --> returns True if leading coeff is one
            [7.13] [UP] is_primitive     --> returns True if content is one
            [7.14] [UP] is_squarefree    --> no factors of multiplicity >= 2
            [7.15] [UP] is_linear        --> all monomials of degree <= 1
  
         [8] Coefficients properties:
  
            [8.1] [UP] has_integer_coeffs  --> for all i, c_i in Z
            [8.2] [UP] has_rational_coeffs --> for all i, c_i in Q
            [8.3] [UP] has_radical_coeffs  --> for all i, c_i in R
            [8.4] [UP] has_complex_coeffs  --> for all i, c_i in C
            [8.5] [UP] has_symbolic_coeffs --> otherwise
  
         [9] Special functionality:
  
            [9.1]  [--] as_monic      --> divides all coefficients by LC
            [9.2]  [--] as_integer    --> returns polynomial with integer coeffs
  
            [9.3]  [-P] content       --> returns GCD of polynomial coefficients
            [9.4]  [--] as_primitive  --> returns content and primitive part
  
            [9.5]  [U-] as_squarefree --> returns square-free part of a polynomial
  
            [9.6]  [U-] as_reduced    --> extract GCD of polynomial monomials
  
            [9.7]  [--] map_coeffs    --> applies a function to all coefficients
            [9.8]  [U-] coeff         --> returns coefficient of the given monomial
  
            [9.9]  [U-] unify_with    --> returns polys with a common set of symbols
  
            [9.10] [U-] diff          --> efficient polynomial differentiation
            [9.11] [U-] integrate     --> efficient polynomial integration
  
            [9.12] [U-] invert        --> computes multiplicative inverse in K[t]
  
         [10] Operations on terms:
  
            [10.1] [U-] add_term --> add a term to a polynomial
            [10.2] [U-] sub_term --> subtract a term from a polynomial
  
            [10.3] [--] mul_term --> multiply a polynomial with a term
            [10.4] [--] div_term --> divide a polynomial by a term
  
            [10.5] [U-] remove_lead_term --> remove leading term (up to order)
            [10.6] [U-] remove_last_term --> remove last term (up to order)
  
         [11] Substitution and evaluation:
  
            [11.1] [U-] __call__    --> evaluates poly a the given point
            [11.2] [U-] evaluate    --> evaluates poly for specific vars
            [11.3] [--] _eval_subs  --> efficiently substitute variables
  
         [12] Comparison functions:
  
            [12.1] [U-] __eq__      --> P1 == P, up to order of monomials
            [12.2] [U-] __ne__      --> P1 != P, up to order of monomials
            [12.3] [U-] __nonzero__ --> check if zero polynomial
  
         [13] String representation functions:
  
            [13.1] [U-] repr      --> returns represantation of 'self'
            [13.3] [U-] str       --> returns pretty representation
  
         [14] Other (derived from Basic) functionality:
  
            [14.1] [U-] _eval_is_polynomial --> checks if poly is a poly
  
         For general information on polynomials and algorithms, refer to:
  
         [1] D. E. Knuth, The Art of Computer Programming: Seminumerical
             Algorithms, v.2, Addison-Wesley Professional, 1997
  
         [2] J. von zur Gathen, J. Gerhard, Modern Computer Algebra,
             Second Edition, Cambridge University Press, 2003
  
         [3] K. Geddes, S. R. Czapor, G. Labahn, Algorithms for
             Computer Algebra, Springer, 1992
  
         [4] D. Bini, V. Y.Pan, Polynomial and Matrix Computations:
             Volume 1: Fundamental Algorithms, Birkhauser Boston, 1994
  
         [5] R. T. Moenck, Practical Fast Polynomial Multiplication,
             Proceedings of the third ACM symposium on Symbolic and
             algebraic computation, ACM Press, Yorktown Heights,
             New York, USA, 1976, pp. 136-148
  
         [6] M. Bronstein, Symbolic Integration I: Transcendental
             Functions, Second Edition, Springer-Verlang, 2005
  
         See also documentation of particular methods.
  
      """
  
      def __new__(cls, poly, *symbols, **flags):
          N, terms = len(symbols), {}
          coeffs, monoms = (), ()
  
          order = flags.pop('order', 'grlex')
  
          if N == 0:
              if isinstance(poly, Poly):
                  if order != poly.order:
                      symbols = poly.symbols
                      N = len(poly.symbols)
                      poly = poly.as_dict()
                  else:
                      return poly
              else:
                  raise SymbolsError("No symbols were given")
  
          stamp = frozenset(symbols)
  
          if len(stamp) != N:
              raise SymbolsError("Duplicate symbols: %s" % (symbols,))
  
          symbols = tuple(map(sympify, symbols))
  
          if any(not s.is_Symbol for s in symbols):
              raise SymbolsError("Invalid symbols: %s" % (symbols,))
  
          if any(not s.is_commutative for s in symbols):
              raise SymbolsError("Non-commutative symbols: %s" % (symbols,))
  
          # { M1: c1, M2: c2, ... }
          if type(poly) is dict:
              if not poly:
                  coeffs = (S.Zero,)
                  monoms = ((0,) * N,)
              else:
                  terms = poly
  
          # ((c1, c2, ...), (M1, M2, ...))
          elif type(poly) is tuple:
              if not poly:
                  coeffs = (S.Zero,)
                  monoms = ((0,) * N,)
              else:
                  if len(poly) == 2:
                      coeffs, monoms = poly[0], poly[1]
  
                      if isinstance(coeffs, (tuple, list)):
                          if not (coeffs or monoms):
                              coeffs = (S.Zero,)
                              monoms = ((0,) * N,)
                      else:
                          coeffs = (coeffs,)
                          monoms = (monoms,)
                  else:
                      raise PolynomialError("Invalid arguments," \
                          " should get '(coeffs, monoms)' tuple")
  
          # [ (c1,M1), (c2,M2), ... ]
          elif type(poly) is list:
              if not poly:
                  coeffs = (S.Zero,)
                  monoms = ((0,) * N,)
              else:
                  for coeff, monom in poly:
                      coeff = sympify(coeff)
  
                      if coeff.has_any_symbols(*symbols):
                          raise CoefficientError("%s coefficient is dependent" \
                              " of polynomial's symbols %s" % (coeff, symbols))
  
                      if type(monom) is int:
                          monom = (monom,)
  
                      if len(monom) != N:
                          raise PolynomialError("Dimension of %s monomial does" \
                              " not match the number of symbols (%d)" % (monom, N))
  
                      if not coeff:
                          continue
  
                      if terms.has_key(monom):
                          coeff += terms[monom]
  
                          if not coeff:
                              del terms[monom]
                              continue
  
                      terms[monom] = coeff
  
                  if not terms:
                      coeffs = (S.Zero,)
                      monoms = ((0,) * N,)
          elif isinstance(poly, Poly):
              if stamp <= poly.stamp:
                  if symbols == poly.symbols:
                      if N > 1 and order != poly.order:
                          terms = poly.as_dict()
                      else:
                          return poly
                  else:
                      terms = Poly._permute(poly, *symbols)
              else:
                  if any(coeff.has_any_symbols(*symbols) for coeff in poly.coeffs):
                      terms = Poly._decompose(poly.as_basic(), *symbols)
                  else:
                      K = len(poly.symbols)
  
                      if symbols[:K] == poly.symbols:
                          coeffs, T = poly.coeffs, (0,) * (N - K)
                          monoms = tuple(M + T for M in poly.monoms)
  
                          if order != poly.order:
                              terms = dict(zip(monoms, coeffs))
                      else:
                          terms = Poly._permute(poly, *symbols)
          else:
              terms = Poly._decompose(poly, *symbols)
  
          if terms:
              if N == 1 and type(terms.keys()[0]) is not tuple:
                  keys = tuple(reversed(sorted(terms.keys())))
  
                  coeffs = [ terms[i] for i in keys ]
                  monoms = [ (i,) for i in keys ]
              else:
                  f = monomial_cmp(order)
  
                  monoms = terms.keys()
                  monoms.sort(f, reverse=True)
  
                  coeffs = [ terms[monom] for monom in monoms ]
  
          args = (tuple(coeffs), tuple(monoms),
                  symbols, order, stamp, None)
  
          return Basic.__new__(cls, *args, **flags)
  
      def __getnewargs__(self):
          coeffs, monoms, symbols, order = self.args[:4]
          return ((coeffs,monoms),) + symbols
  
      @staticmethod
      def _permute(poly, *symbols):
          terms, N = {}, len(symbols)
  
          if symbols == poly.symbols[:N]:
              symbols = poly.symbols[N:]
  
              for coeff, monom in poly.iter_terms():
                  coeff *= monomial_as_basic(monom[N:], *symbols)
  
                  monom = monom[:N]
  
                  if terms.has_key(monom):
                      coeff += terms[monom]
  
                      if not coeff:
                          del terms[monom]
                          continue
  
                  terms[monom] = coeff
          else:
              new_symbols, indices = list(poly.symbols), []
  
              for s in new_symbols[:]:
                  for i, t in enumerate(symbols):
                      if s is t:
                          new_symbols.remove(s)
                          indices.append(i)
                          break
                  else:
                      indices.append(None)
  
              for coeff, monom in poly.iter_terms():
                  exponents, j = [0] * N, 0
  
                  for exp, i in zip(monom, indices):
                      if i is None:
                          coeff *= new_symbols[j]**exp
                          j += 1
                      else:
                          exponents[i] = exp
  
                  monom = tuple(exponents)
  
                  if terms.has_key(monom):
                      coeff += terms[monom]
  
                      if not coeff:
                          del terms[monom]
                          continue
  
                  terms[monom] = coeff
  
          return terms
  
      @staticmethod
      def _decompose(poly, *symbols):
          """Converts basic expression to dictionary representation.
  
             This procedure will expand given expression and scan it,
             extracting coefficients and monomials according to a set
             of symbols, if possible. At this point there is no need
             for specifying ordering information. If conversion is
             not possible an exception is raised.
  
             This method should be left for internal use. To convert
             expression to a polynomial use Polynomial constructor
             explicitly or Basic.as_poly method.
  
             >>> from sympy import *
             >>> x, y = symbols('xy')
  
             >>> Poly._decompose(y*x**2 + 3, x)
             {(0,): 3, (2,): y}
  
             >>> Poly._decompose(y*x**2 + 3, x, y)
             {(0, 0): 3, (2, 1): 1}
  
          """
          result, indices, N = {}, {}, len(symbols)
  
          for i, sym in enumerate(symbols):
              indices[sym] = i
  
          poly = sympify(poly).expand()
  
          if poly.is_Add:
              terms = poly.args
          else:
              if poly.is_Number:
                  return { (0,) * N : poly }
              else:
                  terms = [ poly ]
  
          for term in terms:
              if not term.has_any_symbols(*symbols):
                  coeff, monom = term, (0,) * N
              else:
                  if term.is_Mul:
                      factors = term.args
                  else:
                      factors = [ term ]
  
                  coeff, monom = S.One, [0] * N
  
                  for factor in factors:
                      if factor.has_any_symbols(*symbols):
                          if factor.is_Pow:
                              if factor.exp.is_Integer:
                                  b, e = factor.base, factor.exp.p
  
                                  if b.is_Symbol and e > 0:
                                      monom[indices[b]] += e
                                      continue
                          elif factor.is_Symbol:
                              monom[indices[factor]] += 1
                              continue
  
                          raise PolynomialError("Can't decompose %s" % factor)
                      else:
                          coeff *= factor
  
                  monom = tuple(monom)
  
              if result.has_key(monom):
                  coeff += result[monom]
  
                  if not coeff:
                      del result[monom]
                      continue
  
              result[monom] = coeff
  
          if not result:
              return { (0,) * N : S.One }
          else:
              return result
  
      @staticmethod
      def cancel(f, *symbols):
          """Cancel common factors in a fractional expression.
  
             Given a quotient of polynomials, cancel common factors from
             the numerator and the denominator.  The result is formed in
             an expanded form, even if there was no cancellation.
  
             To improve the  speed of calculations a list of symbols can
             be specified. This can be particularly useful in univariate
             case with additional parameters, to avoid Groebner basis
             computations.
  
             The input fraction can be given as a single expression or
             as a tuple consisting of the numerator and denominator.
  
             >>> from sympy import *
             >>> x,y = symbols('xy')
  
             >>> Poly.cancel((x**2-y**2)/(x-y), x, y)
             x + y
  
             >>> Poly.cancel((x**2-y**2, x-y), x, y)
             x + y
  
          """
          if type(f) is tuple:
              numer, denom = f
          else:
              if f.is_Atom:
                  return f
  
              if not symbols:
                  symbols = f.atoms(Symbol)
  
                  if not symbols:
                      symbols = None
                  else:
                      symbols = sorted(symbols)
  
              numer, denom = f.as_numer_denom()
  
          numer = numer.expand()
  
          if numer is S.Zero:
              return numer
  
          denom = denom.expand()
  
          if denom.is_Number:
              return numer / denom
  
          if symbols is None:
              if denom.is_Add:
                  a, b, c = map(Wild, 'abc')
  
                  r = denom.match(a + b*c**S.Half)
  
                  if r is not None:
                      if r[b] is S.Zero:
                          denom = match[a]
                      else:
                          a, b, c = r[a], r[b], r[c]
  
                          numer *= a-b*c**S.Half
                          numer = numer.expand()
  
                          denom = a**2 - c*b**2
  
              return numer / denom
  
          try:
              p = Poly(numer, *symbols)
              q = Poly(denom, *symbols)
          except PolynomialError:
              return numer / denom
  
          g = sympy.polys.algorithms.poly_gcd(p, q)
  
          p = sympy.polys.algorithms.poly_div(p, g)[0]
          q = sympy.polys.algorithms.poly_div(q, g)[0]
  
          pc = int(p.content)
          qc = int(q.content)
  
          cont = igcd(pc, qc)
  
          if cont != 1:
              p = p.div_term(cont)
              q = q.div_term(cont)
  
          return p.as_basic() / q.as_basic()
  
      def as_basic(self):
          """Converts polynomial to a valid sympy expression.
  
             Takes list representation of a polynomial and constructs
             expression using symbols used during creation of a given
             polynomial. If you wish to have different symbols in the
             resulting expressions, use subs() first.
  
             Note that the result is always returned in expanded form.
  
             >>> from sympy import *
             >>> x,y = symbols('xy')
  
             >>> p = Poly(x**2+3, x)
             >>> p.as_basic()
             3 + x**2
  
          """
          multinomial = dict(zip(self.monoms, self.coeffs))
          return multinomial_as_basic(multinomial, *self.symbols)
  
      def as_dict(self):
          """Convert list representation to dictionary representation.
  
             Due to hashing of immutable expressions in Python, we can
             not store a dictionary directly in a polynomial instance.
             Also it would be inconvenient in case of ordering of
             monomials (unnecessary sorting).
  
             This method creates dictionary efficiently, so it can be
             used in several algorithms which perform best when dicts
             are being used, eg. addition, multiplication etc.
  
             >>> from sympy import *
             >>> x,y = symbols('xy')
  
             >>> p = Poly(x**2+3, x)
             >>> p.as_dict()
             {(0,): 3, (2,): 1}
  
          """
          return dict(zip(self.monoms, self.coeffs))
  
      def as_uv_dict(self):
          """Return dictionary representation with integer keys.
  
             In case of univariate polynomials it is inefficient to
             to handle exponents as singleton tuples, as an example
             see multiplication algorithm.  This method will return
             dictionary representation with all tuples converted to
             explicit integers.
  
             >>> from sympy import *
             >>> x,y = symbols('xy')
  
             >>> p = Poly(x**2+3, x)
             >>> p.as_uv_dict()
             {0: 3, 2: 1}
  
          """
          if self.is_univariate:
              return dict(zip([ M[0] for M in self.monoms ], self.coeffs))
          else:
              raise MultivariatePolyError(self)
  
      @property
      def coeffs(self):
          return self._args[0]
  
      @property
      def monoms(self):
          return self._args[1]
  
      @property
      def symbols(self):
          return self._args[2]
  
      @property
      def order(self):
          return self._args[3]
  
      @property
      def stamp(self):
          return self._args[4]
  
      @property
      def domain(self):
          return self._args[5]
  
      @property
      def args(self):
          return self._args[0:4]
  
      @property
      def flags(self):
          return { 'order' : self.order }
  
      @property
      def length(self):
          return len(self.monoms)
  
      @property
      def degree(self):
          """Returns degree of the leading term. """
          return sum(self.monoms[0])
  
      @property
      def total_degree(self):
          """Returns true total degree of a polynomial. """
          return max([ sum(monom) for monom in self.monoms ])
  
      @property
      def norm(self):
          """Computes oo-norm of polynomial's coefficients. """
          return max([ abs(coeff) for coeff in self.coeffs ])
  
      def iter_basic_args(self):
          for coeff in self.coeffs:
              yield coeff
  
          for i, symbol in enumerate(self.symbols):
              if any(monom[i] != 0 for monom in self.monoms):
                  yield symbol
  
      def iter_coeffs(self):
          for coeff in self.coeffs:
              yield coeff
  
      def iter_monoms(self):
          for monom in self.monoms:
              yield monom
  
      def iter_terms(self):
          for coeff, monom in zip(self.coeffs, self.monoms):
              yield coeff, monom
  
      def iter_all_coeffs(self):
          if self.is_univariate:
              terms = self.as_uv_dict()
  
              for i in xrange(self.degree, -1, -1):
                  if terms.has_key(i):
                      yield terms[i]
                  else:
                      yield S.Zero
          else:
              raise MultivariatePolyError(self)
  
      def iter_all_monoms(self):
          if self.is_univariate:
              for i in xrange(self.degree, -1, -1):
                  yield (i,)
          else:
              raise MultivariatePolyError(self)
  
      def iter_all_terms(self):
          if self.is_univariate:
              terms = self.as_uv_dict()
  
              for i in xrange(self.degree, -1, -1):
                  if terms.has_key(i):
                      yield terms[i], (i,)
                  else:
                      yield S.Zero, (i,)
          else:
              raise MultivariatePolyError(self)
  
      @property
      def lead_coeff(self):
          return self.coeffs[0]
  
      @property
      def lead_monom(self):
          return self.monoms[0]
  
      @property
      def lead_term(self):
          return (self.coeffs[0], self.monoms[0])
  
      LC = lead_coeff
      LM = lead_monom
      LT = lead_term
  
      @property
      def tail_coeff(self):
          return self.coeffs[-1]
  
      @property
      def tail_monom(self):
          return self.monoms[-1]
  
      @property
      def tail_term(self):
          return (self.coeffs[-1], self.monoms[-1])
  
      TC = tail_coeff
      TM = tail_monom
      TT = tail_term
  
      def __abs__(self):
          return self.__class__(([ abs(coeff) for coeff in self.coeffs ],
              self.monoms), *self.symbols, **self.flags)
  
      def __neg__(self):
          return self.__class__(([ -coeff for coeff in self.coeffs ],
              self.monoms), *self.symbols, **self.flags)
  
      def __add__(self, other):
          try:
              self, poly = self.unify_with(other)
          except PolynomialError:
              return self.as_basic() + other.as_basic()
  
          if poly.is_zero:
              return self
  
          if self.is_zero:
              return poly
  
          if poly.is_monomial:
              return self.add_term(*poly.LT)
  
          if self.is_monomial:
              return poly.add_term(*self.LT)
  
          terms = self.as_dict()
  
          for coeff, monom in poly.iter_terms():
              if terms.has_key(monom):
                  coeff += terms[monom]
  
                  coeff = Poly.cancel(coeff)
  
                  if not coeff:
                      del terms[monom]
                      continue
  
              terms[monom] = coeff
  
          return self.__class__(terms, *self.symbols, **self.flags)
  
      def __sub__(self, other):
          try:
              self, poly = self.unify_with(other)
          except PolynomialError:
              return self.as_basic() - other.as_basic()
  
          if poly.is_zero:
              return self
  
          if self.is_zero:
              return -poly
  
          if poly.is_monomial:
              return self.sub_term(*poly.LT)
  
          if self.is_monomial:
              return (-poly).add_term(*self.LT)
  
          terms = self.as_dict()
  
          for coeff, monom in poly.iter_terms():
              if terms.has_key(monom):
                  coeff -= terms[monom]
  
                  coeff = Poly.cancel(coeff)
  
                  if not coeff:
                      del terms[monom]
                      continue
  
              terms[monom] = -coeff
  
          return self.__class__(terms, *self.symbols, **self.flags)
  
      def __mul__(self, other):
          """Efficiently multiply sparse and dense polynomials.
  
             Given polynomials P and Q, if both of them are dense, then
             in univariate case perform dense multiplication, otherwise
             apply Kronecker 'trick', mapping multivariate polynomials
             to univariate ones in last variable and then multiply.
  
             If any of the polynomials is sparse then in both univariate
             and multivariate cases use naive multiplication algorithm.
  
             For more information on implemented algorithms refer to:
  
             [1] R. T. Moenck, Practical Fast Polynomial Multiplication,
                 Proceedings of the third ACM symposium on Symbolic and
                 algebraic computation, ACM Press, Yorktown Heights,
                 New York, USA, 1976, pp. 136-148
          """
          try:
              self, poly = self.unify_with(other)
          except PolynomialError:
              return self.as_basic() * other.as_basic()
  
          if self.is_constant:
              if self.is_zero:
                  return self
              elif self.is_one:
                  return poly
              else:
                  return poly.mul_term(self.LC)
  
          if poly.is_constant:
              if poly.is_zero:
                  return poly
              elif poly.is_one:
                  return self
              else:
                  return self.mul_term(poly.LC)
  
          if self.is_monomial:
              return poly.mul_term(*self.LT)
  
          if poly.is_monomial:
              return self.mul_term(*poly.LT)
  
          if self.is_dense and poly.is_dense:
              if self.is_multivariate:
                  a = monomial_max(*self.monoms)
                  b = monomial_max(*poly.monoms)
  
                  degs, p, q = [1], {}, {}
  
                  for x, y in reversed(zip(a, b)[1:]):
                      degs.insert(0, degs[0]*(x+y+1))
  
                  N = sum([ x*y for x, y in zip(self.LM, degs) ])
                  M = sum([ x*y for x, y in zip(poly.LM, degs) ])
  
                  for coeff, monom in self.iter_terms():
                      p[sum([ x*y for x, y in zip(monom, degs) ])] = coeff
  
                  for coeff, monom in poly.iter_terms():
                      q[sum([ x*y for x, y in zip(monom, degs) ])] = coeff
              else:
                  p, N = self.as_uv_dict(), self.degree
                  q, M = poly.as_uv_dict(), poly.degree
  
              coeffs, monoms = [], []
  
              for k in xrange(N+M, -1, -1):
                  coeff = 0
  
                  for i in xrange(k+1):
                      if p.has_key(i) and q.has_key(k-i):
                          coeff += Poly.cancel(p[i] * q[k-i])
  
                  if coeff:
                      coeffs.append(coeff)
                      monoms.append((k,))
  
              if self.is_multivariate:
                  terms, degs = {}, degs[:-1]
  
                  for i, d in enumerate(monoms):
                      monom, d = [], d[0]
  
                      for j in degs:
                          k, d = divmod(d, j)
                          monom.append(k)
  
                      terms[tuple(monom)+(d,)] = coeffs[i]
              else:
                  terms = coeffs, monoms
          else:
              terms = {}
  
              for coeff_p, M in self.iter_terms():
                  for coeff_q, N in poly.iter_terms():
                      monom = monomial_mul(M, N)
  
                      coeff = coeff_p * coeff_q
                      coeff = Poly.cancel(coeff)
  
                      if terms.has_key(monom):
                          coeff += terms[monom]
  
                          if not coeff:
                              del terms[monom]
                              continue
  
                      terms[monom] = coeff
  
          return self.__class__(terms, *self.symbols, **self.flags)
  
      def __pow__(self, other):
          """Polynomial exponentiation using binary method.
  
             Given a polynomial P and integer exponent N, raise P to power
             N in floor(lg(N)) + 1(N) multiplications,  where lg(N) stands
             for binary logarithm and 1(N) is the number of ones in binary
             representation of N.
  
             For more information on the implemented algorithm refer to:
  
             [1] D. E. Knuth, The Art of Computer Programming: Seminumerical
                 Algorithms, v.2, Addison-Wesley Professional, 1997, pp. 461
  
          """
          other = sympify(other)
  
          if isinstance(other, Poly):
              if other.is_constant:
                  other = sympify(other.lead_coeff)
              else:
                  return self.as_basic()**other.as_basic()
  
          if other.is_Integer:
              if other is S.One:
                  return self
  
              P = self.__class__(S.One, *self.symbols, **self.flags)
  
              if other is S.Zero:
                  return P
  
              N, Q = abs(int(other)), self
  
              while True:
                  N, M = N//2, N
  
                  if M & 1:
                      P *= Q
  
                      if N == 0:
                          break
  
                  Q *= Q
  
              if other.is_positive:
                  return P
              else:
                  return Pow(P, S.NegativeOne)
          else:
              return self.as_basic()**other
  
      def __div__(self, other):
          try:
              return sympy.polys.algorithms.poly_div(self, other)[0]
          except PolynomialError:
              return self.as_basic() / other.as_basic()
  
      def __mod__(self, other):
          try:
              return sympy.polys.algorithms.poly_div(self, other)[1]
          except PolynomialError:
              return S.Zero
  
      def __divmod__(self, other):
          try:
              return sympy.polys.algorithms.poly_div(self, other)
          except PolynomialError:
              return (self.as_basic() / other.as_basic(), S.Zero)
  
      __truediv__ = __div__
  
      @property
      def is_zero(self):
          return self.coeffs in ((S.Zero,), (0,))
  
      @property
      def is_one(self):
          return self.coeffs in ((S.One,), (1,)) and \
              all(e == 0 for e in self.monoms[0])
  
      @property
      def is_number(self):
          return self.is_constant and self.coeffs[0].is_number
  
      @property
      def is_constant(self):
          return len(self.monoms) == 1 and \
              all(e == 0 for e in self.monoms[0])
  
      @property
      def is_monomial(self):
          return len(self.monoms) == 1
  
      @property
      def is_univariate(self):
          return len(self.symbols) == 1
  
      @property
      def is_multivariate(self):
          return len(self.symbols) != 1
  
      @property
      def is_homogeneous(self):
          return any(e != 0 for e in self.monoms[-1])
  
      @property
      def is_inhomogeneous(self):
          return all(e == 0 for e in self.monoms[-1])
  
      @property
      def is_sparse(self):
          return not self.is_dense
  
      @property
      def is_dense(self):
          """Returns True if 'self' is dense.
  
             Let 'self' be a polynomial in M variables of a total degree N
             and having L terms (with non-zero coefficients). Let K be the
             number of monomials in M variables of degree at most N.  Then
             'self' is considered dense if it's at least 90% filled.
  
             The total number of monomials is given by (M + N)! / (M! N!),
             where the factorials are estimated using improved Stirling's
             approximation:
  
                         n! = sqrt((2 n + 1/3) pi) n**n exp(-n)
  
             For more information on this subject refer to:
  
               http://mathworld.wolfram.com/StirlingsApproximation.html
  
          """
          def enum(M, N):
              U = float(M+N)
              V = (6*M+1)*(6*N+1)
              S = 2.4*math.sqrt(U/V)
              A = math.pow(U/M, M)
              B = math.pow(U/N, N)
              return S * A * B
  
          V, N = len(self.symbols), self.total_degree
  
          return (self.length / enum(V, N)) > 0.9
  
      @property
      def is_primitive(self):
          return self.content in (S.One, 1)
  
      @property
      def is_monic(self):
          return self.lead_coeff in (S.One, 1)
  
      @property
      def is_squarefree(self):
          """Returns True if 'self' has no factors of multiplicity >= 2.
  
             >>> from sympy import *
             >>> x,y = symbols('xy')
  
             >>> Poly(x-1, x).is_squarefree
             True
  
             >>> Poly((x-1)**2, x).is_squarefree
             False
  
          """
          return self.as_squarefree() is self
  
      @property
      def is_linear(self):
          return all(sum(monom) <= 1 for monom in self.iter_monoms())
  
      @property
      def has_integer_coeffs(self):
          return self.domain == 'Z'
  
      @property
      def has_rational_coeffs(self):
          return self.domain == 'Q'
  
      @property
      def has_radical_coeffs(self):
          return self.domain == 'R'
  
      @property
      def has_complex_coeffs(self):
          return self.domain == 'C'
  
      @property
      def has_symbolic_coeffs(self):
          return self.domain == 'S'
  
      def as_monic(self):
          """Returns a monic polynomial.
  
             >>> from sympy import *
             >>> x,y = symbols('xy')
  
             >>> Poly(x**2 + 4, x).as_monic()
             Poly(x**2 + 4, x)
  
             >>> Poly(2*x**2 + 4, x).as_monic()
             Poly(x**2 + 2, x)
  
             >>> Poly(y*x**2 + y**2 + y, x).as_monic()
             Poly(x**2 + 1 + y, x)
  
          """
          LC = self.lead_coeff
  
          if not self or LC is S.One:
              return self
  
          coeffs = [ coeff / LC for coeff in self.coeffs ]
  
          for i, coeff in enumerate(coeffs):
              coeffs[i] = Poly.cancel(coeff)
  
          return self.__class__((coeffs, self.monoms),
              *self.symbols, **self.flags)
  
      def as_integer(self):
          """Makes all coefficients integers if possible.
  
             Given a polynomial with rational coefficients,  returns a tuple
             consisting of a common denominator of those coefficients and an
             integer polynomial. Otherwise an exception is being raised.
  
             >>> from sympy import *
             >>> x,y = symbols("xy")
  
             >>> Poly(3*x**2 + x, x).as_integer()
             (1, Poly(3*x**2 + x, x))
  
             >>> Poly(3*x**2 + x/2, x).as_integer()
             (2, Poly(6*x**2 + x, x))
  
          """
          denom, coeffs = 1, []
  
          for coeff in self.iter_coeffs():
              if coeff.is_Rational:
                  coeffs.append(coeff)
  
                  if not coeff.is_Integer:
                      denom = ilcm(denom, coeff.q)
              elif coeff.is_Real and int(coeff) == coeff:
                  coeffs.append(Integer(int(coeff)))
              else:
                  raise CoefficientError("%s is not a rational number" % coeff)
  
          denom = sympify(denom)
  
          if denom is not S(1):
              coeffs = [ coeff * denom for coeff in self.iter_coeffs() ]
  
          return denom, self.__class__((coeffs, self.monoms),
              *self.symbols, **self.flags)
  
      @property
      def content(self):
          """Returns integer GCD of all the coefficients.
  
             >>> from sympy import *
             >>> x,y,z = symbols('xyz')
  
             >>> p = Poly(2*x + 5*x*y, x, y)
             >>> p.content
             1
  
             >>> p = Poly(2*x + 4*x*y, x, y)
             >>> p.content
             2
  
             >>> p = Poly(2*x + z*x*y, x, y)
             >>> p.content
             1
  
          """
          if self.is_zero:
              return S.One
          else:
              content = 0
  
              for coeff in self.coeffs:
                  if coeff.is_Rational:
                      content = igcd(content, coeff.p)
                  else:
                      return S.One
  
              return Integer(content)
  
      def as_primitive(self):
          """Returns content and primitive part of a polynomial.
  
             >>> from sympy import *
             >>> x,y = symbols('xy')
  
             >>> p = Poly(4*x**2 + 2*x, x)
             >>> p.as_primitive()
             (2, Poly(2*x**2 + x, x))
  
          """
          content = self.content
  
          if content is S.Zero or content is S.One:
              return content, self
          else:
              coeffs = [ coeff / content for coeff in self.coeffs ]
  
              return content, self.__class__((coeffs,
                  self.monoms), *self.symbols, **self.flags)
  
      def as_squarefree(self):
          """Returns square-free part of a polynomial.
  
             >>> from sympy import *
             >>> x,y = symbols('xy')
  
             >>> Poly((x-1)**2, x).as_squarefree()
             Poly(x - 1, x)
  
          """
          f, A = self, sympy.polys.algorithms
  
          if f.is_univariate:
              F = f.as_primitive()[1]
  
              h = f.diff()
  
              g = A.poly_gcd(F, h)
              r = A.poly_div(f, g)
  
              return r[0]
          else:
              raise MultivariatePolyError(f)
  
      def as_reduced(self):
          """Remove GCD of monomials from 'self'.
  
             >>> from sympy import *
             >>> x,y = symbols('xy')
  
             >>> Poly(x**3 + x, x).as_reduced()
             ((1,), Poly(x**2 + 1, x))
  
             >>> Poly(x**3*y+x**2*y**2, x, y).as_reduced()
             ((2, 1), Poly(x + y, x, y))
  
          """
          if self.is_inhomogeneous:
              return (0,)*len(self.symbols), self
          else:
              if self.is_univariate:
                  gcd = self.monoms[-1]
  
                  terms = self.coeffs, [ (monom - gcd[0],)
                      for (monom,) in self.monoms ]
              else:
                  gcd = monomial_min(*self.monoms)
  
                  if all(not n for n in gcd):
                      return gcd, self
  
                  terms = {}
  
                  for coeff, monom in self.iter_terms():
                      terms[monomial_div(monom, gcd)] = coeff
  
              return gcd, Poly(terms, *self.symbols, **self.flags)
  
      def map_coeffs(self, f, *args, **kwargs):
          """Apply a function to all the coefficients.
  
             >>> from sympy import *
             >>> x,y,u,v = symbols('xyuv')
  
             >>> p = Poly(x**2 + 2*x*y, x, y)
             >>> q = p.map_coeffs(lambda c: 2*c)
  
             >>> q.as_basic()
             4*x*y + 2*x**2
  
             >>> p = Poly(u*x**2 + v*x*y, x, y)
             >>> q = p.map_coeffs(expand, complex=True)
  
             >>> q.as_basic()
             x**2*(I*im(u) + re(u)) + x*y*(I*im(v) + re(v))
  
          """
          terms = []
  
          for coeff, monom in self.iter_terms():
              coeff = f(coeff, *args, **kwargs)
  
              if coeff.has_any_symbols(*self.symbols):
                  raise CoefficientError("%s coefficient is dependent" \
                      " of polynomial's symbols %s" % (coeff, self.symbols))
              elif coeff is not S.Zero:
                  terms.append((coeff, monom))
  
          return self.__class__(terms, *self.symbols, **self.flags)
  
      def coeff(self, *monom):
          """Returns coefficient of a particular monomial.
  
             >>> from sympy import *
             >>> x,y = symbols('xy')
  
             >>> p = Poly(3*x**2*y + 4*x*y**2 + 1, x, y)
  
             >>> p.coeff(2, 1)
             3
             >>> p.coeff(1, 2)
             4
             >>> p.coeff(1, 1)
             0
  
             If len(monom) == 0 then returns coeff of the constant term:
  
             >>> p.coeff(0, 0)
             1
             >>> p.coeff()
             1
  
          """
          if not monom:
              if all(e == 0 for e in self.monoms[-1]):
                  return self.coeffs[-1]
          else:
              for i in xrange(len(self.monoms)):
                  if self.monoms[i] == monom:
                      return self.coeffs[i]
  
          return S.Zero
  
      def unify_with(self, other):
          """Unify 'self' with a polynomial or a set of polynomials.
  
             This method will return polynomials of the same type,  dominated
             by 'self', with a common set of symbols (which is a union of all
             symbols from all polynomials) and with common monomial ordering.
  
             You can pass a polynomial  or an expression to this method, or
             a list or a tuple of polynomials or expressions. If any of the
             inputs would be an expression then it will be converted to a
             polynomial.
  
          """
          symbols, stamp = self.symbols, self.stamp
          flags, cls = self.flags, self.__class__
  
          if isinstance(other, (tuple, list, set)):
              for poly in other:
                  if isinstance(poly, Poly):
                      stamp |= poly.stamp
                  #elif atoms:
                  #    stamp |= poly.atoms(Symbol)
  
              if not (stamp <= self.stamp):
                  symbols = tuple(sorted(stamp))
                  self = cls(self, *symbols, **flags)
  
              other = other.__class__( cls(poly,
                  *symbols, **flags) for poly in other )
          else:
              other = sympify(other)
  
              if isinstance(other, Poly):
                  stamp |= other.stamp
              #elif atoms:
              #    stamp |= other.atoms(Symbol)
  
              if not (stamp <= self.stamp):
                  symbols = tuple(sorted(stamp))
                  self = cls(self, *symbols, **flags)
  
              other = cls(other, *symbols, **flags)
  
          return self, other
  
      def diff(self, *symbols):
          """Efficiently differentiate polynomials.
  
             Differentiate a polynomial with respect to a set of symbols.  If
             a symbol is polynomial's variable, then monomial differentiation
             is performed and coefficients updated with integer factors. In
             other case each of the coefficients is differentiated.
  
             Additionally, for each of the symbols you can specify a single
             positive integer which will indicate the number of times to
             perform differentiation using this symbol.
  
             >>> from sympy import *
             >>> x,y,z = symbols('xyz')
  
             >>> p = Poly(x**2*y + z**2, x, y)
  
             >>> p.diff(x)
             Poly(2*x*y, x, y)
  
             >>> p.diff(x, 2)
             Poly(2*y, x, y)
  
             >>> p.diff(x, 2, y)
             Poly(2, x, y)
  
             >>> p.diff(z)
             Poly(2*z, x, y)
  
          """
          if self.is_zero:
              return self
  
          new_symbols = []
  
          for s in symbols:
              s = sympify(s)
  
              if s.is_Symbol:
                  new_symbols.append((s, 1))
              elif s.is_Integer:
                  sym, _ = new_symbols.pop()
                  new_symbols.append((sym, int(s)))
              else:
                  raise TypeError
  
          if not new_symbols:
              if self.is_univariate:
                  new_symbols = [(self.symbols[0], 1)]
              else:
                  return self
  
          indices, symbols = {}, self.stamp
  
          for i, s in enumerate(self.symbols):
              indices[s] = i
  
          poly = self.as_dict()
  
          for s, k in new_symbols:
              new_poly = {}
  
              if s in symbols:
                  i = indices[s]
  
                  for M, coeff in poly.iteritems():
                      n = M[i] - k
  
                      if n >= 0:
                          monom = M[:i]+(n,)+M[i+1:]
  
                          for j in xrange(n, M[i]):
                              coeff *= j+1
  
                          if new_poly.has_key(monom):
                              coeff += new_poly[monom]
  
                              if not coeff:
                                  del new_poly[monom]
                                  continue
  
                          new_poly[monom] = coeff
              elif not isinstance(self, IntegerPoly):
                  for monom, coeff in poly.iteritems():
                      if coeff.has_any_symbols(s):
                          coeff = coeff.diff(*([s]*k))
  
                          if coeff:
                              new_poly[monom] = coeff
  
              poly = new_poly
  
              if not poly:
                  break
  
          return self.__class__(poly, *self.symbols, **self.flags)
  
      def integrate(self, *symbols):
          """Efficiently integrate polynomials.
  
             Integrate a polynomial with respect a set of symbols. If a
             symbol is polynomial's variable, then monomial integration
             is performed and coefficients updated with integer factors.
             In other case each of the coefficients is integrated.
  
             Additionally, for each of the symbols you can specify a
             single positive integer which will indicate the number
             of times to perform integration using this symbol.
  
             >>> from sympy import *
             >>> x,y,z = symbols('xyz')
  
             >>> p = Poly(x**2*y + z**2, x, y)
  
             >>> p.integrate(x)
             Poly(1/3*x**3*y + z**2*x, x, y)
  
             >>> p.integrate(x, 2)
             Poly(1/12*x**4*y + z**2/2*x**2, x, y)
  
             >>> p.integrate(x, 2, y)
             Poly(1/24*x**4*y**2 + z**2/2*x**2*y, x, y)
  
             >>> p.integrate(z)
             Poly(z*x**2*y + z**3/3, x, y)
  
          """
          if self.is_zero:
              return self
  
          new_symbols = []
  
          for s in symbols:
              s = sympify(s)
  
              if s.is_Symbol:
                  new_symbols.append((s, 1))
              elif s.is_Integer:
                  sym, _ = new_symbols.pop()
                  new_symbols.append((sym, int(s)))
              else:
                  raise TypeError
  
          indices, symbols = {}, self.stamp
  
          for i, s in enumerate(self.symbols):
              indices[s] = i
  
          poly = self.as_dict()
  
          for s, k in new_symbols:
              new_poly = {}
  
              if s in symbols:
                  i = indices[s]
  
                  for M, coeff in poly.iteritems():
                      n = M[i] + k
  
                      monom = M[:i]+(n,)+M[i+1:]
  
                      for j in xrange(M[i], n):
                          coeff /= j+1
  
                      if new_poly.has_key(monom):
                          coeff += new_poly[monom]
  
                          if not coeff:
                              del new_poly[monom]
                              continue
  
                      new_poly[monom] = coeff
              else:
                  for M, coeff in poly.iteritems():
                      new_poly[M] = C.Integral(coeff, *([s]*k)).doit()
  
              if not new_poly:
                  break
              else:
                  poly = new_poly
  
          return self.__class__(poly, *self.symbols, **self.flags)
  
      def invert(self, modulus):
          """Compute multiplicative inverse in K[x]. """
          self, modulus = self.unify_with(modulus)
  
          if self.is_multivariate:
              raise PolynomialError("only univariate polynomials can be inverted")
  
          s, g = sympy.polys.algorithms.poly_half_gcdex(self, modulus)
  
          if g.is_one:
              return sympy.polys.algorithms.poly_div(s, modulus)[1]
          else:
              raise ZeroDivisionError("polynomial division")
  
  
      def add_term(self, coeff, monom):
          """Efficiently add a single term to 'self'.
  
             The list of monomials is sorted at initialization time, this
             motivates usage of binary search algorithm to find an index
             of an existing monomial or a suitable place for a new one.
             This gives O(lg(n)) complexity, where 'n' is the initial
             number of terms, superior to naive approach.
  
             For more information on the implemented algorithm refer to:
  
             [1] D. E. Knuth, The Art of Computer Programming: Sorting
                 and Searching, v.1, Addison-Wesley Professional, 1998
  
          """
          if self.is_zero:
              coeffs = (coeff,)
              monoms = (monom,)
          else:
              coeffs = list(self.coeffs)
              monoms = list(self.monoms)
  
              compare = monomial_cmp(self.order)
  
              if compare(monom, monoms[0]) > 0:
                  coeffs.insert(0, coeff)
                  monoms.insert(0, monom)
              elif compare(monom, monoms[-1]) < 0:
                  coeffs.append(coeff)
                  monoms.append(monom)
              else:
                  lo, hi = 0, len(monoms)-1
  
                  while lo <= hi:
                      i = (lo + hi) // 2
  
                      k = compare(monom, monoms[i])
  
                      if not k:
                          coeff += coeffs[i]
  
                          if coeff:
                              coeffs[i] = coeff
                          else:
                              del coeffs[i]
                              del monoms[i]
  
                          break
                      else:
                          if k > 0:
                              hi = i - 1
                          else:
                              lo = i + 1
                  else:
                      coeffs.insert(i, coeff)
                      monoms.insert(i, monom)
  
          return self.__class__((coeffs, monoms),
              *self.symbols, **self.flags)
  
      def sub_term(self, coeff, monom):
          """Efficiently subtract a single term from 'self'.
  
             The list of monomials is sorted at initialization time, this
             motivates usage of binary search algorithm to find an index
             of an existing monomial or a suitable place for a new one.
             This gives O(lg(n)) complexity, where 'n' is the initial
             number of terms, superior to naive approach.
  
             For more information on the implemented algorithm refer to:
  
             [1] D. E. Knuth, The Art of Computer Programming: Sorting
                 and Searching, v.2, Addison-Wesley Professional, 1998
  
          """
          if self.is_zero:
              coeffs = (-coeff,)
              monoms = (monom,)
          else:
              coeffs = list(self.coeffs)
              monoms = list(self.monoms)
  
              compare = monomial_cmp(self.order)
  
              if compare(monom, monoms[0]) > 0:
                  coeffs.insert(0, -coeff)
                  monoms.insert(0, monom)
              elif compare(monom, monoms[-1]) < 0:
                  coeffs.append(-coeff)
                  monoms.append(monom)
              else:
                  lo, hi = 0, len(monoms)-1
  
                  while lo <= hi:
                      i = (lo + hi) // 2
  
                      k = compare(monom, monoms[i])
  
                      if not k:
                          coeff -= coeffs[i]
  
                          if coeff:
                              coeffs[i] = -coeff
                          else:
                              del coeffs[i]
                              del monoms[i]
  
                          break
                      else:
                          if k > 0:
                              hi = i - 1
                          else:
                              lo = i + 1
                  else:
                      coeffs.insert(i, -coeff)
                      monoms.insert(i, monom)
  
          return self.__class__((coeffs, monoms),
              *self.symbols, **self.flags)
  
      def mul_term(self, coeff, monom=None):
          """Efficiently multiply 'self' with a single term. """
          coeff = sympify(coeff)
  
          if coeff is S.Zero:
              terms = ()
          else:
              if coeff is S.One:
                  if monom is None:
                      return self
                  else:
                      coeffs = self.coeffs
              else:
                  coeffs = [ p_coeff * coeff for p_coeff in self.coeffs ]
  
                  for i, coeff in enumerate(coeffs):
                      coeffs[i] = Poly.cancel(coeff)
  
              if monom is None:
                  monoms = self.monoms
              else:
                  monoms = [ monomial_mul(p_monom, monom)
                      for p_monom in self.monoms ]
  
              terms = coeffs, monoms
  
          return self.__class__(terms, *self.symbols, **self.flags)
  
      def div_term(self, coeff, monom=None):
          """Efficiently divide 'self' by a single term. """
          coeff = sympify(coeff)
  
          if coeff is S.Zero:
              raise ZeroDivisionError
  
          if coeff is S.One:
              if monom is None:
                  return self
              else:
                  coeffs = self.coeffs
          else:
              coeffs = [ p_coeff / coeff for p_coeff in self.coeffs ]
  
              for i, coeff in enumerate(coeffs):
                  coeffs[i] = Poly.cancel(coeff)
  
          if monom is None:
              monoms = self.monoms
          else:
              monoms = [ monomial_div(p_monom, monom)
                  for p_monom in self.monoms ]
  
              if any(monom is None for monom in monoms):
                  raise PolynomialError("%s monomial must divide" \
                      " exactly the given polynomial" % (monom,))
  
          return self.__class__((coeffs, monoms),
              *self.symbols, **self.flags)
  
      def kill_lead_term(self):
          """Removes leading term from 'self'. """
          terms = self.coeffs[1:], self.monoms[1:]
          return self.__class__(terms, *self.symbols, **self.flags)
  
      def kill_last_term(self):
          """Removes last term from 'self'. """
          terms = self.coeffs[:-1], self.monoms[:-1]
          return self.__class__(terms, *self.symbols, **self.flags)
  
      def __call__(self, *point):
          """Efficiently evaluate polynomial at a given point.
  
             Evaluation is always done using Horner scheme.  In multivariate
             case a greedy algorithm is used to obtain a sequence of partial
             evaluations which minimizes the total number of multiplications
             required to perform this evaluation. This strategy is efficient
             for most of multivariate polynomials.
  
             Note that evaluation is done for all variables, which means the
             dimension of the given point must match the number of symbols.
  
             If you wish to efficiently evaluate polynomial for a subset of
             symbols use 'evaluate' method instead. Alternatively one can
             use Basic.subs() for this purpose.
  
             >>> from sympy import *
             >>> x,y = symbols('xy')
  
             >>> p = Poly(x**2 + 2*x + 1, x)
  
             >>> p(2)
             9
             >>> p(y)
             1 + y*(2 + y)
  
             For more information on the implemented algorithm refer to:
  
             [1] M. Ceberio, V. Kreinovich, Greedy Algorithms for Optimizing
                 Multivariate Horner Schemes, ACM SIGSAM Bulletin, Volume 38,
                 Issue 1, 2004, pp. 8-15
  
          """
          N = len(self.symbols)
  
          if len(point) != N:
              raise ValueError("Dimension of %s does not match" \
                  " the number of symbols (%d)" % (point, N))
  
          if self.is_univariate:
              terms = self.as_uv_dict()
  
              point, result = point[0], 0
  
              for k in xrange(self.degree, -1, -1):
                  result *= point
  
                  if terms.has_key(k):
                      result += terms[k]
  
              return result
          else:
              def evaluate(terms):
                  count = [0] * N
  
                  for monom in terms.iterkeys():
                      for i, M in enumerate(monom):
                          if M != 0:
                              count[i] += 1
  
                  K = max(count)
  
                  if K <= 1:
                      result = 0
  
                      for monom, coeff in terms.iteritems():
                          for base, exp in zip(point, monom):
                              if exp != 0:
                                  if exp == 1:
                                      coeff *= base
                                  else:
                                      coeff *= base**exp
  
                          result += coeff
  
                      return result
                  else:
                      k, indeps, depend = count.index(K), {}, {}
  
                      n = min([ M[k] for M in terms.iterkeys() if M[k] ])
  
                      for M, coeff in terms.iteritems():
                          if M[k] != 0:
                              depend[M[:k]+(M[k]-n,)+M[k+1:]] = coeff
                          else:
                              indeps[M] = coeff
  
                      result = point[k]**n * evaluate(depend)
  
                      if indeps:
                          return result + evaluate(indeps)
                      else:
                          return result
  
              return evaluate(self.as_dict())
  
      def evaluate(self, pattern):
          """Evaluates polynomial for a given set of symbols. """
          symbols = list(self.symbols)
  
          if type(pattern) is dict:
              pattern = pattern.items()
          elif type(pattern) is tuple:
              pattern = [pattern]
  
          poly = self.as_dict()
  
          for s, value in pattern:
              if s not in self.stamp:
                  raise PolynomialError("%s symbol does" \
                      " not belong to %s" % (s, self.symbols))
              else:
                  i = symbols.index(s)
  
              # 'value' might be int | long
              if isinstance(value, Symbol):
                  if value in self.stamp:
                      raise PolynomialError("%s symbol must" \
                          " not belong to %s" % (s, self.symbols))
  
              terms = {}
  
              for M, coeff in poly.iteritems():
                  monom = M[:i] + M[i+1:]
                  coeff *= value ** M[i]
  
                  coeff = Poly.cancel(coeff)
  
                  if terms.has_key(monom):
                      coeff += terms[monom]
  
                      if not coeff:
                          del terms[monom]
                          continue
                  elif not coeff:
                      continue
  
                  terms[monom] = coeff
  
              del symbols[i]
              poly = terms
  
          if not symbols:
              if not poly:
                  return S.Zero
              else:
                  return poly.popitem()[1]
          else:
              return self.__class__(poly, *symbols, **self.flags)
  
      def _eval_subs(self, old, new):
          symbols = list(self.symbols)
  
          if old in self.stamp:
              terms, i = {}, symbols.index(old)
  
              if new.is_Symbol:
                  if new in self.stamp:
                      j = symbols.index(new)
  
                      for coeff, M in self.iter_terms():
                          n, m = M[i], M[j]
  
                          if i < j:
                              monom = M[:i]+M[i+1:j]+(m+n,)+M[j+1:]
                          else:
                              monom = M[:j]+(m+n,)+M[j+1:i]+M[i+1:]
  
                          if terms.has_key(monom):
                              coeff += terms[monom]
  
                              if not coeff:
                                  del terms[monom]
                                  continue
  
                          terms[monom] = coeff
  
                      del symbols[i]
  
                      return self.__class__(terms, *symbols, **self.flags)
                  else:
                      for coeff in self.coeffs:
                          if coeff.has_any_symbols(new):
                              break
                      else:
                          symbols[i], terms = new, (self.coeffs, self.monoms)
                          return self.__class__(terms, *symbols, **self.flags)
              elif new.is_number:
                  if len(symbols) == 1:
                      return self(new)
                  else:
                      return self.evaluate((old, new))
          elif not new.has_any_symbols(*symbols):
              coeffs = [ sympify(coeff)._eval_subs(old, new) for coeff in self.coeffs ]
              return self.__class__(zip(coeffs, self.monoms), *symbols, **self.flags)
  
          result = self.as_basic()._eval_subs(old, new)
  
          try:
              return self.__class__(result, *symbols, **self.flags)
          except PolynomialError:
              return result
  
      def __eq__(self, other):
          """Compare polynomials up to order of symbols and monomials. """
          try:
              poly = self.__class__(other, *self.symbols, **self.flags)
          except PolynomialError:
              return False
  
          if self.length != poly.length:
              return False
  
          if hash(self.monoms) != hash(poly.monoms):
              return False
  
          if hash(self.coeffs) != hash(poly.coeffs):
              for a, b in zip(self.coeffs, poly.coeffs):
                  if a != b:
                      return False
  
          return True
  
      def __ne__(self, other):
          """Compare polynomials up to order of symbols and monomials. """
          try:
              poly = self.__class__(other, *self.symbols, **self.flags)
          except PolynomialError:
              return True
  
          if self.length != poly.length:
              return True
  
          if hash(self.monoms) != hash(poly.monoms):
              return True
  
          if hash(self.coeffs) != hash(poly.coeffs):
              for a, b in zip(self.coeffs, poly.coeffs):
                  if a != b:
                      return True
  
          return False
  
      def __nonzero__(self):
          return self.coeffs not in ((S.Zero,), (0,))
  
      def _eval_is_polynomial(self, symbols):
          try:
              self.__class__(self, *symbols, **self.flags)
          except PolynomialError:
              return False
  
          return True
  
  class IntegerPoly(Poly):
      pass
  
  def multinomial_as_basic(multinomial, *symbols):
      """
      Converts the multinomial to Add/Mul/Pow instances.
  
      multinomial is a dict of {powers: coefficient} pairs, powers is a tuple of
      python integers, coefficient is a python integer.
      """
      l = []
      for powers, k in multinomial.iteritems():
          term = [k]
          for i,e in enumerate(powers):
              term.append(Pow(symbols[i], powers[i]))
          l.append(Mul(*term))
      result = Add(*l)
      return result
  
  