text: """User-friendly public interface to polynomial functions. """
  
  from __future__ import print_function, division
  
  from functools import wraps, reduce
  from operator import mul
  
  from sympy.core import (
      S, Basic, Expr, I, Integer, Add, Mul, Dummy, Tuple
  )
  from sympy.core.basic import preorder_traversal
  from sympy.core.compatibility import iterable, ordered
  from sympy.core.decorators import _sympifyit
  from sympy.core.function import Derivative
  from sympy.core.mul import _keep_coeff
  from sympy.core.relational import Relational
  from sympy.core.symbol import Symbol
  from sympy.core.sympify import sympify, _sympify
  from sympy.logic.boolalg import BooleanAtom
  from sympy.polys import polyoptions as options
  from sympy.polys.constructor import construct_domain
  from sympy.polys.domains import FF, QQ, ZZ
  from sympy.polys.fglmtools import matrix_fglm
  from sympy.polys.groebnertools import groebner as _groebner
  from sympy.polys.monomials import Monomial
  from sympy.polys.orderings import monomial_key
  from sympy.polys.polyclasses import DMP
  from sympy.polys.polyerrors import (
      OperationNotSupported, DomainError,
      CoercionFailed, UnificationFailed,
      GeneratorsNeeded, PolynomialError,
      MultivariatePolynomialError,
      ExactQuotientFailed,
      PolificationFailed,
      ComputationFailed,
      GeneratorsError,
  )
  from sympy.polys.polyutils import (
      basic_from_dict,
      _sort_gens,
      _unify_gens,
      _dict_reorder,
      _dict_from_expr,
      _parallel_dict_from_expr,
  )
  from sympy.polys.rationaltools import together
  from sympy.polys.rootisolation import dup_isolate_real_roots_list
  from sympy.utilities import group, sift, public, filldedent
  from sympy.utilities.exceptions import SymPyDeprecationWarning
  
  # Required to avoid errors
  import sympy.polys
  
  import mpmath
  from mpmath.libmp.libhyper import NoConvergence
  from sympy.abc import x, y,z#for issue 6322
  from sympy.parsing.sympy_parser import parse_expr
  
  
  
  def _polifyit(func):
      @wraps(func)
      def wrapper(f, g):
          g = _sympify(g)
          if isinstance(g, Poly):
              return func(f, g)
          elif isinstance(g, Expr):
              try:
                  g = f.from_expr(g, *f.gens)
              except PolynomialError:
                  if g.is_Matrix:
                      return NotImplemented
                  expr_method = getattr(f.as_expr(), func.__name__)
                  result = expr_method(g)
                  if result is not NotImplemented:
                      SymPyDeprecationWarning(
                          feature="Mixing Poly with non-polynomial expressions in binary operations",
                          issue=18613,
                          deprecated_since_version="1.6",
                          useinstead="the as_expr or as_poly method to convert types").warn()
                  return result
              else:
                  return func(f, g)
          else:
              return NotImplemented
      return wrapper
  
  
  
  @public
  class Poly(Basic):
      """
      Generic class for representing and operating on polynomial expressions.
  
      Poly is a subclass of Basic rather than Expr but instances can be
      converted to Expr with the ``as_expr`` method.
  
      Examples
      ========
  
      >>> from sympy import Poly
      >>> from sympy.abc import x, y
  
      Create a univariate polynomial:
  
      >>> Poly(x*(x**2 + x - 1)**2)
      Poly(x**5 + 2*x**4 - x**3 - 2*x**2 + x, x, domain='ZZ')
  
      Create a univariate polynomial with specific domain:
  
      >>> from sympy import sqrt
      >>> Poly(x**2 + 2*x + sqrt(3), domain='R')
      Poly(1.0*x**2 + 2.0*x + 1.73205080756888, x, domain='RR')
  
      Create a multivariate polynomial:
  
      >>> Poly(y*x**2 + x*y + 1)
      Poly(x**2*y + x*y + 1, x, y, domain='ZZ')
  
      Create a univariate polynomial, where y is a constant:
  
      >>> Poly(y*x**2 + x*y + 1,x)
      Poly(y*x**2 + y*x + 1, x, domain='ZZ[y]')
  
      You can evaluate the above polynomial as a function of y:
  
      >>> Poly(y*x**2 + x*y + 1,x).eval(2)
      6*y + 1
  
      See Also
      ========
  
      sympy.core.expr.Expr
  
      """
  
      __slots__ = ('rep', 'gens')
  
      is_commutative = True
      is_Poly = True
      _op_priority = 10.001
  
      def __new__(cls, rep, *gens, **args):
          """Create a new polynomial instance out of something useful. """
          opt = options.build_options(gens, args)
  
          if 'order' in opt:
              raise NotImplementedError("'order' keyword is not implemented yet")
  
          if iterable(rep, exclude=str):
              if isinstance(rep, dict):
                  return cls._from_dict(rep, opt)
              else:
                  return cls._from_list(list(rep), opt)
          else:
              rep = sympify(rep)
  
              if rep.is_Poly:
                  return cls._from_poly(rep, opt)
              else:
                  return cls._from_expr(rep, opt)
  
      # Poly does not pass its args to Basic.__new__ to be stored in _args so we
      # have to emulate them here with an args property that derives from rep
      # and gens which are instance attributes. This also means we need to
      # define _hashable_content. The _hashable_content is rep and gens but args
      # uses expr instead of rep (expr is the Basic version of rep). Passing
      # expr in args means that Basic methods like subs should work. Using rep
      # otherwise means that Poly can remain more efficient than Basic by
      # avoiding creating a Basic instance just to be hashable.
  
      @classmethod
      def new(cls, rep, *gens):
          """Construct :class:`Poly` instance from raw representation. """
          if not isinstance(rep, DMP):
              raise PolynomialError(
                  "invalid polynomial representation: %s" % rep)
          elif rep.lev != len(gens) - 1:
              raise PolynomialError("invalid arguments: %s, %s" % (rep, gens))
  
          obj = Basic.__new__(cls)
          obj.rep = rep
          obj.gens = gens
  
          return obj
  
      @property
      def expr(self):
          return basic_from_dict(self.rep.to_sympy_dict(), *self.gens)
  
      @property
      def args(self):
          return (self.expr,) + self.gens
  
      def _hashable_content(self):
          return (self.rep,) + self.gens
  
      @classmethod
      def from_dict(cls, rep, *gens, **args):
          """Construct a polynomial from a ``dict``. """
          opt = options.build_options(gens, args)
          return cls._from_dict(rep, opt)
  
      @classmethod
      def from_list(cls, rep, *gens, **args):
          """Construct a polynomial from a ``list``. """
          opt = options.build_options(gens, args)
          return cls._from_list(rep, opt)
  
      @classmethod
      def from_poly(cls, rep, *gens, **args):
          """Construct a polynomial from a polynomial. """
          opt = options.build_options(gens, args)
          return cls._from_poly(rep, opt)
  
      @classmethod
      def from_expr(cls, rep, *gens, **args):
          """Construct a polynomial from an expression. """
          opt = options.build_options(gens, args)
          return cls._from_expr(rep, opt)
  
      @classmethod
      def _from_dict(cls, rep, opt):
          """Construct a polynomial from a ``dict``. """
          gens = opt.gens
  
          if not gens:
              raise GeneratorsNeeded(
                  "can't initialize from 'dict' without generators")
  
          level = len(gens) - 1
          domain = opt.domain
  
          if domain is None:
              domain, rep = construct_domain(rep, opt=opt)
          else:
              for monom, coeff in rep.items():
                  rep[monom] = domain.convert(coeff)
  
          return cls.new(DMP.from_dict(rep, level, domain), *gens)
  
      @classmethod
      def _from_list(cls, rep, opt):
          """Construct a polynomial from a ``list``. """
          gens = opt.gens
  
          if not gens:
              raise GeneratorsNeeded(
                  "can't initialize from 'list' without generators")
          elif len(gens) != 1:
              raise MultivariatePolynomialError(
                  "'list' representation not supported")
  
          level = len(gens) - 1
          domain = opt.domain
  
          if domain is None:
              domain, rep = construct_domain(rep, opt=opt)
          else:
              rep = list(map(domain.convert, rep))
  
          return cls.new(DMP.from_list(rep, level, domain), *gens)
  
      @classmethod
      def _from_poly(cls, rep, opt):
          """Construct a polynomial from a polynomial. """
          if cls != rep.__class__:
              rep = cls.new(rep.rep, *rep.gens)
  
          gens = opt.gens
          field = opt.field
          domain = opt.domain
  
          if gens and rep.gens != gens:
              if set(rep.gens) != set(gens):
                  return cls._from_expr(rep.as_expr(), opt)
              else:
                  rep = rep.reorder(*gens)
  
          if 'domain' in opt and domain:
              rep = rep.set_domain(domain)
          elif field is True:
              rep = rep.to_field()
  
          return rep
  
      @classmethod
      def _from_expr(cls, rep, opt):
          """Construct a polynomial from an expression. """
          rep, opt = _dict_from_expr(rep, opt)
          return cls._from_dict(rep, opt)
  
      def __hash__(self):
          return super(Poly, self).__hash__()
  
      @property
      def free_symbols(self):
          """
          Free symbols of a polynomial expression.
  
          Examples
          ========
  
          >>> from sympy import Poly
          >>> from sympy.abc import x, y, z
  
          >>> Poly(x**2 + 1).free_symbols
          {x}
          >>> Poly(x**2 + y).free_symbols
          {x, y}
          >>> Poly(x**2 + y, x).free_symbols
          {x, y}
          >>> Poly(x**2 + y, x, z).free_symbols
          {x, y}
  
          """
          symbols = set()
          gens = self.gens
          for i in range(len(gens)):
              for monom in self.monoms():
                  if monom[i]:
                      symbols |= gens[i].free_symbols
                      break
  
          return symbols | self.free_symbols_in_domain
  
      @property
      def free_symbols_in_domain(self):
          """
          Free symbols of the domain of ``self``.
  
          Examples
          ========
  
          >>> from sympy import Poly
          >>> from sympy.abc import x, y
  
          >>> Poly(x**2 + 1).free_symbols_in_domain
          set()
          >>> Poly(x**2 + y).free_symbols_in_domain
          set()
          >>> Poly(x**2 + y, x).free_symbols_in_domain
          {y}
  
          """
          domain, symbols = self.rep.dom, set()
  
          if domain.is_Composite:
              for gen in domain.symbols:
                  symbols |= gen.free_symbols
          elif domain.is_EX:
              for coeff in self.coeffs():
                  symbols |= coeff.free_symbols
  
          return symbols
  
      @property
      def gen(self):
          """
          Return the principal generator.
  
          Examples
          ========
  
          >>> from sympy import Poly
          >>> from sympy.abc import x
  
          >>> Poly(x**2 + 1, x).gen
          x
  
          """
          return self.gens[0]
  
      @property
      def domain(self):
          """Get the ground domain of ``self``. """
          return self.get_domain()
  
      @property
      def zero(self):
          """Return zero polynomial with ``self``'s properties. """
          return self.new(self.rep.zero(self.rep.lev, self.rep.dom), *self.gens)
  
      @property
      def one(self):
          """Return one polynomial with ``self``'s properties. """
          return self.new(self.rep.one(self.rep.lev, self.rep.dom), *self.gens)
  
      @property
      def unit(self):
          """Return unit polynomial with ``self``'s properties. """
          return self.new(self.rep.unit(self.rep.lev, self.rep.dom), *self.gens)
  
      def unify(f, g):
          """
          Make ``f`` and ``g`` belong to the same domain.
  
          Examples
          ========
  
          >>> from sympy import Poly
          >>> from sympy.abc import x
  
          >>> f, g = Poly(x/2 + 1), Poly(2*x + 1)
  
          >>> f
          Poly(1/2*x + 1, x, domain='QQ')
          >>> g
          Poly(2*x + 1, x, domain='ZZ')
  
          >>> F, G = f.unify(g)
  
          >>> F
          Poly(1/2*x + 1, x, domain='QQ')
          >>> G
          Poly(2*x + 1, x, domain='QQ')
  
          """
          _, per, F, G = f._unify(g)
          return per(F), per(G)
  
      def _unify(f, g):
          g = sympify(g)
  
          if not g.is_Poly:
              try:
                  return f.rep.dom, f.per, f.rep, f.rep.per(f.rep.dom.from_sympy(g))
              except CoercionFailed:
                  raise UnificationFailed("can't unify %s with %s" % (f, g))
  
          if isinstance(f.rep, DMP) and isinstance(g.rep, DMP):
              gens = _unify_gens(f.gens, g.gens)
  
              dom, lev = f.rep.dom.unify(g.rep.dom, gens), len(gens) - 1
  
              if f.gens != gens:
                  f_monoms, f_coeffs = _dict_reorder(
                      f.rep.to_dict(), f.gens, gens)
  
                  if f.rep.dom != dom:
                      f_coeffs = [dom.convert(c, f.rep.dom) for c in f_coeffs]
  
                  F = DMP(dict(list(zip(f_monoms, f_coeffs))), dom, lev)
              else:
                  F = f.rep.convert(dom)
  
              if g.gens != gens:
                  g_monoms, g_coeffs = _dict_reorder(
                      g.rep.to_dict(), g.gens, gens)
  
                  if g.rep.dom != dom:
                      g_coeffs = [dom.convert(c, g.rep.dom) for c in g_coeffs]
  
                  G = DMP(dict(list(zip(g_monoms, g_coeffs))), dom, lev)
              else:
                  G = g.rep.convert(dom)
          else:
              raise UnificationFailed("can't unify %s with %s" % (f, g))
  
          cls = f.__class__
  
          def per(rep, dom=dom, gens=gens, remove=None):
              if remove is not None:
                  gens = gens[:remove] + gens[remove + 1:]
  
                  if not gens:
                      return dom.to_sympy(rep)
  
              return cls.new(rep, *gens)
  
          return dom, per, F, G
  
      def per(f, rep, gens=None, remove=None):
          """
          Create a Poly out of the given representation.
  
          Examples
          ========
  
          >>> from sympy import Poly, ZZ
          >>> from sympy.abc import x, y
  
          >>> from sympy.polys.polyclasses import DMP
  
          >>> a = Poly(x**2 + 1)
  
          >>> a.per(DMP([ZZ(1), ZZ(1)], ZZ), gens=[y])
          Poly(y + 1, y, domain='ZZ')
  
          """
          if gens is None:
              gens = f.gens
  
          if remove is not None:
              gens = gens[:remove] + gens[remove + 1:]
  
              if not gens:
                  return f.rep.dom.to_sympy(rep)
  
          return f.__class__.new(rep, *gens)
  
      def set_domain(f, domain):
          """Set the ground domain of ``f``. """
          opt = options.build_options(f.gens, {'domain': domain})
          return f.per(f.rep.convert(opt.domain))
  
      def get_domain(f):
          """Get the ground domain of ``f``. """
          return f.rep.dom
  
      def set_modulus(f, modulus):
          """
          Set the modulus of ``f``.
  
          Examples
          ========
  
          >>> from sympy import Poly
          >>> from sympy.abc import x
  
          >>> Poly(5*x**2 + 2*x - 1, x).set_modulus(2)
          Poly(x**2 + 1, x, modulus=2)
  
          """
          modulus = options.Modulus.preprocess(modulus)
          return f.set_domain(FF(modulus))
  
      def get_modulus(f):
          """
          Get the modulus of ``f``.
  
          Examples
          ========
  
          >>> from sympy import Poly
          >>> from sympy.abc import x
  
          >>> Poly(x**2 + 1, modulus=2).get_modulus()
          2
  
          """
          domain = f.get_domain()
  
          if domain.is_FiniteField:
              return Integer(domain.characteristic())
          else:
              raise PolynomialError("not a polynomial over a Galois field")
  
      def _eval_subs(f, old, new):
          """Internal implementation of :func:`subs`. """
          if old in f.gens:
              if new.is_number:
                  return f.eval(old, new)
              else:
                  try:
                      return f.replace(old, new)
                  except PolynomialError:
                      pass
  
          return f.as_expr().subs(old, new)
  
      def exclude(f):
          """
          Remove unnecessary generators from ``f``.
  
          Examples
          ========
  
          >>> from sympy import Poly
          >>> from sympy.abc import a, b, c, d, x
  
          >>> Poly(a + x, a, b, c, d, x).exclude()
          Poly(a + x, a, x, domain='ZZ')
  
          """
          J, new = f.rep.exclude()
          gens = []
  
          for j in range(len(f.gens)):
              if j not in J:
                  gens.append(f.gens[j])
  
          return f.per(new, gens=gens)
  
      def replace(f, x, y=None, *_ignore):
          # XXX this does not match Basic's signature
          """
          Replace ``x`` with ``y`` in generators list.
  
          Examples
          ========
  
          >>> from sympy import Poly
          >>> from sympy.abc import x, y
  
          >>> Poly(x**2 + 1, x).replace(x, y)
          Poly(y**2 + 1, y, domain='ZZ')
  
          """
          if y is None:
              if f.is_univariate:
                  x, y = f.gen, x
              else:
                  raise PolynomialError(
                      "syntax supported only in univariate case")
  
          if x == y or x not in f.gens:
              return f
  
          if x in f.gens and y not in f.gens:
              dom = f.get_domain()
  
              if not dom.is_Composite or y not in dom.symbols:
                  gens = list(f.gens)
                  gens[gens.index(x)] = y
                  return f.per(f.rep, gens=gens)
  
          raise PolynomialError("can't replace %s with %s in %s" % (x, y, f))
  
      def match(f, *args, **kwargs):
          """Match expression from Poly. See Basic.match()"""
          return f.as_expr().match(*args, **kwargs)
  
      def reorder(f, *gens, **args):
          """
          Efficiently apply new order of generators.
  
          Examples
          ========
  
          >>> from sympy import Poly
          >>> from sympy.abc import x, y
  
          >>> Poly(x**2 + x*y**2, x, y).reorder(y, x)
          Poly(y**2*x + x**2, y, x, domain='ZZ')
  
          """
          opt = options.Options((), args)
  
          if not gens:
              gens = _sort_gens(f.gens, opt=opt)
          elif set(f.gens) != set(gens):
              raise PolynomialError(
                  "generators list can differ only up to order of elements")
  
          rep = dict(list(zip(*_dict_reorder(f.rep.to_dict(), f.gens, gens))))
  
          return f.per(DMP(rep, f.rep.dom, len(gens) - 1), gens=gens)
  
      def ltrim(f, gen):
          """
          Remove dummy generators from ``f`` that are to the left of
          specified ``gen`` in the generators as ordered. When ``gen``
          is an integer, it refers to the generator located at that
          position within the tuple of generators of ``f``.
  
          Examples
          ========
  
          >>> from sympy import Poly
          >>> from sympy.abc import x, y, z
  
          >>> Poly(y**2 + y*z**2, x, y, z).ltrim(y)
          Poly(y**2 + y*z**2, y, z, domain='ZZ')
          >>> Poly(z, x, y, z).ltrim(-1)
          Poly(z, z, domain='ZZ')
  
          """
          rep = f.as_dict(native=True)
          j = f._gen_to_level(gen)
  
          terms = {}
  
          for monom, coeff in rep.items():
  
              if any(monom[:j]):
                  # some generator is used in the portion to be trimmed
                  raise PolynomialError("can't left trim %s" % f)
  
              terms[monom[j:]] = coeff
  
          gens = f.gens[j:]
  
          return f.new(DMP.from_dict(terms, len(gens) - 1, f.rep.dom), *gens)
  
      def has_only_gens(f, *gens):
          """
          Return ``True`` if ``Poly(f, *gens)`` retains ground domain.
  
          Examples
          ========
  
          >>> from sympy import Poly
          >>> from sympy.abc import x, y, z
  
          >>> Poly(x*y + 1, x, y, z).has_only_gens(x, y)
          True
          >>> Poly(x*y + z, x, y, z).has_only_gens(x, y)
          False
  
          """
          indices = set()
  
          for gen in gens:
              try:
                  index = f.gens.index(gen)
              except ValueError:
                  raise GeneratorsError(
                      "%s doesn't have %s as generator" % (f, gen))
              else:
                  indices.add(index)
  
          for monom in f.monoms():
              for i, elt in enumerate(monom):
                  if i not in indices and elt:
                      return False
  
          return True
  
      def to_ring(f):
          """
          Make the ground domain a ring.
  
          Examples
          ========
  
          >>> from sympy import Poly, QQ
          >>> from sympy.abc import x
  
          >>> Poly(x**2 + 1, domain=QQ).to_ring()
          Poly(x**2 + 1, x, domain='ZZ')
  
          """
          if hasattr(f.rep, 'to_ring'):
              result = f.rep.to_ring()
          else:  # pragma: no cover
              raise OperationNotSupported(f, 'to_ring')
  
          return f.per(result)
  
      def to_field(f):
          """
          Make the ground domain a field.
  
          Examples
          ========
  
          >>> from sympy import Poly, ZZ
          >>> from sympy.abc import x
  
          >>> Poly(x**2 + 1, x, domain=ZZ).to_field()
          Poly(x**2 + 1, x, domain='QQ')
  
          """
          if hasattr(f.rep, 'to_field'):
              result = f.rep.to_field()
          else:  # pragma: no cover
              raise OperationNotSupported(f, 'to_field')
  
          return f.per(result)
  
      def to_exact(f):
          """
          Make the ground domain exact.
  
          Examples
          ========
  
          >>> from sympy import Poly, RR
          >>> from sympy.abc import x
  
          >>> Poly(x**2 + 1.0, x, domain=RR).to_exact()
          Poly(x**2 + 1, x, domain='QQ')
  
          """
          if hasattr(f.rep, 'to_exact'):
              result = f.rep.to_exact()
          else:  # pragma: no cover
              raise OperationNotSupported(f, 'to_exact')
  
          return f.per(result)
  
      def retract(f, field=None):
          """
          Recalculate the ground domain of a polynomial.
  
          Examples
          ========
  
          >>> from sympy import Poly
          >>> from sympy.abc import x
  
          >>> f = Poly(x**2 + 1, x, domain='QQ[y]')
          >>> f
          Poly(x**2 + 1, x, domain='QQ[y]')
  
          >>> f.retract()
          Poly(x**2 + 1, x, domain='ZZ')
          >>> f.retract(field=True)
          Poly(x**2 + 1, x, domain='QQ')
  
          """
          dom, rep = construct_domain(f.as_dict(zero=True),
              field=field, composite=f.domain.is_Composite or None)
          return f.from_dict(rep, f.gens, domain=dom)
  
      def slice(f, x, m, n=None):
          """Take a continuous subsequence of terms of ``f``. """
          if n is None:
              j, m, n = 0, x, m
          else:
              j = f._gen_to_level(x)
  
          m, n = int(m), int(n)
  
          if hasattr(f.rep, 'slice'):
              result = f.rep.slice(m, n, j)
          else:  # pragma: no cover
              raise OperationNotSupported(f, 'slice')
  
          return f.per(result)
  
      def coeffs(f, order=None):
          """
          Returns all non-zero coefficients from ``f`` in lex order.
  
          Examples
          ========
  
          >>> from sympy import Poly
          >>> from sympy.abc import x
  
          >>> Poly(x**3 + 2*x + 3, x).coeffs()
          [1, 2, 3]
  
          See Also
          ========
          all_coeffs
          coeff_monomial
          nth
  
          """
          return [f.rep.dom.to_sympy(c) for c in f.rep.coeffs(order=order)]
  
      def monoms(f, order=None):
          """
          Returns all non-zero monomials from ``f`` in lex order.
  
          Examples
          ========
  
          >>> from sympy import Poly
          >>> from sympy.abc import x, y
  
          >>> Poly(x**2 + 2*x*y**2 + x*y + 3*y, x, y).monoms()
          [(2, 0), (1, 2), (1, 1), (0, 1)]
  
          See Also
          ========
          all_monoms
  
          """
          return f.rep.monoms(order=order)
  
      def terms(f, order=None):
          """
          Returns all non-zero terms from ``f`` in lex order.
  
          Examples
          ========
  
          >>> from sympy import Poly
          >>> from sympy.abc import x, y
  
          >>> Poly(x**2 + 2*x*y**2 + x*y + 3*y, x, y).terms()
          [((2, 0), 1), ((1, 2), 2), ((1, 1), 1), ((0, 1), 3)]
  
          See Also
          ========
          all_terms
  
          """
          return [(m, f.rep.dom.to_sympy(c)) for m, c in f.rep.terms(order=order)]
  
      def all_coeffs(f):
          """
          Returns all coefficients from a univariate polynomial ``f``.
  
          Examples
          ========
  
          >>> from sympy import Poly
          >>> from sympy.abc import x
  
          >>> Poly(x**3 + 2*x - 1, x).all_coeffs()
          [1, 0, 2, -1]
  
          """
          return [f.rep.dom.to_sympy(c) for c in f.rep.all_coeffs()]
  
      def all_monoms(f):
          """
          Returns all monomials from a univariate polynomial ``f``.
  
          Examples
          ========
  
          >>> from sympy import Poly
          >>> from sympy.abc import x
  
          >>> Poly(x**3 + 2*x - 1, x).all_monoms()
          [(3,), (2,), (1,), (0,)]
  
          See Also
          ========
          all_terms
  
          """
          return f.rep.all_monoms()
  
      def all_terms(f):
          """
          Returns all terms from a univariate polynomial ``f``.
  
          Examples
          ========
  
          >>> from sympy import Poly
          >>> from sympy.abc import x
  
          >>> Poly(x**3 + 2*x - 1, x).all_terms()
          [((3,), 1), ((2,), 0), ((1,), 2), ((0,), -1)]
  
          """
          return [(m, f.rep.dom.to_sympy(c)) for m, c in f.rep.all_terms()]
  
      def termwise(f, func, *gens, **args):
          """
          Apply a function to all terms of ``f``.
  
          Examples
          ========
  
          >>> from sympy import Poly
          >>> from sympy.abc import x
  
          >>> def func(k, coeff):
          ...     k = k[0]
          ...     return coeff//10**(2-k)
  
          >>> Poly(x**2 + 20*x + 400).termwise(func)
          Poly(x**2 + 2*x + 4, x, domain='ZZ')
  
          """
          terms = {}
  
          for monom, coeff in f.terms():
              result = func(monom, coeff)
  
              if isinstance(result, tuple):
                  monom, coeff = result
              else:
                  coeff = result
  
              if coeff:
                  if monom not in terms:
                      terms[monom] = coeff
                  else:
                      raise PolynomialError(
                          "%s monomial was generated twice" % monom)
  
          return f.from_dict(terms, *(gens or f.gens), **args)
  
      def length(f):
          """
          Returns the number of non-zero terms in ``f``.
  
          Examples
          ========
  
          >>> from sympy import Poly
          >>> from sympy.abc import x
  
          >>> Poly(x**2 + 2*x - 1).length()
          3
  
          """
          return len(f.as_dict())
  
      def as_dict(f, native=False, zero=False):
          """
          Switch to a ``dict`` representation.
  
          Examples
          ========
  
          >>> from sympy import Poly
          >>> from sympy.abc import x, y
  
          >>> Poly(x**2 + 2*x*y**2 - y, x, y).as_dict()
          {(0, 1): -1, (1, 2): 2, (2, 0): 1}
  
          """
          if native:
              return f.rep.to_dict(zero=zero)
          else:
              return f.rep.to_sympy_dict(zero=zero)
  
      def as_list(f, native=False):
          """Switch to a ``list`` representation. """
          if native:
              return f.rep.to_list()
          else:
              return f.rep.to_sympy_list()
  
      def as_expr(f, *gens):
          """
          Convert a Poly instance to an Expr instance.
  
          Examples
          ========
  
          >>> from sympy import Poly
          >>> from sympy.abc import x, y
  
          >>> f = Poly(x**2 + 2*x*y**2 - y, x, y)
  
          >>> f.as_expr()
          x**2 + 2*x*y**2 - y
          >>> f.as_expr({x: 5})
          10*y**2 - y + 25
          >>> f.as_expr(5, 6)
          379
  
          """
          if not gens:
              return f.expr
  
          if len(gens) == 1 and isinstance(gens[0], dict):
              mapping = gens[0]
              gens = list(f.gens)
  
              for gen, value in mapping.items():
                  try:
                      index = gens.index(gen)
                  except ValueError:
                      raise GeneratorsError(
                          "%s doesn't have %s as generator" % (f, gen))
                  else:
                      gens[index] = value
  
          return basic_from_dict(f.rep.to_sympy_dict(), *gens)
  
      def as_poly(self, *gens, **args):
          """Converts ``self`` to a polynomial or returns ``None``.
  
          >>> from sympy import sin
          >>> from sympy.abc import x, y
  
          >>> print((x**2 + x*y).as_poly())
          Poly(x**2 + x*y, x, y, domain='ZZ')
  
          >>> print((x**2 + x*y).as_poly(x, y))
          Poly(x**2 + x*y, x, y, domain='ZZ')
  
          >>> print((x**2 + sin(y)).as_poly(x, y))
          None
  
          """
          try:
              poly = Poly(self, *gens, **args)
  
              if not poly.is_Poly:
                  return None
              else:
                  return poly
          except PolynomialError:
              return None
  
      def lift(f):
          """
          Convert algebraic coefficients to rationals.
  
          Examples
          ========
  
          >>> from sympy import Poly, I
          >>> from sympy.abc import x
  
          >>> Poly(x**2 + I*x + 1, x, extension=I).lift()
          Poly(x**4 + 3*x**2 + 1, x, domain='QQ')
  
          """
          if hasattr(f.rep, 'lift'):
              result = f.rep.lift()
          else:  # pragma: no cover
              raise OperationNotSupported(f, 'lift')
  
          return f.per(result)
  
      def deflate(f):
          """
          Reduce degree of ``f`` by mapping ``x_i**m`` to ``y_i``.
  
          Examples
          ========
  
          >>> from sympy import Poly
          >>> from sympy.abc import x, y
  
          >>> Poly(x**6*y**2 + x**3 + 1, x, y).deflate()
          ((3, 2), Poly(x**2*y + x + 1, x, y, domain='ZZ'))
  
          """
          if hasattr(f.rep, 'deflate'):
              J, result = f.rep.deflate()
          else:  # pragma: no cover
              raise OperationNotSupported(f, 'deflate')
  
          return J, f.per(result)
  
      def inject(f, front=False):
          """
          Inject ground domain generators into ``f``.
  
          Examples
          ========
  
          >>> from sympy import Poly
          >>> from sympy.abc import x, y
  
          >>> f = Poly(x**2*y + x*y**3 + x*y + 1, x)
  
          >>> f.inject()
          Poly(x**2*y + x*y**3 + x*y + 1, x, y, domain='ZZ')
          >>> f.inject(front=True)
          Poly(y**3*x + y*x**2 + y*x + 1, y, x, domain='ZZ')
  
          """
          dom = f.rep.dom
  
          if dom.is_Numerical:
              return f
          elif not dom.is_Poly:
              raise DomainError("can't inject generators over %s" % dom)
  
          if hasattr(f.rep, 'inject'):
              result = f.rep.inject(front=front)
          else:  # pragma: no cover
              raise OperationNotSupported(f, 'inject')
  
          if front:
              gens = dom.symbols + f.gens
          else:
              gens = f.gens + dom.symbols
  
          return f.new(result, *gens)
  
      def eject(f, *gens):
          """
          Eject selected generators into the ground domain.
  
          Examples
          ========
  
          >>> from sympy import Poly
          >>> from sympy.abc import x, y
  
          >>> f = Poly(x**2*y + x*y**3 + x*y + 1, x, y)
  
          >>> f.eject(x)
          Poly(x*y**3 + (x**2 + x)*y + 1, y, domain='ZZ[x]')
          >>> f.eject(y)
          Poly(y*x**2 + (y**3 + y)*x + 1, x, domain='ZZ[y]')
  
          """
          dom = f.rep.dom
  
          if not dom.is_Numerical:
              raise DomainError("can't eject generators over %s" % dom)
  
          k = len(gens)
  
          if f.gens[:k] == gens:
              _gens, front = f.gens[k:], True
          elif f.gens[-k:] == gens:
              _gens, front = f.gens[:-k], False
          else:
              raise NotImplementedError(
                  "can only eject front or back generators")
  
          dom = dom.inject(*gens)
  
          if hasattr(f.rep, 'eject'):
              result = f.rep.eject(dom, front=front)
          else:  # pragma: no cover
              raise OperationNotSupported(f, 'eject')
  
          return f.new(result, *_gens)
  
      def terms_gcd(f):
          """
          Remove GCD of terms from the polynomial ``f``.
  
          Examples
          ========
  
          >>> from sympy import Poly
          >>> from sympy.abc import x, y
  
          >>> Poly(x**6*y**2 + x**3*y, x, y).terms_gcd()
          ((3, 1), Poly(x**3*y + 1, x, y, domain='ZZ'))
  
          """
          if hasattr(f.rep, 'terms_gcd'):
              J, result = f.rep.terms_gcd()
          else:  # pragma: no cover
              raise OperationNotSupported(f, 'terms_gcd')
  
          return J, f.per(result)
  
      def add_ground(f, coeff):
          """
          Add an element of the ground domain to ``f``.
  
          Examples
          ========
  
          >>> from sympy import Poly
          >>> from sympy.abc import x
  
          >>> Poly(x + 1).add_ground(2)
          Poly(x + 3, x, domain='ZZ')
  
          """
          if hasattr(f.rep, 'add_ground'):
              result = f.rep.add_ground(coeff)
          else:  # pragma: no cover
              raise OperationNotSupported(f, 'add_ground')
  
          return f.per(result)
  
      def sub_ground(f, coeff):
          """
          Subtract an element of the ground domain from ``f``.
  
          Examples
          ========
  
          >>> from sympy import Poly
          >>> from sympy.abc import x
  
          >>> Poly(x + 1).sub_ground(2)
          Poly(x - 1, x, domain='ZZ')
  
          """
          if hasattr(f.rep, 'sub_ground'):
              result = f.rep.sub_ground(coeff)
          else:  # pragma: no cover
              raise OperationNotSupported(f, 'sub_ground')
  
          return f.per(result)
  
      def mul_ground(f, coeff):
          """
          Multiply ``f`` by a an element of the ground domain.
  
          Examples
          ========
  
          >>> from sympy import Poly
          >>> from sympy.abc import x
  
          >>> Poly(x + 1).mul_ground(2)
          Poly(2*x + 2, x, domain='ZZ')
  
          """
          if hasattr(f.rep, 'mul_ground'):
              result = f.rep.mul_ground(coeff)
          else:  # pragma: no cover
              raise OperationNotSupported(f, 'mul_ground')
  
          return f.per(result)
  
      def quo_ground(f, coeff):
          """
          Quotient of ``f`` by a an element of the ground domain.
  
          Examples
          ========
  
          >>> from sympy import Poly
          >>> from sympy.abc import x
  
          >>> Poly(2*x + 4).quo_ground(2)
          Poly(x + 2, x, domain='ZZ')
  
          >>> Poly(2*x + 3).quo_ground(2)
          Poly(x + 1, x, domain='ZZ')
  
          """
          if hasattr(f.rep, 'quo_ground'):
              result = f.rep.quo_ground(coeff)
          else:  # pragma: no cover
              raise OperationNotSupported(f, 'quo_ground')
  
          return f.per(result)
  
      def exquo_ground(f, coeff):
          """
          Exact quotient of ``f`` by a an element of the ground domain.
  
          Examples
          ========
  
          >>> from sympy import Poly
          >>> from sympy.abc import x
  
          >>> Poly(2*x + 4).exquo_ground(2)
          Poly(x + 2, x, domain='ZZ')
  
          >>> Poly(2*x + 3).exquo_ground(2)
          Traceback (most recent call last):
          ...
          ExactQuotientFailed: 2 does not divide 3 in ZZ
  
          """
          if hasattr(f.rep, 'exquo_ground'):
              result = f.rep.exquo_ground(coeff)
          else:  # pragma: no cover
              raise OperationNotSupported(f, 'exquo_ground')
  
          return f.per(result)
  
      def abs(f):
          """
          Make all coefficients in ``f`` positive.
  
          Examples
          ========
  
          >>> from sympy import Poly
          >>> from sympy.abc import x
  
          >>> Poly(x**2 - 1, x).abs()
          Poly(x**2 + 1, x, domain='ZZ')
  
          """
          if hasattr(f.rep, 'abs'):
              result = f.rep.abs()
          else:  # pragma: no cover
              raise OperationNotSupported(f, 'abs')
  
          return f.per(result)
  
      def neg(f):
          """
          Negate all coefficients in ``f``.
  
          Examples
          ========
  
          >>> from sympy import Poly
          >>> from sympy.abc import x
  
          >>> Poly(x**2 - 1, x).neg()
          Poly(-x**2 + 1, x, domain='ZZ')
  
          >>> -Poly(x**2 - 1, x)
          Poly(-x**2 + 1, x, domain='ZZ')
  
          """
          if hasattr(f.rep, 'neg'):
              result = f.rep.neg()
          else:  # pragma: no cover
              raise OperationNotSupported(f, 'neg')
  
          return f.per(result)
  
      def add(f, g):
          """
          Add two polynomials ``f`` and ``g``.
  
          Examples
          ========
  
          >>> from sympy import Poly
          >>> from sympy.abc import x
  
          >>> Poly(x**2 + 1, x).add(Poly(x - 2, x))
          Poly(x**2 + x - 1, x, domain='ZZ')
  
          >>> Poly(x**2 + 1, x) + Poly(x - 2, x)
          Poly(x**2 + x - 1, x, domain='ZZ')
  
          """
          g = sympify(g)
  
          if not g.is_Poly:
              return f.add_ground(g)
  
          _, per, F, G = f._unify(g)
  
          if hasattr(f.rep, 'add'):
              result = F.add(G)
          else:  # pragma: no cover
              raise OperationNotSupported(f, 'add')
  
          return per(result)
  
      def sub(f, g):
          """
          Subtract two polynomials ``f`` and ``g``.
  
          Examples
          ========
  
          >>> from sympy import Poly
          >>> from sympy.abc import x
  
          >>> Poly(x**2 + 1, x).sub(Poly(x - 2, x))
          Poly(x**2 - x + 3, x, domain='ZZ')
  
          >>> Poly(x**2 + 1, x) - Poly(x - 2, x)
          Poly(x**2 - x + 3, x, domain='ZZ')
  
          """
          g = sympify(g)
  
          if not g.is_Poly:
              return f.sub_ground(g)
  
          _, per, F, G = f._unify(g)
  
          if hasattr(f.rep, 'sub'):
              result = F.sub(G)
          else:  # pragma: no cover
              raise OperationNotSupported(f, 'sub')
  
          return per(result)
  
      def mul(f, g):
          """
          Multiply two polynomials ``f`` and ``g``.
  
          Examples
          ========
  
          >>> from sympy import Poly
          >>> from sympy.abc import x
  
          >>> Poly(x**2 + 1, x).mul(Poly(x - 2, x))
          Poly(x**3 - 2*x**2 + x - 2, x, domain='ZZ')
  
          >>> Poly(x**2 + 1, x)*Poly(x - 2, x)
          Poly(x**3 - 2*x**2 + x - 2, x, domain='ZZ')
  
          """
          g = sympify(g)
  
          if not g.is_Poly:
              return f.mul_ground(g)
  
          _, per, F, G = f._unify(g)
  
          if hasattr(f.rep, 'mul'):
              result = F.mul(G)
          else:  # pragma: no cover
              raise OperationNotSupported(f, 'mul')
  
          return per(result)
  
      def sqr(f):
          """
          Square a polynomial ``f``.
  
          Examples
          ========
  
          >>> from sympy import Poly
          >>> from sympy.abc import x
  
          >>> Poly(x - 2, x).sqr()
          Poly(x**2 - 4*x + 4, x, domain='ZZ')
  
          >>> Poly(x - 2, x)**2
          Poly(x**2 - 4*x + 4, x, domain='ZZ')
  
          """
          if hasattr(f.rep, 'sqr'):
              result = f.rep.sqr()
          else:  # pragma: no cover
              raise OperationNotSupported(f, 'sqr')
  
          return f.per(result)
  
      def pow(f, n):
          """
          Raise ``f`` to a non-negative power ``n``.
  
          Examples
          ========
  
          >>> from sympy import Poly
          >>> from sympy.abc import x
  
          >>> Poly(x - 2, x).pow(3)
          Poly(x**3 - 6*x**2 + 12*x - 8, x, domain='ZZ')
  
          >>> Poly(x - 2, x)**3
          Poly(x**3 - 6*x**2 + 12*x - 8, x, domain='ZZ')
  
          """
          n = int(n)
  
          if hasattr(f.rep, 'pow'):
              result = f.rep.pow(n)
          else:  # pragma: no cover
              raise OperationNotSupported(f, 'pow')
  
          return f.per(result)
  
      def pdiv(f, g):
          """
          Polynomial pseudo-division of ``f`` by ``g``.
  
          Examples
          ========
  
          >>> from sympy import Poly
          >>> from sympy.abc import x
  
          >>> Poly(x**2 + 1, x).pdiv(Poly(2*x - 4, x))
          (Poly(2*x + 4, x, domain='ZZ'), Poly(20, x, domain='ZZ'))
  
          """
          _, per, F, G = f._unify(g)
  
          if hasattr(f.rep, 'pdiv'):
              q, r = F.pdiv(G)
          else:  # pragma: no cover
              raise OperationNotSupported(f, 'pdiv')
  
          return per(q), per(r)
  
      def prem(f, g):
          """
          Polynomial pseudo-remainder of ``f`` by ``g``.
  
          Caveat: The function prem(f, g, x) can be safely used to compute
            in Z[x] _only_ subresultant polynomial remainder sequences (prs's).
  
            To safely compute Euclidean and Sturmian prs's in Z[x]
            employ anyone of the corresponding functions found in
            the module sympy.polys.subresultants_qq_zz. The functions
            in the module with suffix _pg compute prs's in Z[x] employing
            rem(f, g, x), whereas the functions with suffix _amv
            compute prs's in Z[x] employing rem_z(f, g, x).
  
            The function rem_z(f, g, x) differs from prem(f, g, x) in that
            to compute the remainder polynomials in Z[x] it premultiplies
            the divident times the absolute value of the leading coefficient
            of the divisor raised to the power degree(f, x) - degree(g, x) + 1.
  
  
          Examples
          ========
  
          >>> from sympy import Poly
          >>> from sympy.abc import x
  
          >>> Poly(x**2 + 1, x).prem(Poly(2*x - 4, x))
          Poly(20, x, domain='ZZ')
  
          """
          _, per, F, G = f._unify(g)
  
          if hasattr(f.rep, 'prem'):
              result = F.prem(G)
          else:  # pragma: no cover
              raise OperationNotSupported(f, 'prem')
  
          return per(result)
  
      def pquo(f, g):
          """
          Polynomial pseudo-quotient of ``f`` by ``g``.
  
          See the Caveat note in the function prem(f, g).
  
          Examples
          ========
  
          >>> from sympy import Poly
          >>> from sympy.abc import x
  
          >>> Poly(x**2 + 1, x).pquo(Poly(2*x - 4, x))
          Poly(2*x + 4, x, domain='ZZ')
  
          >>> Poly(x**2 - 1, x).pquo(Poly(2*x - 2, x))
          Poly(2*x + 2, x, domain='ZZ')
  
          """
          _, per, F, G = f._unify(g)
  
          if hasattr(f.rep, 'pquo'):
              result = F.pquo(G)
          else:  # pragma: no cover
              raise OperationNotSupported(f, 'pquo')
  
          return per(result)
  
      def pexquo(f, g):
          """
          Polynomial exact pseudo-quotient of ``f`` by ``g``.
  
          Examples
          ========
  
          >>> from sympy import Poly
          >>> from sympy.abc import x
  
          >>> Poly(x**2 - 1, x).pexquo(Poly(2*x - 2, x))
          Poly(2*x + 2, x, domain='ZZ')
  
          >>> Poly(x**2 + 1, x).pexquo(Poly(2*x - 4, x))
          Traceback (most recent call last):
          ...
          ExactQuotientFailed: 2*x - 4 does not divide x**2 + 1
  
          """
          _, per, F, G = f._unify(g)
  
          if hasattr(f.rep, 'pexquo'):
              try:
                  result = F.pexquo(G)
              except ExactQuotientFailed as exc:
                  raise exc.new(f.as_expr(), g.as_expr())
          else:  # pragma: no cover
              raise OperationNotSupported(f, 'pexquo')
  
          return per(result)
  
      def div(f, g, auto=True):
          """
          Polynomial division with remainder of ``f`` by ``g``.
  
          Examples
          ========
  
          >>> from sympy import Poly
          >>> from sympy.abc import x
  
          >>> Poly(x**2 + 1, x).div(Poly(2*x - 4, x))
          (Poly(1/2*x + 1, x, domain='QQ'), Poly(5, x, domain='QQ'))
  
          >>> Poly(x**2 + 1, x).div(Poly(2*x - 4, x), auto=False)
          (Poly(0, x, domain='ZZ'), Poly(x**2 + 1, x, domain='ZZ'))
  
          """
          dom, per, F, G = f._unify(g)
          retract = False
  
          if auto and dom.is_Ring and not dom.is_Field:
              F, G = F.to_field(), G.to_field()
              retract = True
  
          if hasattr(f.rep, 'div'):
              q, r = F.div(G)
          else:  # pragma: no cover
              raise OperationNotSupported(f, 'div')
  
          if retract:
              try:
                  Q, R = q.to_ring(), r.to_ring()
              except CoercionFailed:
                  pass
              else:
                  q, r = Q, R
  
          return per(q), per(r)
  
      def rem(f, g, auto=True):
          """
          Computes the polynomial remainder of ``f`` by ``g``.
  
          Examples
          ========
  
          >>> from sympy import Poly
          >>> from sympy.abc import x
  
          >>> Poly(x**2 + 1, x).rem(Poly(2*x - 4, x))
          Poly(5, x, domain='ZZ')
  
          >>> Poly(x**2 + 1, x).rem(Poly(2*x - 4, x), auto=False)
          Poly(x**2 + 1, x, domain='ZZ')
  
          """
          dom, per, F, G = f._unify(g)
          retract = False
  
          if auto and dom.is_Ring and not dom.is_Field:
              F, G = F.to_field(), G.to_field()
              retract = True
  
          if hasattr(f.rep, 'rem'):
              r = F.rem(G)
          else:  # pragma: no cover
              raise OperationNotSupported(f, 'rem')
  
          if retract:
              try:
                  r = r.to_ring()
              except CoercionFailed:
                  pass
  
          return per(r)
  
      def quo(f, g, auto=True):
          """
          Computes polynomial quotient of ``f`` by ``g``.
  
          Examples
          ========
  
          >>> from sympy import Poly
          >>> from sympy.abc import x
  
          >>> Poly(x**2 + 1, x).quo(Poly(2*x - 4, x))
          Poly(1/2*x + 1, x, domain='QQ')
  
          >>> Poly(x**2 - 1, x).quo(Poly(x - 1, x))
          Poly(x + 1, x, domain='ZZ')
  
          """
          dom, per, F, G = f._unify(g)
          retract = False
  
          if auto and dom.is_Ring and not dom.is_Field:
              F, G = F.to_field(), G.to_field()
              retract = True
  
          if hasattr(f.rep, 'quo'):
              q = F.quo(G)
          else:  # pragma: no cover
              raise OperationNotSupported(f, 'quo')
  
          if retract:
              try:
                  q = q.to_ring()
              except CoercionFailed:
                  pass
  
          return per(q)
  
      def exquo(f, g, auto=True):
          """
          Computes polynomial exact quotient of ``f`` by ``g``.
  
          Examples
          ========
  
          >>> from sympy import Poly
          >>> from sympy.abc import x
  
          >>> Poly(x**2 - 1, x).exquo(Poly(x - 1, x))
          Poly(x + 1, x, domain='ZZ')
  
          >>> Poly(x**2 + 1, x).exquo(Poly(2*x - 4, x))
          Traceback (most recent call last):
          ...
          ExactQuotientFailed: 2*x - 4 does not divide x**2 + 1
  
          """
          dom, per, F, G = f._unify(g)
          retract = False
  
          if auto and dom.is_Ring and not dom.is_Field:
              F, G = F.to_field(), G.to_field()
              retract = True
  
          if hasattr(f.rep, 'exquo'):
              try:
                  q = F.exquo(G)
              except ExactQuotientFailed as exc:
                  raise exc.new(f.as_expr(), g.as_expr())
          else:  # pragma: no cover
              raise OperationNotSupported(f, 'exquo')
  
          if retract:
              try:
                  q = q.to_ring()
              except CoercionFailed:
                  pass
  
          return per(q)
  
      def _gen_to_level(f, gen):
          """Returns level associated with the given generator. """
          if isinstance(gen, int):
              length = len(f.gens)
  
              if -length <= gen < length:
                  if gen < 0:
                      return length + gen
                  else:
                      return gen
              else:
                  raise PolynomialError("-%s <= gen < %s expected, got %s" %
                                        (length, length, gen))
          else:
              try:
                  return f.gens.index(sympify(gen))
              except ValueError:
                  raise PolynomialError(
                      "a valid generator expected, got %s" % gen)
  
      def degree(f, gen=0):
          """
          Returns degree of ``f`` in ``x_j``.
  
          The degree of 0 is negative infinity.
  
          Examples
          ========
  
          >>> from sympy import Poly
          >>> from sympy.abc import x, y
  
          >>> Poly(x**2 + y*x + 1, x, y).degree()
          2
          >>> Poly(x**2 + y*x + y, x, y).degree(y)
          1
          >>> Poly(0, x).degree()
          -oo
  
          """
          j = f._gen_to_level(gen)
  
          if hasattr(f.rep, 'degree'):
              return f.rep.degree(j)
          else:  # pragma: no cover
              raise OperationNotSupported(f, 'degree')
  
      def degree_list(f):
          """
          Returns a list of degrees of ``f``.
  
          Examples
          ========
  
          >>> from sympy import Poly
          >>> from sympy.abc import x, y
  
          >>> Poly(x**2 + y*x + 1, x, y).degree_list()
          (2, 1)
  
          """
          if hasattr(f.rep, 'degree_list'):
              return f.rep.degree_list()
          else:  # pragma: no cover
              raise OperationNotSupported(f, 'degree_list')
  
      def total_degree(f):
          """
          Returns the total degree of ``f``.
  
          Examples
          ========
  
          >>> from sympy import Poly
          >>> from sympy.abc import x, y
  
          >>> Poly(x**2 + y*x + 1, x, y).total_degree()
          2
          >>> Poly(x + y**5, x, y).total_degree()
          5
  
          """
          if hasattr(f.rep, 'total_degree'):
              return f.rep.total_degree()
          else:  # pragma: no cover
              raise OperationNotSupported(f, 'total_degree')
  
      def homogenize(f, s):
          """
          Returns the homogeneous polynomial of ``f``.
  
          A homogeneous polynomial is a polynomial whose all monomials with
          non-zero coefficients have the same total degree. If you only
          want to check if a polynomial is homogeneous, then use
          :func:`Poly.is_homogeneous`. If you want not only to check if a
          polynomial is homogeneous but also compute its homogeneous order,
          then use :func:`Poly.homogeneous_order`.
  
          Examples
          ========
  
          >>> from sympy import Poly
          >>> from sympy.abc import x, y, z
  
          >>> f = Poly(x**5 + 2*x**2*y**2 + 9*x*y**3)
          >>> f.homogenize(z)
          Poly(x**5 + 2*x**2*y**2*z + 9*x*y**3*z, x, y, z, domain='ZZ')
  
          """
          if not isinstance(s, Symbol):
              raise TypeError("``Symbol`` expected, got %s" % type(s))
          if s in f.gens:
              i = f.gens.index(s)
              gens = f.gens
          else:
              i = len(f.gens)
              gens = f.gens + (s,)
          if hasattr(f.rep, 'homogenize'):
              return f.per(f.rep.homogenize(i), gens=gens)
          raise OperationNotSupported(f, 'homogeneous_order')
  
      def homogeneous_order(f):
          """
          Returns the homogeneous order of ``f``.
  
          A homogeneous polynomial is a polynomial whose all monomials with
          non-zero coefficients have the same total degree. This degree is
          the homogeneous order of ``f``. If you only want to check if a
          polynomial is homogeneous, then use :func:`Poly.is_homogeneous`.
  
          Examples
          ========
  
          >>> from sympy import Poly
          >>> from sympy.abc import x, y
  
          >>> f = Poly(x**5 + 2*x**3*y**2 + 9*x*y**4)
          >>> f.homogeneous_order()
          5
  
          """
          if hasattr(f.rep, 'homogeneous_order'):
              return f.rep.homogeneous_order()
          else:  # pragma: no cover
              raise OperationNotSupported(f, 'homogeneous_order')
  
      def LC(f, order=None):
          """
          Returns the leading coefficient of ``f``.
  
          Examples
          ========
  
          >>> from sympy import Poly
          >>> from sympy.abc import x
  
          >>> Poly(4*x**3 + 2*x**2 + 3*x, x).LC()
          4
  
          """
          if order is not None:
              return f.coeffs(order)[0]
  
          if hasattr(f.rep, 'LC'):
              result = f.rep.LC()
          else:  # pragma: no cover
              raise OperationNotSupported(f, 'LC')
  
          return f.rep.dom.to_sympy(result)
  
      def TC(f):
          """
          Returns the trailing coefficient of ``f``.
  
          Examples
          ========
  
          >>> from sympy import Poly
          >>> from sympy.abc import x
  
          >>> Poly(x**3 + 2*x**2 + 3*x, x).TC()
          0
  
          """
          if hasattr(f.rep, 'TC'):
              result = f.rep.TC()
          else:  # pragma: no cover
              raise OperationNotSupported(f, 'TC')
  
          return f.rep.dom.to_sympy(result)
  
      def EC(f, order=None):
          """
          Returns the last non-zero coefficient of ``f``.
  
          Examples
          ========
  
          >>> from sympy import Poly
          >>> from sympy.abc import x
  
          >>> Poly(x**3 + 2*x**2 + 3*x, x).EC()
          3
  
          """
          if hasattr(f.rep, 'coeffs'):
              return f.coeffs(order)[-1]
          else:  # pragma: no cover
              raise OperationNotSupported(f, 'EC')
  
      def coeff_monomial(f, monom):
          """
          Returns the coefficient of ``monom`` in ``f`` if there, else None.
  
          Examples
          ========
  
          >>> from sympy import Poly, exp
          >>> from sympy.abc import x, y
  
          >>> p = Poly(24*x*y*exp(8) + 23*x, x, y)
  
          >>> p.coeff_monomial(x)
          23
          >>> p.coeff_monomial(y)
          0
          >>> p.coeff_monomial(x*y)
          24*exp(8)
  
          Note that ``Expr.coeff()`` behaves differently, collecting terms
          if possible; the Poly must be converted to an Expr to use that
          method, however:
  
          >>> p.as_expr().coeff(x)
          24*y*exp(8) + 23
          >>> p.as_expr().coeff(y)
          24*x*exp(8)
          >>> p.as_expr().coeff(x*y)
          24*exp(8)
  
          See Also
          ========
          nth: more efficient query using exponents of the monomial's generators
  
          """
          return f.nth(*Monomial(monom, f.gens).exponents)
  
      def nth(f, *N):
          """
          Returns the ``n``-th coefficient of ``f`` where ``N`` are the
          exponents of the generators in the term of interest.
  
          Examples
          ========
  
          >>> from sympy import Poly, sqrt
          >>> from sympy.abc import x, y
  
          >>> Poly(x**3 + 2*x**2 + 3*x, x).nth(2)
          2
          >>> Poly(x**3 + 2*x*y**2 + y**2, x, y).nth(1, 2)
          2
          >>> Poly(4*sqrt(x)*y)
          Poly(4*y*(sqrt(x)), y, sqrt(x), domain='ZZ')
          >>> _.nth(1, 1)
          4
  
          See Also
          ========
          coeff_monomial
  
          """
          if hasattr(f.rep, 'nth'):
              if len(N) != len(f.gens):
                  raise ValueError('exponent of each generator must be specified')
              result = f.rep.nth(*list(map(int, N)))
          else:  # pragma: no cover
              raise OperationNotSupported(f, 'nth')
  
          return f.rep.dom.to_sympy(result)
  
      def coeff(f, x, n=1, right=False):
          # the semantics of coeff_monomial and Expr.coeff are different;
          # if someone is working with a Poly, they should be aware of the
          # differences and chose the method best suited for the query.
          # Alternatively, a pure-polys method could be written here but
          # at this time the ``right`` keyword would be ignored because Poly
          # doesn't work with non-commutatives.
          raise NotImplementedError(
              'Either convert to Expr with `as_expr` method '
              'to use Expr\'s coeff method or else use the '
              '`coeff_monomial` method of Polys.')
  
      def LM(f, order=None):
          """
          Returns the leading monomial of ``f``.
  
          The Leading monomial signifies the monomial having
          the highest power of the principal generator in the
          expression f.
  
          Examples
          ========
  
          >>> from sympy import Poly
          >>> from sympy.abc import x, y
  
          >>> Poly(4*x**2 + 2*x*y**2 + x*y + 3*y, x, y).LM()
          x**2*y**0
  
          """
          return Monomial(f.monoms(order)[0], f.gens)
  
      def EM(f, order=None):
          """
          Returns the last non-zero monomial of ``f``.
  
          Examples
          ========
  
          >>> from sympy import Poly
          >>> from sympy.abc import x, y
  
          >>> Poly(4*x**2 + 2*x*y**2 + x*y + 3*y, x, y).EM()
          x**0*y**1
  
          """
          return Monomial(f.monoms(order)[-1], f.gens)
  
      def LT(f, order=None):
          """
          Returns the leading term of ``f``.
  
          The Leading term signifies the term having
          the highest power of the principal generator in the
          expression f along with its coefficient.
  
          Examples
          ========
  
          >>> from sympy import Poly
          >>> from sympy.abc import x, y
  
          >>> Poly(4*x**2 + 2*x*y**2 + x*y + 3*y, x, y).LT()
          (x**2*y**0, 4)
  
          """
          monom, coeff = f.terms(order)[0]
          return Monomial(monom, f.gens), coeff
  
      def ET(f, order=None):
          """
          Returns the last non-zero term of ``f``.
  
          Examples
          ========
  
          >>> from sympy import Poly
          >>> from sympy.abc import x, y
  
          >>> Poly(4*x**2 + 2*x*y**2 + x*y + 3*y, x, y).ET()
          (x**0*y**1, 3)
  
          """
          monom, coeff = f.terms(order)[-1]
          return Monomial(monom, f.gens), coeff
  
      def max_norm(f):
          """
          Returns maximum norm of ``f``.
  
          Examples
          ========
  
          >>> from sympy import Poly
          >>> from sympy.abc import x
  
          >>> Poly(-x**2 + 2*x - 3, x).max_norm()
          3
  
          """
          if hasattr(f.rep, 'max_norm'):
              result = f.rep.max_norm()
          else:  # pragma: no cover
              raise OperationNotSupported(f, 'max_norm')
  
          return f.rep.dom.to_sympy(result)
  
      def l1_norm(f):
          """
          Returns l1 norm of ``f``.
  
          Examples
          ========
  
          >>> from sympy import Poly
          >>> from sympy.abc import x
  
          >>> Poly(-x**2 + 2*x - 3, x).l1_norm()
          6
  
          """
          if hasattr(f.rep, 'l1_norm'):
              result = f.rep.l1_norm()
          else:  # pragma: no cover
              raise OperationNotSupported(f, 'l1_norm')
  
          return f.rep.dom.to_sympy(result)
  
      def clear_denoms(self, convert=False):
          """
          Clear denominators, but keep the ground domain.
  
          Examples
          ========
  
          >>> from sympy import Poly, S, QQ
          >>> from sympy.abc import x
  
          >>> f = Poly(x/2 + S(1)/3, x, domain=QQ)
  
          >>> f.clear_denoms()
          (6, Poly(3*x + 2, x, domain='QQ'))
          >>> f.clear_denoms(convert=True)
          (6, Poly(3*x + 2, x, domain='ZZ'))
  
          """
          f = self
  
          if not f.rep.dom.is_Field:
              return S.One, f
  
          dom = f.get_domain()
          if dom.has_assoc_Ring:
              dom = f.rep.dom.get_ring()
  
          if hasattr(f.rep, 'clear_denoms'):
              coeff, result = f.rep.clear_denoms()
          else:  # pragma: no cover
              raise OperationNotSupported(f, 'clear_denoms')
  
          coeff, f = dom.to_sympy(coeff), f.per(result)
  
          if not convert or not dom.has_assoc_Ring:
              return coeff, f
          else:
              return coeff, f.to_ring()
  
      def rat_clear_denoms(self, g):
          """
          Clear denominators in a rational function ``f/g``.
  
          Examples
          ========
  
          >>> from sympy import Poly
          >>> from sympy.abc import x, y
  
          >>> f = Poly(x**2/y + 1, x)
          >>> g = Poly(x**3 + y, x)
  
          >>> p, q = f.rat_clear_denoms(g)
  
          >>> p
          Poly(x**2 + y, x, domain='ZZ[y]')
          >>> q
          Poly(y*x**3 + y**2, x, domain='ZZ[y]')
  
          """
          f = self
  
          dom, per, f, g = f._unify(g)
  
          f = per(f)
          g = per(g)
  
          if not (dom.is_Field and dom.has_assoc_Ring):
              return f, g
  
          a, f = f.clear_denoms(convert=True)
          b, g = g.clear_denoms(convert=True)
  
          f = f.mul_ground(b)
          g = g.mul_ground(a)
  
          return f, g
  
      def integrate(self, *specs, **args):
          """
          Computes indefinite integral of ``f``.
  
          Examples
          ========
  
          >>> from sympy import Poly
          >>> from sympy.abc import x, y
  
          >>> Poly(x**2 + 2*x + 1, x).integrate()
          Poly(1/3*x**3 + x**2 + x, x, domain='QQ')
  
          >>> Poly(x*y**2 + x, x, y).integrate((0, 1), (1, 0))
          Poly(1/2*x**2*y**2 + 1/2*x**2, x, y, domain='QQ')
  
          """
          f = self
  
          if args.get('auto', True) and f.rep.dom.is_Ring:
              f = f.to_field()
  
          if hasattr(f.rep, 'integrate'):
              if not specs:
                  return f.per(f.rep.integrate(m=1))
  
              rep = f.rep
  
              for spec in specs:
                  if type(spec) is tuple:
                      gen, m = spec
                  else:
                      gen, m = spec, 1
  
                  rep = rep.integrate(int(m), f._gen_to_level(gen))
  
              return f.per(rep)
          else:  # pragma: no cover
              raise OperationNotSupported(f, 'integrate')
  
      def diff(f, *specs, **kwargs):
          """
          Computes partial derivative of ``f``.
  
          Examples
          ========
  
          >>> from sympy import Poly
          >>> from sympy.abc import x, y
  
          >>> Poly(x**2 + 2*x + 1, x).diff()
          Poly(2*x + 2, x, domain='ZZ')
  
          >>> Poly(x*y**2 + x, x, y).diff((0, 0), (1, 1))
          Poly(2*x*y, x, y, domain='ZZ')
  
          """
          if not kwargs.get('evaluate', True):
              return Derivative(f, *specs, **kwargs)
  
          if hasattr(f.rep, 'diff'):
              if not specs:
                  return f.per(f.rep.diff(m=1))
  
              rep = f.rep
  
              for spec in specs:
                  if type(spec) is tuple:
                      gen, m = spec
                  else:
                      gen, m = spec, 1
  
                  rep = rep.diff(int(m), f._gen_to_level(gen))
  
              return f.per(rep)
          else:  # pragma: no cover
              raise OperationNotSupported(f, 'diff')
  
      _eval_derivative = diff
  
      def eval(self, x, a=None, auto=True):
          """
          Evaluate ``f`` at ``a`` in the given variable.
  
          Examples
          ========
  
          >>> from sympy import Poly
          >>> from sympy.abc import x, y, z
  
          >>> Poly(x**2 + 2*x + 3, x).eval(2)
          11
  
          >>> Poly(2*x*y + 3*x + y + 2, x, y).eval(x, 2)
          Poly(5*y + 8, y, domain='ZZ')
  
          >>> f = Poly(2*x*y + 3*x + y + 2*z, x, y, z)
  
          >>> f.eval({x: 2})
          Poly(5*y + 2*z + 6, y, z, domain='ZZ')
          >>> f.eval({x: 2, y: 5})
          Poly(2*z + 31, z, domain='ZZ')
          >>> f.eval({x: 2, y: 5, z: 7})
          45
  
          >>> f.eval((2, 5))
          Poly(2*z + 31, z, domain='ZZ')
          >>> f(2, 5)
          Poly(2*z + 31, z, domain='ZZ')
  
          """
          f = self
  
          if a is None:
              if isinstance(x, dict):
                  mapping = x
  
                  for gen, value in mapping.items():
                      f = f.eval(gen, value)
  
                  return f
              elif isinstance(x, (tuple, list)):
                  values = x
  
                  if len(values) > len(f.gens):
                      raise ValueError("too many values provided")
  
                  for gen, value in zip(f.gens, values):
                      f = f.eval(gen, value)
  
                  return f
              else:
                  j, a = 0, x
          else:
              j = f._gen_to_level(x)
  
          if not hasattr(f.rep, 'eval'):  # pragma: no cover
              raise OperationNotSupported(f, 'eval')
  
          try:
              result = f.rep.eval(a, j)
          except CoercionFailed:
              if not auto:
                  raise DomainError("can't evaluate at %s in %s" % (a, f.rep.dom))
              else:
                  a_domain, [a] = construct_domain([a])
                  new_domain = f.get_domain().unify_with_symbols(a_domain, f.gens)
  
                  f = f.set_domain(new_domain)
                  a = new_domain.convert(a, a_domain)
  
                  result = f.rep.eval(a, j)
  
          return f.per(result, remove=j)
  
      def __call__(f, *values):
          """
          Evaluate ``f`` at the give values.
  
          Examples
          ========
  
          >>> from sympy import Poly
          >>> from sympy.abc import x, y, z
  
          >>> f = Poly(2*x*y + 3*x + y + 2*z, x, y, z)
  
          >>> f(2)
          Poly(5*y + 2*z + 6, y, z, domain='ZZ')
          >>> f(2, 5)
          Poly(2*z + 31, z, domain='ZZ')
          >>> f(2, 5, 7)
          45
  
          """
          return f.eval(values)
  
      def half_gcdex(f, g, auto=True):
          """
          Half extended Euclidean algorithm of ``f`` and ``g``.
  
          Returns ``(s, h)`` such that ``h = gcd(f, g)`` and ``s*f = h (mod g)``.
  
          Examples
          ========
  
          >>> from sympy import Poly
          >>> from sympy.abc import x
  
          >>> f = x**4 - 2*x**3 - 6*x**2 + 12*x + 15
          >>> g = x**3 + x**2 - 4*x - 4
  
          >>> Poly(f).half_gcdex(Poly(g))
          (Poly(-1/5*x + 3/5, x, domain='QQ'), Poly(x + 1, x, domain='QQ'))
  
          """
          dom, per, F, G = f._unify(g)
  
          if auto and dom.is_Ring:
              F, G = F.to_field(), G.to_field()
  
          if hasattr(f.rep, 'half_gcdex'):
              s, h = F.half_gcdex(G)
          else:  # pragma: no cover
              raise OperationNotSupported(f, 'half_gcdex')
  
          return per(s), per(h)
  
      def gcdex(f, g, auto=True):
          """
          Extended Euclidean algorithm of ``f`` and ``g``.
  
          Returns ``(s, t, h)`` such that ``h = gcd(f, g)`` and ``s*f + t*g = h``.
  
          Examples
          ========
  
          >>> from sympy import Poly
          >>> from sympy.abc import x
  
          >>> f = x**4 - 2*x**3 - 6*x**2 + 12*x + 15
          >>> g = x**3 + x**2 - 4*x - 4
  
          >>> Poly(f).gcdex(Poly(g))
          (Poly(-1/5*x + 3/5, x, domain='QQ'),
           Poly(1/5*x**2 - 6/5*x + 2, x, domain='QQ'),
           Poly(x + 1, x, domain='QQ'))
  
          """
          dom, per, F, G = f._unify(g)
  
          if auto and dom.is_Ring:
              F, G = F.to_field(), G.to_field()
  
          if hasattr(f.rep, 'gcdex'):
              s, t, h = F.gcdex(G)
          else:  # pragma: no cover
              raise OperationNotSupported(f, 'gcdex')
  
          return per(s), per(t), per(h)
  
      def invert(f, g, auto=True):
          """
          Invert ``f`` modulo ``g`` when possible.
  
          Examples
          ========
  
          >>> from sympy import Poly
          >>> from sympy.abc import x
  
          >>> Poly(x**2 - 1, x).invert(Poly(2*x - 1, x))
          Poly(-4/3, x, domain='QQ')
  
          >>> Poly(x**2 - 1, x).invert(Poly(x - 1, x))
          Traceback (most recent call last):
          ...
          NotInvertible: zero divisor
  
          """
          dom, per, F, G = f._unify(g)
  
          if auto and dom.is_Ring:
              F, G = F.to_field(), G.to_field()
  
          if hasattr(f.rep, 'invert'):
              result = F.invert(G)
          else:  # pragma: no cover
              raise OperationNotSupported(f, 'invert')
  
          return per(result)
  
      def revert(f, n):
          """
          Compute ``f**(-1)`` mod ``x**n``.
  
          Examples
          ========
  
          >>> from sympy import Poly
          >>> from sympy.abc import x
  
          >>> Poly(1, x).revert(2)
          Poly(1, x, domain='ZZ')
  
          >>> Poly(1 + x, x).revert(1)
          Poly(1, x, domain='ZZ')
  
          >>> Poly(x**2 - 1, x).revert(1)
          Traceback (most recent call last):
          ...
          NotReversible: only unity is reversible in a ring
  
          >>> Poly(1/x, x).revert(1)
          Traceback (most recent call last):
          ...
          PolynomialError: 1/x contains an element of the generators set
  
          """
          if hasattr(f.rep, 'revert'):
              result = f.rep.revert(int(n))
          else:  # pragma: no cover
              raise OperationNotSupported(f, 'revert')
  
          return f.per(result)
  
      def subresultants(f, g):
          """
          Computes the subresultant PRS of ``f`` and ``g``.
  
          Examples
          ========
  
          >>> from sympy import Poly
          >>> from sympy.abc import x
  
          >>> Poly(x**2 + 1, x).subresultants(Poly(x**2 - 1, x))
          [Poly(x**2 + 1, x, domain='ZZ'),
           Poly(x**2 - 1, x, domain='ZZ'),
           Poly(-2, x, domain='ZZ')]
  
          """
          _, per, F, G = f._unify(g)
  
          if hasattr(f.rep, 'subresultants'):
              result = F.subresultants(G)
          else:  # pragma: no cover
              raise OperationNotSupported(f, 'subresultants')
  
          return list(map(per, result))
  
      def resultant(f, g, includePRS=False):
          """
          Computes the resultant of ``f`` and ``g`` via PRS.
  
          If includePRS=True, it includes the subresultant PRS in the result.
          Because the PRS is used to calculate the resultant, this is more
          efficient than calling :func:`subresultants` separately.
  
          Examples
          ========
  
          >>> from sympy import Poly
          >>> from sympy.abc import x
  
          >>> f = Poly(x**2 + 1, x)
  
          >>> f.resultant(Poly(x**2 - 1, x))
          4
          >>> f.resultant(Poly(x**2 - 1, x), includePRS=True)
          (4, [Poly(x**2 + 1, x, domain='ZZ'), Poly(x**2 - 1, x, domain='ZZ'),
               Poly(-2, x, domain='ZZ')])
  
          """
          _, per, F, G = f._unify(g)
  
          if hasattr(f.rep, 'resultant'):
              if includePRS:
                  result, R = F.resultant(G, includePRS=includePRS)
              else:
                  result = F.resultant(G)
          else:  # pragma: no cover
              raise OperationNotSupported(f, 'resultant')
  
          if includePRS:
              return (per(result, remove=0), list(map(per, R)))
          return per(result, remove=0)
  
      def discriminant(f):
          """
          Computes the discriminant of ``f``.
  
          Examples
          ========
  
          >>> from sympy import Poly
          >>> from sympy.abc import x
  
          >>> Poly(x**2 + 2*x + 3, x).discriminant()
          -8
  
          """
          if hasattr(f.rep, 'discriminant'):
              result = f.rep.discriminant()
          else:  # pragma: no cover
              raise OperationNotSupported(f, 'discriminant')
  
          return f.per(result, remove=0)
  
      def dispersionset(f, g=None):
          r"""Compute the *dispersion set* of two polynomials.
  
          For two polynomials `f(x)` and `g(x)` with `\deg f > 0`
          and `\deg g > 0` the dispersion set `\operatorname{J}(f, g)` is defined as:
  
          .. math::
              \operatorname{J}(f, g)
              & := \{a \in \mathbb{N}_0 | \gcd(f(x), g(x+a)) \neq 1\} \\
              &  = \{a \in \mathbb{N}_0 | \deg \gcd(f(x), g(x+a)) \geq 1\}
  
          For a single polynomial one defines `\operatorname{J}(f) := \operatorname{J}(f, f)`.
  
          Examples
          ========
  
          >>> from sympy import poly
          >>> from sympy.polys.dispersion import dispersion, dispersionset
          >>> from sympy.abc import x
  
          Dispersion set and dispersion of a simple polynomial:
  
          >>> fp = poly((x - 3)*(x + 3), x)
          >>> sorted(dispersionset(fp))
          [0, 6]
          >>> dispersion(fp)
          6
  
          Note that the definition of the dispersion is not symmetric:
  
          >>> fp = poly(x**4 - 3*x**2 + 1, x)
          >>> gp = fp.shift(-3)
          >>> sorted(dispersionset(fp, gp))
          [2, 3, 4]
          >>> dispersion(fp, gp)
          4
          >>> sorted(dispersionset(gp, fp))
          []
          >>> dispersion(gp, fp)
          -oo
  
          Computing the dispersion also works over field extensions:
  
          >>> from sympy import sqrt
          >>> fp = poly(x**2 + sqrt(5)*x - 1, x, domain='QQ<sqrt(5)>')
          >>> gp = poly(x**2 + (2 + sqrt(5))*x + sqrt(5), x, domain='QQ<sqrt(5)>')
          >>> sorted(dispersionset(fp, gp))
          [2]
          >>> sorted(dispersionset(gp, fp))
          [1, 4]
  
          We can even perform the computations for polynomials
          having symbolic coefficients:
  
          >>> from sympy.abc import a
          >>> fp = poly(4*x**4 + (4*a + 8)*x**3 + (a**2 + 6*a + 4)*x**2 + (a**2 + 2*a)*x, x)
          >>> sorted(dispersionset(fp))
          [0, 1]
  
          See Also
          ========
  
          dispersion
  
          References
          ==========
  
          1. [ManWright94]_
          2. [Koepf98]_
          3. [Abramov71]_
          4. [Man93]_
          """
          from sympy.polys.dispersion import dispersionset
          return dispersionset(f, g)
  
      def dispersion(f, g=None):
          r"""Compute the *dispersion* of polynomials.
  
          For two polynomials `f(x)` and `g(x)` with `\deg f > 0`
          and `\deg g > 0` the dispersion `\operatorname{dis}(f, g)` is defined as:
  
          .. math::
              \operatorname{dis}(f, g)
              & := \max\{ J(f,g) \cup \{0\} \} \\
              &  = \max\{ \{a \in \mathbb{N} | \gcd(f(x), g(x+a)) \neq 1\} \cup \{0\} \}
  
          and for a single polynomial `\operatorname{dis}(f) := \operatorname{dis}(f, f)`.
  
          Examples
          ========
  
          >>> from sympy import poly
          >>> from sympy.polys.dispersion import dispersion, dispersionset
          >>> from sympy.abc import x
  
          Dispersion set and dispersion of a simple polynomial:
  
          >>> fp = poly((x - 3)*(x + 3), x)
          >>> sorted(dispersionset(fp))
          [0, 6]
          >>> dispersion(fp)
          6
  
          Note that the definition of the dispersion is not symmetric:
  
          >>> fp = poly(x**4 - 3*x**2 + 1, x)
          >>> gp = fp.shift(-3)
          >>> sorted(dispersionset(fp, gp))
          [2, 3, 4]
          >>> dispersion(fp, gp)
          4
          >>> sorted(dispersionset(gp, fp))
          []
          >>> dispersion(gp, fp)
          -oo
  
          Computing the dispersion also works over field extensions:
  
          >>> from sympy import sqrt
          >>> fp = poly(x**2 + sqrt(5)*x - 1, x, domain='QQ<sqrt(5)>')
          >>> gp = poly(x**2 + (2 + sqrt(5))*x + sqrt(5), x, domain='QQ<sqrt(5)>')
          >>> sorted(dispersionset(fp, gp))
          [2]
          >>> sorted(dispersionset(gp, fp))
          [1, 4]
  
          We can even perform the computations for polynomials
          having symbolic coefficients:
  
          >>> from sympy.abc import a
          >>> fp = poly(4*x**4 + (4*a + 8)*x**3 + (a**2 + 6*a + 4)*x**2 + (a**2 + 2*a)*x, x)
          >>> sorted(dispersionset(fp))
          [0, 1]
  
          See Also
          ========
  
          dispersionset
  
          References
          ==========
  
          1. [ManWright94]_
          2. [Koepf98]_
          3. [Abramov71]_
          4. [Man93]_
          """
          from sympy.polys.dispersion import dispersion
          return dispersion(f, g)
  
      def cofactors(f, g):
          """
          Returns the GCD of ``f`` and ``g`` and their cofactors.
  
          Returns polynomials ``(h, cff, cfg)`` such that ``h = gcd(f, g)``, and
          ``cff = quo(f, h)`` and ``cfg = quo(g, h)`` are, so called, cofactors
          of ``f`` and ``g``.
  
          Examples
          ========
  
          >>> from sympy import Poly
          >>> from sympy.abc import x
  
          >>> Poly(x**2 - 1, x).cofactors(Poly(x**2 - 3*x + 2, x))
          (Poly(x - 1, x, domain='ZZ'),
           Poly(x + 1, x, domain='ZZ'),
           Poly(x - 2, x, domain='ZZ'))
  
          """
          _, per, F, G = f._unify(g)
  
          if hasattr(f.rep, 'cofactors'):
              h, cff, cfg = F.cofactors(G)
          else:  # pragma: no cover
              raise OperationNotSupported(f, 'cofactors')
  
          return per(h), per(cff), per(cfg)
  
      def gcd(f, g):
          """
          Returns the polynomial GCD of ``f`` and ``g``.
  
          Examples
          ========
  
          >>> from sympy import Poly
          >>> from sympy.abc import x
  
          >>> Poly(x**2 - 1, x).gcd(Poly(x**2 - 3*x + 2, x))
          Poly(x - 1, x, domain='ZZ')
  
          """
          _, per, F, G = f._unify(g)
  
          if hasattr(f.rep, 'gcd'):
              result = F.gcd(G)
          else:  # pragma: no cover
              raise OperationNotSupported(f, 'gcd')
  
          return per(result)
  
      def lcm(f, g):
          """
          Returns polynomial LCM of ``f`` and ``g``.
  
          Examples
          ========
  
          >>> from sympy import Poly
          >>> from sympy.abc import x
  
          >>> Poly(x**2 - 1, x).lcm(Poly(x**2 - 3*x + 2, x))
          Poly(x**3 - 2*x**2 - x + 2, x, domain='ZZ')
  
          """
          _, per, F, G = f._unify(g)
  
          if hasattr(f.rep, 'lcm'):
              result = F.lcm(G)
          else:  # pragma: no cover
              raise OperationNotSupported(f, 'lcm')
  
          return per(result)
  
      def trunc(f, p):
          """
          Reduce ``f`` modulo a constant ``p``.
  
          Examples
          ========
  
          >>> from sympy import Poly
          >>> from sympy.abc import x
  
          >>> Poly(2*x**3 + 3*x**2 + 5*x + 7, x).trunc(3)
          Poly(-x**3 - x + 1, x, domain='ZZ')
  
          """
          p = f.rep.dom.convert(p)
  
          if hasattr(f.rep, 'trunc'):
              result = f.rep.trunc(p)
          else:  # pragma: no cover
              raise OperationNotSupported(f, 'trunc')
  
          return f.per(result)
  
      def monic(self, auto=True):
          """
          Divides all coefficients by ``LC(f)``.
  
          Examples
          ========
  
          >>> from sympy import Poly, ZZ
          >>> from sympy.abc import x
  
          >>> Poly(3*x**2 + 6*x + 9, x, domain=ZZ).monic()
          Poly(x**2 + 2*x + 3, x, domain='QQ')
  
          >>> Poly(3*x**2 + 4*x + 2, x, domain=ZZ).monic()
          Poly(x**2 + 4/3*x + 2/3, x, domain='QQ')
  
          """
          f = self
  
          if auto and f.rep.dom.is_Ring:
              f = f.to_field()
  
          if hasattr(f.rep, 'monic'):
              result = f.rep.monic()
          else:  # pragma: no cover
              raise OperationNotSupported(f, 'monic')
  
          return f.per(result)
  
      def content(f):
          """
          Returns the GCD of polynomial coefficients.
  
          Examples
          ========
  
          >>> from sympy import Poly
          >>> from sympy.abc import x
  
          >>> Poly(6*x**2 + 8*x + 12, x).content()
          2
  
          """
          if hasattr(f.rep, 'content'):
              result = f.rep.content()
          else:  # pragma: no cover
              raise OperationNotSupported(f, 'content')
  
          return f.rep.dom.to_sympy(result)
  
      def primitive(f):
          """
          Returns the content and a primitive form of ``f``.
  
          Examples
          ========
  
          >>> from sympy import Poly
          >>> from sympy.abc import x
  
          >>> Poly(2*x**2 + 8*x + 12, x).primitive()
          (2, Poly(x**2 + 4*x + 6, x, domain='ZZ'))
  
          """
          if hasattr(f.rep, 'primitive'):
              cont, result = f.rep.primitive()
          else:  # pragma: no cover
              raise OperationNotSupported(f, 'primitive')
  
          return f.rep.dom.to_sympy(cont), f.per(result)
  
      def compose(f, g):
          """
          Computes the functional composition of ``f`` and ``g``.
  
          Examples
          ========
  
          >>> from sympy import Poly
          >>> from sympy.abc import x
  
          >>> Poly(x**2 + x, x).compose(Poly(x - 1, x))
          Poly(x**2 - x, x, domain='ZZ')
  
          """
          _, per, F, G = f._unify(g)
  
          if hasattr(f.rep, 'compose'):
              result = F.compose(G)
          else:  # pragma: no cover
              raise OperationNotSupported(f, 'compose')
  
          return per(result)
  
      def decompose(f):
          """
          Computes a functional decomposition of ``f``.
  
          Examples
          ========
  
          >>> from sympy import Poly
          >>> from sympy.abc import x
  
          >>> Poly(x**4 + 2*x**3 - x - 1, x, domain='ZZ').decompose()
          [Poly(x**2 - x - 1, x, domain='ZZ'), Poly(x**2 + x, x, domain='ZZ')]
  
          """
          if hasattr(f.rep, 'decompose'):
              result = f.rep.decompose()
          else:  # pragma: no cover
              raise OperationNotSupported(f, 'decompose')
  
          return list(map(f.per, result))
  
      def shift(f, a):
          """
          Efficiently compute Taylor shift ``f(x + a)``.
  
          Examples
          ========
  
          >>> from sympy import Poly
          >>> from sympy.abc import x
  
          >>> Poly(x**2 - 2*x + 1, x).shift(2)
          Poly(x**2 + 2*x + 1, x, domain='ZZ')
  
          """
          if hasattr(f.rep, 'shift'):
              result = f.rep.shift(a)
          else:  # pragma: no cover
              raise OperationNotSupported(f, 'shift')
  
          return f.per(result)
  
      def transform(f, p, q):
          """
          Efficiently evaluate the functional transformation ``q**n * f(p/q)``.
  
  
          Examples
          ========
  
          >>> from sympy import Poly
          >>> from sympy.abc import x
  
          >>> Poly(x**2 - 2*x + 1, x).transform(Poly(x + 1, x), Poly(x - 1, x))
          Poly(4, x, domain='ZZ')
  
          """
          P, Q = p.unify(q)
          F, P = f.unify(P)
          F, Q = F.unify(Q)
  
          if hasattr(F.rep, 'transform'):
              result = F.rep.transform(P.rep, Q.rep)
          else:  # pragma: no cover
              raise OperationNotSupported(F, 'transform')
  
          return F.per(result)
  
      def sturm(self, auto=True):
          """
          Computes the Sturm sequence of ``f``.
  
          Examples
          ========
  
          >>> from sympy import Poly
          >>> from sympy.abc import x
  
          >>> Poly(x**3 - 2*x**2 + x - 3, x).sturm()
          [Poly(x**3 - 2*x**2 + x - 3, x, domain='QQ'),
           Poly(3*x**2 - 4*x + 1, x, domain='QQ'),
           Poly(2/9*x + 25/9, x, domain='QQ'),
           Poly(-2079/4, x, domain='QQ')]
  
          """
          f = self
  
          if auto and f.rep.dom.is_Ring:
              f = f.to_field()
  
          if hasattr(f.rep, 'sturm'):
              result = f.rep.sturm()
          else:  # pragma: no cover
              raise OperationNotSupported(f, 'sturm')
  
          return list(map(f.per, result))
  
      def gff_list(f):
          """
          Computes greatest factorial factorization of ``f``.
  
          Examples
          ========
  
          >>> from sympy import Poly
          >>> from sympy.abc import x
  
          >>> f = x**5 + 2*x**4 - x**3 - 2*x**2
  
          >>> Poly(f).gff_list()
          [(Poly(x, x, domain='ZZ'), 1), (Poly(x + 2, x, domain='ZZ'), 4)]
  
          """
          if hasattr(f.rep, 'gff_list'):
              result = f.rep.gff_list()
          else:  # pragma: no cover
              raise OperationNotSupported(f, 'gff_list')
  
          return [(f.per(g), k) for g, k in result]
  
      def norm(f):
          """
          Computes the product, ``Norm(f)``, of the conjugates of
          a polynomial ``f`` defined over a number field ``K``.
  
          Examples
          ========
  
          >>> from sympy import Poly, sqrt
          >>> from sympy.abc import x
  
          >>> a, b = sqrt(2), sqrt(3)
  
          A polynomial over a quadratic extension.
          Two conjugates x - a and x + a.
  
          >>> f = Poly(x - a, x, extension=a)
          >>> f.norm()
          Poly(x**2 - 2, x, domain='QQ')
  
          A polynomial over a quartic extension.
          Four conjugates x - a, x - a, x + a and x + a.
  
          >>> f = Poly(x - a, x, extension=(a, b))
          >>> f.norm()
          Poly(x**4 - 4*x**2 + 4, x, domain='QQ')
  
          """
          if hasattr(f.rep, 'norm'):
              r = f.rep.norm()
          else:  # pragma: no cover
              raise OperationNotSupported(f, 'norm')
  
          return f.per(r)
  
      def sqf_norm(f):
          """
          Computes square-free norm of ``f``.
  
          Returns ``s``, ``f``, ``r``, such that ``g(x) = f(x-sa)`` and
          ``r(x) = Norm(g(x))`` is a square-free polynomial over ``K``,
          where ``a`` is the algebraic extension of the ground domain.
  
          Examples
          ========
  
          >>> from sympy import Poly, sqrt
          >>> from sympy.abc import x
  
          >>> s, f, r = Poly(x**2 + 1, x, extension=[sqrt(3)]).sqf_norm()
  
          >>> s
          1
          >>> f
          Poly(x**2 - 2*sqrt(3)*x + 4, x, domain='QQ<sqrt(3)>')
          >>> r
          Poly(x**4 - 4*x**2 + 16, x, domain='QQ')
  
          """
          if hasattr(f.rep, 'sqf_norm'):
              s, g, r = f.rep.sqf_norm()
          else:  # pragma: no cover
              raise OperationNotSupported(f, 'sqf_norm')
  
          return s, f.per(g), f.per(r)
  
      def sqf_part(f):
          """
          Computes square-free part of ``f``.
  
          Examples
          ========
  
          >>> from sympy import Poly
          >>> from sympy.abc import x
  
          >>> Poly(x**3 - 3*x - 2, x).sqf_part()
          Poly(x**2 - x - 2, x, domain='ZZ')
  
          """
          if hasattr(f.rep, 'sqf_part'):
              result = f.rep.sqf_part()
          else:  # pragma: no cover
              raise OperationNotSupported(f, 'sqf_part')
  
          return f.per(result)
  
      def sqf_list(f, all=False):
          """
          Returns a list of square-free factors of ``f``.
  
          Examples
          ========
  
          >>> from sympy import Poly
          >>> from sympy.abc import x
  
          >>> f = 2*x**5 + 16*x**4 + 50*x**3 + 76*x**2 + 56*x + 16
  
          >>> Poly(f).sqf_list()
          (2, [(Poly(x + 1, x, domain='ZZ'), 2),
               (Poly(x + 2, x, domain='ZZ'), 3)])
  
          >>> Poly(f).sqf_list(all=True)
          (2, [(Poly(1, x, domain='ZZ'), 1),
               (Poly(x + 1, x, domain='ZZ'), 2),
               (Poly(x + 2, x, domain='ZZ'), 3)])
  
          """
          if hasattr(f.rep, 'sqf_list'):
              coeff, factors = f.rep.sqf_list(all)
          else:  # pragma: no cover
              raise OperationNotSupported(f, 'sqf_list')
  
          return f.rep.dom.to_sympy(coeff), [(f.per(g), k) for g, k in factors]
  
      def sqf_list_include(f, all=False):
          """
          Returns a list of square-free factors of ``f``.
  
          Examples
          ========
  
          >>> from sympy import Poly, expand
          >>> from sympy.abc import x
  
          >>> f = expand(2*(x + 1)**3*x**4)
          >>> f
          2*x**7 + 6*x**6 + 6*x**5 + 2*x**4
  
          >>> Poly(f).sqf_list_include()
          [(Poly(2, x, domain='ZZ'), 1),
           (Poly(x + 1, x, domain='ZZ'), 3),
           (Poly(x, x, domain='ZZ'), 4)]
  
          >>> Poly(f).sqf_list_include(all=True)
          [(Poly(2, x, domain='ZZ'), 1),
           (Poly(1, x, domain='ZZ'), 2),
           (Poly(x + 1, x, domain='ZZ'), 3),
           (Poly(x, x, domain='ZZ'), 4)]
  
          """
          if hasattr(f.rep, 'sqf_list_include'):
              factors = f.rep.sqf_list_include(all)
          else:  # pragma: no cover
              raise OperationNotSupported(f, 'sqf_list_include')
  
          return [(f.per(g), k) for g, k in factors]
  
      def factor_list(f):
          """
          Returns a list of irreducible factors of ``f``.
  
          Examples
          ========
  
          >>> from sympy import Poly
          >>> from sympy.abc import x, y
  
          >>> f = 2*x**5 + 2*x**4*y + 4*x**3 + 4*x**2*y + 2*x + 2*y
  
          >>> Poly(f).factor_list()
          (2, [(Poly(x + y, x, y, domain='ZZ'), 1),
               (Poly(x**2 + 1, x, y, domain='ZZ'), 2)])
  
          """
          if hasattr(f.rep, 'factor_list'):
              try:
                  coeff, factors = f.rep.factor_list()
              except DomainError:
                  return S.One, [(f, 1)]
          else:  # pragma: no cover
              raise OperationNotSupported(f, 'factor_list')
  
          return f.rep.dom.to_sympy(coeff), [(f.per(g), k) for g, k in factors]
  
      def factor_list_include(f):
          """
          Returns a list of irreducible factors of ``f``.
  
          Examples
          ========
  
          >>> from sympy import Poly
          >>> from sympy.abc import x, y
  
          >>> f = 2*x**5 + 2*x**4*y + 4*x**3 + 4*x**2*y + 2*x + 2*y
  
          >>> Poly(f).factor_list_include()
          [(Poly(2*x + 2*y, x, y, domain='ZZ'), 1),
           (Poly(x**2 + 1, x, y, domain='ZZ'), 2)]
  
          """
          if hasattr(f.rep, 'factor_list_include'):
              try:
                  factors = f.rep.factor_list_include()
              except DomainError:
                  return [(f, 1)]
          else:  # pragma: no cover
              raise OperationNotSupported(f, 'factor_list_include')
  
          return [(f.per(g), k) for g, k in factors]
  
      def intervals(f, all=False, eps=None, inf=None, sup=None, fast=False, sqf=False):
          """
          Compute isolating intervals for roots of ``f``.
  
          For real roots the Vincent-Akritas-Strzebonski (VAS) continued fractions method is used.
  
          References
          ==========
          .. [#] Alkiviadis G. Akritas and Adam W. Strzebonski: A Comparative Study of Two Real Root
              Isolation Methods . Nonlinear Analysis: Modelling and Control, Vol. 10, No. 4, 297-304, 2005.
          .. [#] Alkiviadis G. Akritas, Adam W. Strzebonski and Panagiotis S. Vigklas: Improving the
              Performance of the Continued Fractions Method Using new Bounds of Positive Roots. Nonlinear
              Analysis: Modelling and Control, Vol. 13, No. 3, 265-279, 2008.
  
          Examples
          ========
  
          >>> from sympy import Poly
          >>> from sympy.abc import x
  
          >>> Poly(x**2 - 3, x).intervals()
          [((-2, -1), 1), ((1, 2), 1)]
          >>> Poly(x**2 - 3, x).intervals(eps=1e-2)
          [((-26/15, -19/11), 1), ((19/11, 26/15), 1)]
  
          """
          if eps is not None:
              eps = QQ.convert(eps)
  
              if eps <= 0:
                  raise ValueError("'eps' must be a positive rational")
  
          if inf is not None:
              inf = QQ.convert(inf)
          if sup is not None:
              sup = QQ.convert(sup)
  
          if hasattr(f.rep, 'intervals'):
              result = f.rep.intervals(
                  all=all, eps=eps, inf=inf, sup=sup, fast=fast, sqf=sqf)
          else:  # pragma: no cover
              raise OperationNotSupported(f, 'intervals')
  
          if sqf:
              def _real(interval):
                  s, t = interval
                  return (QQ.to_sympy(s), QQ.to_sympy(t))
  
              if not all:
                  return list(map(_real, result))
  
              def _complex(rectangle):
                  (u, v), (s, t) = rectangle
                  return (QQ.to_sympy(u) + I*QQ.to_sympy(v),
                          QQ.to_sympy(s) + I*QQ.to_sympy(t))
  
              real_part, complex_part = result
  
              return list(map(_real, real_part)), list(map(_complex, complex_part))
          else:
              def _real(interval):
                  (s, t), k = interval
                  return ((QQ.to_sympy(s), QQ.to_sympy(t)), k)
  
              if not all:
                  return list(map(_real, result))
  
              def _complex(rectangle):
                  ((u, v), (s, t)), k = rectangle
                  return ((QQ.to_sympy(u) + I*QQ.to_sympy(v),
                           QQ.to_sympy(s) + I*QQ.to_sympy(t)), k)
  
              real_part, complex_part = result
  
              return list(map(_real, real_part)), list(map(_complex, complex_part))
  
      def refine_root(f, s, t, eps=None, steps=None, fast=False, check_sqf=False):
          """
          Refine an isolating interval of a root to the given precision.
  
          Examples
          ========
  
          >>> from sympy import Poly
          >>> from sympy.abc import x
  
          >>> Poly(x**2 - 3, x).refine_root(1, 2, eps=1e-2)
          (19/11, 26/15)
  
          """
          if check_sqf and not f.is_sqf:
              raise PolynomialError("only square-free polynomials supported")
  
          s, t = QQ.convert(s), QQ.convert(t)
  
          if eps is not None:
              eps = QQ.convert(eps)
  
              if eps <= 0:
                  raise ValueError("'eps' must be a positive rational")
  
          if steps is not None:
              steps = int(steps)
          elif eps is None:
              steps = 1
  
          if hasattr(f.rep, 'refine_root'):
              S, T = f.rep.refine_root(s, t, eps=eps, steps=steps, fast=fast)
          else:  # pragma: no cover
              raise OperationNotSupported(f, 'refine_root')
  
          return QQ.to_sympy(S), QQ.to_sympy(T)
  
      def count_roots(f, inf=None, sup=None):
          """
          Return the number of roots of ``f`` in ``[inf, sup]`` interval.
  
          Examples
          ========
  
          >>> from sympy import Poly, I
          >>> from sympy.abc import x
  
          >>> Poly(x**4 - 4, x).count_roots(-3, 3)
          2
          >>> Poly(x**4 - 4, x).count_roots(0, 1 + 3*I)
          1
  
          """
          inf_real, sup_real = True, True
  
          if inf is not None:
              inf = sympify(inf)
  
              if inf is S.NegativeInfinity:
                  inf = None
              else:
                  re, im = inf.as_real_imag()
  
                  if not im:
                      inf = QQ.convert(inf)
                  else:
                      inf, inf_real = list(map(QQ.convert, (re, im))), False
  
          if sup is not None:
              sup = sympify(sup)
  
              if sup is S.Infinity:
                  sup = None
              else:
                  re, im = sup.as_real_imag()
  
                  if not im:
                      sup = QQ.convert(sup)
                  else:
                      sup, sup_real = list(map(QQ.convert, (re, im))), False
  
          if inf_real and sup_real:
              if hasattr(f.rep, 'count_real_roots'):
                  count = f.rep.count_real_roots(inf=inf, sup=sup)
              else:  # pragma: no cover
                  raise OperationNotSupported(f, 'count_real_roots')
          else:
              if inf_real and inf is not None:
                  inf = (inf, QQ.zero)
  
              if sup_real and sup is not None:
                  sup = (sup, QQ.zero)
  
              if hasattr(f.rep, 'count_complex_roots'):
                  count = f.rep.count_complex_roots(inf=inf, sup=sup)
              else:  # pragma: no cover
                  raise OperationNotSupported(f, 'count_complex_roots')
  
          return Integer(count)
  
      def root(f, index, radicals=True):
          """
          Get an indexed root of a polynomial.
  
          Examples
          ========
  
          >>> from sympy import Poly
          >>> from sympy.abc import x
  
          >>> f = Poly(2*x**3 - 7*x**2 + 4*x + 4)
  
          >>> f.root(0)
          -1/2
          >>> f.root(1)
          2
          >>> f.root(2)
          2
          >>> f.root(3)
          Traceback (most recent call last):
          ...
          IndexError: root index out of [-3, 2] range, got 3
  
          >>> Poly(x**5 + x + 1).root(0)
          CRootOf(x**3 - x**2 + 1, 0)
  
          """
          return sympy.polys.rootoftools.rootof(f, index, radicals=radicals)
  
      def real_roots(f, multiple=True, radicals=True):
          """
          Return a list of real roots with multiplicities.
  
          Examples
          ========
  
          >>> from sympy import Poly
          >>> from sympy.abc import x
  
          >>> Poly(2*x**3 - 7*x**2 + 4*x + 4).real_roots()
          [-1/2, 2, 2]
          >>> Poly(x**3 + x + 1).real_roots()
          [CRootOf(x**3 + x + 1, 0)]
  
          """
          reals = sympy.polys.rootoftools.CRootOf.real_roots(f, radicals=radicals)
  
          if multiple:
              return reals
          else:
              return group(reals, multiple=False)
  
      def all_roots(f, multiple=True, radicals=True):
          """
          Return a list of real and complex roots with multiplicities.
  
          Examples
          ========
  
          >>> from sympy import Poly
          >>> from sympy.abc import x
  
          >>> Poly(2*x**3 - 7*x**2 + 4*x + 4).all_roots()
          [-1/2, 2, 2]
          >>> Poly(x**3 + x + 1).all_roots()
          [CRootOf(x**3 + x + 1, 0),
           CRootOf(x**3 + x + 1, 1),
           CRootOf(x**3 + x + 1, 2)]
  
          """
          roots = sympy.polys.rootoftools.CRootOf.all_roots(f, radicals=radicals)
  
          if multiple:
              return roots
          else:
              return group(roots, multiple=False)
  
      def nroots(f, n=15, maxsteps=50, cleanup=True):
          """
          Compute numerical approximations of roots of ``f``.
  
          Parameters
          ==========
  
          n ... the number of digits to calculate
          maxsteps ... the maximum number of iterations to do
  
          If the accuracy `n` cannot be reached in `maxsteps`, it will raise an
          exception. You need to rerun with higher maxsteps.
  
          Examples
          ========
  
          >>> from sympy import Poly
          >>> from sympy.abc import x
  
          >>> Poly(x**2 - 3).nroots(n=15)
          [-1.73205080756888, 1.73205080756888]
          >>> Poly(x**2 - 3).nroots(n=30)
          [-1.73205080756887729352744634151, 1.73205080756887729352744634151]
  
          """
          from sympy.functions.elementary.complexes import sign
          if f.is_multivariate:
              raise MultivariatePolynomialError(
                  "can't compute numerical roots of %s" % f)
  
          if f.degree() <= 0:
              return []
  
          # For integer and rational coefficients, convert them to integers only
          # (for accuracy). Otherwise just try to convert the coefficients to
          # mpmath.mpc and raise an exception if the conversion fails.
          if f.rep.dom is ZZ:
              coeffs = [int(coeff) for coeff in f.all_coeffs()]
          elif f.rep.dom is QQ:
              denoms = [coeff.q for coeff in f.all_coeffs()]
              from sympy.core.numbers import ilcm
              fac = ilcm(*denoms)
              coeffs = [int(coeff*fac) for coeff in f.all_coeffs()]
          else:
              coeffs = [coeff.evalf(n=n).as_real_imag()
                      for coeff in f.all_coeffs()]
              try:
                  coeffs = [mpmath.mpc(*coeff) for coeff in coeffs]
              except TypeError:
                  raise DomainError("Numerical domain expected, got %s" % \
                          f.rep.dom)
  
          dps = mpmath.mp.dps
          mpmath.mp.dps = n
  
          try:
              # We need to add extra precision to guard against losing accuracy.
              # 10 times the degree of the polynomial seems to work well.
              roots = mpmath.polyroots(coeffs, maxsteps=maxsteps,
                      cleanup=cleanup, error=False, extraprec=f.degree()*10)
  
              # Mpmath puts real roots first, then complex ones (as does all_roots)
              # so we make sure this convention holds here, too.
              roots = list(map(sympify,
                  sorted(roots, key=lambda r: (1 if r.imag else 0, r.real, abs(r.imag), sign(r.imag)))))
          except NoConvergence:
              raise NoConvergence(
                  'convergence to root failed; try n < %s or maxsteps > %s' % (
                  n, maxsteps))
          finally:
              mpmath.mp.dps = dps
  
          return roots
  
      def ground_roots(f):
          """
          Compute roots of ``f`` by factorization in the ground domain.
  
          Examples
          ========
  
          >>> from sympy import Poly
          >>> from sympy.abc import x
  
          >>> Poly(x**6 - 4*x**4 + 4*x**3 - x**2).ground_roots()
          {0: 2, 1: 2}
  
          """
          if f.is_multivariate:
              raise MultivariatePolynomialError(
                  "can't compute ground roots of %s" % f)
  
          roots = {}
  
          for factor, k in f.factor_list()[1]:
              if factor.is_linear:
                  a, b = factor.all_coeffs()
                  roots[-b/a] = k
  
          return roots
  
      def nth_power_roots_poly(f, n):
          """
          Construct a polynomial with n-th powers of roots of ``f``.
  
          Examples
          ========
  
          >>> from sympy import Poly
          >>> from sympy.abc import x
  
          >>> f = Poly(x**4 - x**2 + 1)
  
          >>> f.nth_power_roots_poly(2)
          Poly(x**4 - 2*x**3 + 3*x**2 - 2*x + 1, x, domain='ZZ')
          >>> f.nth_power_roots_poly(3)
          Poly(x**4 + 2*x**2 + 1, x, domain='ZZ')
          >>> f.nth_power_roots_poly(4)
          Poly(x**4 + 2*x**3 + 3*x**2 + 2*x + 1, x, domain='ZZ')
          >>> f.nth_power_roots_poly(12)
          Poly(x**4 - 4*x**3 + 6*x**2 - 4*x + 1, x, domain='ZZ')
  
          """
          if f.is_multivariate:
              raise MultivariatePolynomialError(
                  "must be a univariate polynomial")
  
          N = sympify(n)
  
          if N.is_Integer and N >= 1:
              n = int(N)
          else:
              raise ValueError("'n' must an integer and n >= 1, got %s" % n)
  
          x = f.gen
          t = Dummy('t')
  
          r = f.resultant(f.__class__.from_expr(x**n - t, x, t))
  
          return r.replace(t, x)
  
      def cancel(f, g, include=False):
          """
          Cancel common factors in a rational function ``f/g``.
  
          Examples
          ========
  
          >>> from sympy import Poly
          >>> from sympy.abc import x
  
          >>> Poly(2*x**2 - 2, x).cancel(Poly(x**2 - 2*x + 1, x))
          (1, Poly(2*x + 2, x, domain='ZZ'), Poly(x - 1, x, domain='ZZ'))
  
          >>> Poly(2*x**2 - 2, x).cancel(Poly(x**2 - 2*x + 1, x), include=True)
          (Poly(2*x + 2, x, domain='ZZ'), Poly(x - 1, x, domain='ZZ'))
  
          """
          dom, per, F, G = f._unify(g)
  
          if hasattr(F, 'cancel'):
              result = F.cancel(G, include=include)
          else:  # pragma: no cover
              raise OperationNotSupported(f, 'cancel')
  
          if not include:
              if dom.has_assoc_Ring:
                  dom = dom.get_ring()
  
              cp, cq, p, q = result
  
              cp = dom.to_sympy(cp)
              cq = dom.to_sympy(cq)
  
              return cp/cq, per(p), per(q)
          else:
              return tuple(map(per, result))
  
      @property
      def is_zero(f):
          """
          Returns ``True`` if ``f`` is a zero polynomial.
  
          Examples
          ========
  
          >>> from sympy import Poly
          >>> from sympy.abc import x
  
          >>> Poly(0, x).is_zero
          True
          >>> Poly(1, x).is_zero
          False
  
          """
          return f.rep.is_zero
  
      @property
      def is_one(f):
          """
          Returns ``True`` if ``f`` is a unit polynomial.
  
          Examples
          ========
  
          >>> from sympy import Poly
          >>> from sympy.abc import x
  
          >>> Poly(0, x).is_one
          False
          >>> Poly(1, x).is_one
          True
  
          """
          return f.rep.is_one
  
      @property
      def is_sqf(f):
          """
          Returns ``True`` if ``f`` is a square-free polynomial.
  
          Examples
          ========
  
          >>> from sympy import Poly
          >>> from sympy.abc import x
  
          >>> Poly(x**2 - 2*x + 1, x).is_sqf
          False
          >>> Poly(x**2 - 1, x).is_sqf
          True
  
          """
          return f.rep.is_sqf
  
      @property
      def is_monic(f):
          """
          Returns ``True`` if the leading coefficient of ``f`` is one.
  
          Examples
          ========
  
          >>> from sympy import Poly
          >>> from sympy.abc import x
  
          >>> Poly(x + 2, x).is_monic
          True
          >>> Poly(2*x + 2, x).is_monic
          False
  
          """
          return f.rep.is_monic
  
      @property
      def is_primitive(f):
          """
          Returns ``True`` if GCD of the coefficients of ``f`` is one.
  
          Examples
          ========
  
          >>> from sympy import Poly
          >>> from sympy.abc import x
  
          >>> Poly(2*x**2 + 6*x + 12, x).is_primitive
          False
          >>> Poly(x**2 + 3*x + 6, x).is_primitive
          True
  
          """
          return f.rep.is_primitive
  
      @property
      def is_ground(f):
          """
          Returns ``True`` if ``f`` is an element of the ground domain.
  
          Examples
          ========
  
          >>> from sympy import Poly
          >>> from sympy.abc import x, y
  
          >>> Poly(x, x).is_ground
          False
          >>> Poly(2, x).is_ground
          True
          >>> Poly(y, x).is_ground
          True
  
          """
          return f.rep.is_ground
  
      @property
      def is_linear(f):
          """
          Returns ``True`` if ``f`` is linear in all its variables.
  
          Examples
          ========
  
          >>> from sympy import Poly
          >>> from sympy.abc import x, y
  
          >>> Poly(x + y + 2, x, y).is_linear
          True
          >>> Poly(x*y + 2, x, y).is_linear
          False
  
          """
          return f.rep.is_linear
  
      @property
      def is_quadratic(f):
          """
          Returns ``True`` if ``f`` is quadratic in all its variables.
  
          Examples
          ========
  
          >>> from sympy import Poly
          >>> from sympy.abc import x, y
  
          >>> Poly(x*y + 2, x, y).is_quadratic
          True
          >>> Poly(x*y**2 + 2, x, y).is_quadratic
          False
  
          """
          return f.rep.is_quadratic
  
      @property
      def is_monomial(f):
          """
          Returns ``True`` if ``f`` is zero or has only one term.
  
          Examples
          ========
  
          >>> from sympy import Poly
          >>> from sympy.abc import x
  
          >>> Poly(3*x**2, x).is_monomial
          True
          >>> Poly(3*x**2 + 1, x).is_monomial
          False
  
          """
          return f.rep.is_monomial
  
      @property
      def is_homogeneous(f):
          """
          Returns ``True`` if ``f`` is a homogeneous polynomial.
  
          A homogeneous polynomial is a polynomial whose all monomials with
          non-zero coefficients have the same total degree. If you want not
          only to check if a polynomial is homogeneous but also compute its
          homogeneous order, then use :func:`Poly.homogeneous_order`.
  
          Examples
          ========
  
          >>> from sympy import Poly
          >>> from sympy.abc import x, y
  
          >>> Poly(x**2 + x*y, x, y).is_homogeneous
          True
          >>> Poly(x**3 + x*y, x, y).is_homogeneous
          False
  
          """
          return f.rep.is_homogeneous
  
      @property
      def is_irreducible(f):
          """
          Returns ``True`` if ``f`` has no factors over its domain.
  
          Examples
          ========
  
          >>> from sympy import Poly
          >>> from sympy.abc import x
  
          >>> Poly(x**2 + x + 1, x, modulus=2).is_irreducible
          True
          >>> Poly(x**2 + 1, x, modulus=2).is_irreducible
          False
  
          """
          return f.rep.is_irreducible
  
      @property
      def is_univariate(f):
          """
          Returns ``True`` if ``f`` is a univariate polynomial.
  
          Examples
          ========
  
          >>> from sympy import Poly
          >>> from sympy.abc import x, y
  
          >>> Poly(x**2 + x + 1, x).is_univariate
          True
          >>> Poly(x*y**2 + x*y + 1, x, y).is_univariate
          False
          >>> Poly(x*y**2 + x*y + 1, x).is_univariate
          True
          >>> Poly(x**2 + x + 1, x, y).is_univariate
          False
  
          """
          return len(f.gens) == 1
  
      @property
      def is_multivariate(f):
          """
          Returns ``True`` if ``f`` is a multivariate polynomial.
  
          Examples
          ========
  
          >>> from sympy import Poly
          >>> from sympy.abc import x, y
  
          >>> Poly(x**2 + x + 1, x).is_multivariate
          False
          >>> Poly(x*y**2 + x*y + 1, x, y).is_multivariate
          True
          >>> Poly(x*y**2 + x*y + 1, x).is_multivariate
          False
          >>> Poly(x**2 + x + 1, x, y).is_multivariate
          True
  
          """
          return len(f.gens) != 1
  
      @property
      def is_cyclotomic(f):
          """
          Returns ``True`` if ``f`` is a cyclotomic polnomial.
  
          Examples
          ========
  
          >>> from sympy import Poly
          >>> from sympy.abc import x
  
          >>> f = x**16 + x**14 - x**10 + x**8 - x**6 + x**2 + 1
  
          >>> Poly(f).is_cyclotomic
          False
  
          >>> g = x**16 + x**14 - x**10 - x**8 - x**6 + x**2 + 1
  
          >>> Poly(g).is_cyclotomic
          True
  
          """
          return f.rep.is_cyclotomic
  
      def __abs__(f):
          return f.abs()
  
      def __neg__(f):
          return f.neg()
  
      @_polifyit
      def __add__(f, g):
          return f.add(g)
  
      @_polifyit
      def __radd__(f, g):
          return g.add(f)
  
      @_polifyit
      def __sub__(f, g):
          return f.sub(g)
  
      @_polifyit
      def __rsub__(f, g):
          return g.sub(f)
  
      @_polifyit
      def __mul__(f, g):
          return f.mul(g)
  
      @_polifyit
      def __rmul__(f, g):
          return g.mul(f)
  
      @_sympifyit('n', NotImplemented)
      def __pow__(f, n):
          if n.is_Integer and n >= 0:
              return f.pow(n)
          else:
              return NotImplemented
  
      @_polifyit
      def __divmod__(f, g):
          return f.div(g)
  
      @_polifyit
      def __rdivmod__(f, g):
          return g.div(f)
  
      @_polifyit
      def __mod__(f, g):
          return f.rem(g)
  
      @_polifyit
      def __rmod__(f, g):
          return g.rem(f)
  
      @_polifyit
      def __floordiv__(f, g):
          return f.quo(g)
  
      @_polifyit
      def __rfloordiv__(f, g):
          return g.quo(f)
  
      @_sympifyit('g', NotImplemented)
      def __div__(f, g):
          return f.as_expr()/g.as_expr()
  
      @_sympifyit('g', NotImplemented)
      def __rdiv__(f, g):
          return g.as_expr()/f.as_expr()
  
      __truediv__ = __div__
      __rtruediv__ = __rdiv__
  
      @_sympifyit('other', NotImplemented)
      def __eq__(self, other):
          f, g = self, other
  
          if not g.is_Poly:
              try:
                  g = f.__class__(g, f.gens, domain=f.get_domain())
              except (PolynomialError, DomainError, CoercionFailed):
                  return False
  
          if f.gens != g.gens:
              return False
  
          if f.rep.dom != g.rep.dom:
              return False
  
          return f.rep == g.rep
  
      @_sympifyit('g', NotImplemented)
      def __ne__(f, g):
          return not f == g
  
      def __nonzero__(f):
          return not f.is_zero
  
      __bool__ = __nonzero__
  
      def eq(f, g, strict=False):
          if not strict:
              return f == g
          else:
              return f._strict_eq(sympify(g))
  
      def ne(f, g, strict=False):
          return not f.eq(g, strict=strict)
  
      def _strict_eq(f, g):
          return isinstance(g, f.__class__) and f.gens == g.gens and f.rep.eq(g.rep, strict=True)
  
  
  @public
  class PurePoly(Poly):
      """Class for representing pure polynomials. """
  
      def _hashable_content(self):
          """Allow SymPy to hash Poly instances. """
          return (self.rep,)
  
      def __hash__(self):
          return super(PurePoly, self).__hash__()
  
      @property
      def free_symbols(self):
          """
          Free symbols of a polynomial.
  
          Examples
          ========
  
          >>> from sympy import PurePoly
          >>> from sympy.abc import x, y
  
          >>> PurePoly(x**2 + 1).free_symbols
          set()
          >>> PurePoly(x**2 + y).free_symbols
          set()
          >>> PurePoly(x**2 + y, x).free_symbols
          {y}
  
          """
          return self.free_symbols_in_domain
  
      @_sympifyit('other', NotImplemented)
      def __eq__(self, other):
          f, g = self, other
  
          if not g.is_Poly:
              try:
                  g = f.__class__(g, f.gens, domain=f.get_domain())
              except (PolynomialError, DomainError, CoercionFailed):
                  return False
  
          if len(f.gens) != len(g.gens):
              return False
  
          if f.rep.dom != g.rep.dom:
              try:
                  dom = f.rep.dom.unify(g.rep.dom, f.gens)
              except UnificationFailed:
                  return False
  
              f = f.set_domain(dom)
              g = g.set_domain(dom)
  
          return f.rep == g.rep
  
      def _strict_eq(f, g):
          return isinstance(g, f.__class__) and f.rep.eq(g.rep, strict=True)
  
      def _unify(f, g):
          g = sympify(g)
  
          if not g.is_Poly:
              try:
                  return f.rep.dom, f.per, f.rep, f.rep.per(f.rep.dom.from_sympy(g))
              except CoercionFailed:
                  raise UnificationFailed("can't unify %s with %s" % (f, g))
  
          if len(f.gens) != len(g.gens):
              raise UnificationFailed("can't unify %s with %s" % (f, g))
  
          if not (isinstance(f.rep, DMP) and isinstance(g.rep, DMP)):
              raise UnificationFailed("can't unify %s with %s" % (f, g))
  
          cls = f.__class__
          gens = f.gens
  
          dom = f.rep.dom.unify(g.rep.dom, gens)
  
          F = f.rep.convert(dom)
          G = g.rep.convert(dom)
  
          def per(rep, dom=dom, gens=gens, remove=None):
              if remove is not None:
                  gens = gens[:remove] + gens[remove + 1:]
  
                  if not gens:
                      return dom.to_sympy(rep)
  
              return cls.new(rep, *gens)
  
          return dom, per, F, G
  
  
  @public
  def poly_from_expr(expr, *gens, **args):
      """Construct a polynomial from an expression. """
      opt = options.build_options(gens, args)
      return _poly_from_expr(expr, opt)
  
  
  def _poly_from_expr(expr, opt):
      """Construct a polynomial from an expression. """
      orig, expr = expr, sympify(expr)
  
      if not isinstance(expr, Basic):
          raise PolificationFailed(opt, orig, expr)
      elif expr.is_Poly:
          poly = expr.__class__._from_poly(expr, opt)
  
          opt.gens = poly.gens
          opt.domain = poly.domain
  
          if opt.polys is None:
              opt.polys = True
  
          return poly, opt
      elif opt.expand:
          expr = expr.expand()
  
      rep, opt = _dict_from_expr(expr, opt)
      if not opt.gens:
          raise PolificationFailed(opt, orig, expr)
  
      monoms, coeffs = list(zip(*list(rep.items())))
      domain = opt.domain
  
      if domain is None:
          opt.domain, coeffs = construct_domain(coeffs, opt=opt)
      else:
          coeffs = list(map(domain.from_sympy, coeffs))
  
      rep = dict(list(zip(monoms, coeffs)))
      poly = Poly._from_dict(rep, opt)
  
      if opt.polys is None:
          opt.polys = False
  
      return poly, opt
  
  
  @public
  def parallel_poly_from_expr(exprs, *gens, **args):
      """Construct polynomials from expressions. """
      opt = options.build_options(gens, args)
      return _parallel_poly_from_expr(exprs, opt)
  
  
  def _parallel_poly_from_expr(exprs, opt):
      """Construct polynomials from expressions. """
      from sympy.functions.elementary.piecewise import Piecewise
  
      if len(exprs) == 2:
          f, g = exprs
  
          if isinstance(f, Poly) and isinstance(g, Poly):
              f = f.__class__._from_poly(f, opt)
              g = g.__class__._from_poly(g, opt)
  
              f, g = f.unify(g)
  
              opt.gens = f.gens
              opt.domain = f.domain
  
              if opt.polys is None:
                  opt.polys = True
  
              return [f, g], opt
  
      origs, exprs = list(exprs), []
      _exprs, _polys = [], []
  
      failed = False
  
      for i, expr in enumerate(origs):
          expr = sympify(expr)
  
          if isinstance(expr, Basic):
              if expr.is_Poly:
                  _polys.append(i)
              else:
                  _exprs.append(i)
  
                  if opt.expand:
                      expr = expr.expand()
          else:
              failed = True
  
          exprs.append(expr)
  
      if failed:
          raise PolificationFailed(opt, origs, exprs, True)
  
      if _polys:
          # XXX: this is a temporary solution
          for i in _polys:
              exprs[i] = exprs[i].as_expr()
  
      reps, opt = _parallel_dict_from_expr(exprs, opt)
      if not opt.gens:
          raise PolificationFailed(opt, origs, exprs, True)
  
      for k in opt.gens:
          if isinstance(k, Piecewise):
              raise PolynomialError("Piecewise generators do not make sense")
  
      coeffs_list, lengths = [], []
  
      all_monoms = []
      all_coeffs = []
  
      for rep in reps:
          monoms, coeffs = list(zip(*list(rep.items())))
  
          coeffs_list.extend(coeffs)
          all_monoms.append(monoms)
  
          lengths.append(len(coeffs))
  
      domain = opt.domain
  
      if domain is None:
          opt.domain, coeffs_list = construct_domain(coeffs_list, opt=opt)
      else:
          coeffs_list = list(map(domain.from_sympy, coeffs_list))
  
      for k in lengths:
          all_coeffs.append(coeffs_list[:k])
          coeffs_list = coeffs_list[k:]
  
      polys = []
  
      for monoms, coeffs in zip(all_monoms, all_coeffs):
          rep = dict(list(zip(monoms, coeffs)))
          poly = Poly._from_dict(rep, opt)
          polys.append(poly)
  
      if opt.polys is None:
          opt.polys = bool(_polys)
  
      return polys, opt
  
  
  def _update_args(args, key, value):
      """Add a new ``(key, value)`` pair to arguments ``dict``. """
      args = dict(args)
  
      if key not in args:
          args[key] = value
  
      return args
  
  
  @public
  def degree(f, gen=0):
      """
      Return the degree of ``f`` in the given variable.
  
      The degree of 0 is negative infinity.
  
      Examples
      ========
  
      >>> from sympy import degree
      >>> from sympy.abc import x, y
  
      >>> degree(x**2 + y*x + 1, gen=x)
      2
      >>> degree(x**2 + y*x + 1, gen=y)
      1
      >>> degree(0, x)
      -oo
  
      See also
      ========
  
      sympy.polys.polytools.Poly.total_degree
      degree_list
      """
      if type(f)==type((x+y)**2) or type(f)==type(20*(x**40 + y)**10):#to match types <class 'sympy.core.power.Pow'> and 
                                                                  #<class 'sympy.core.mul.Mul'>
      inside=str(f)
      gen1=gen
      result = inside.rindex('**')
      l=inside[:result]
      exponent=int(inside[result+2:])
      expression=parse_expr(l)
  
      return degree(expression, gen=gen1)*exponent
  
      f = sympify(f, strict=True)
      gen_is_Num = sympify(gen, strict=True).is_Number
      if f.is_Poly:
          p = f
          isNum = p.as_expr().is_Number
      else:
          isNum = f.is_Number
          if not isNum:
              if gen_is_Num:
                  p, _ = poly_from_expr(f)
              else:
                  p, _ = poly_from_expr(f, gen)
  
      if isNum:
          return S.Zero if f else S.NegativeInfinity
  
      if not gen_is_Num:
          if f.is_Poly and gen not in p.gens:
              # try recast without explicit gens
              p, _ = poly_from_expr(f.as_expr())
          if gen not in p.gens:
              return S.Zero
      elif not f.is_Poly and len(f.free_symbols) > 1:
          raise TypeError(filldedent('''
           A symbolic generator of interest is required for a multivariate
           expression like func = %s, e.g. degree(func, gen = %s) instead of
           degree(func, gen = %s).
          ''' % (f, next(ordered(f.free_symbols)), gen)))
  
      return Integer(p.degree(gen))
  
  
  @public
  def total_degree(f, *gens):
      """
      Return the total_degree of ``f`` in the given variables.
  
      Examples
      ========
      >>> from sympy import total_degree, Poly
      >>> from sympy.abc import x, y, z
  
      >>> total_degree(1)
      0
      >>> total_degree(x + x*y)
      2
      >>> total_degree(x + x*y, x)
      1
  
      If the expression is a Poly and no variables are given
      then the generators of the Poly will be used:
  
      >>> p = Poly(x + x*y, y)
      >>> total_degree(p)
      1
  
      To deal with the underlying expression of the Poly, convert
      it to an Expr:
  
      >>> total_degree(p.as_expr())
      2
  
      This is done automatically if any variables are given:
  
      >>> total_degree(p, x)
      1
  
      See also
      ========
      degree
      """
  
      p = sympify(f)
      if p.is_Poly:
          p = p.as_expr()
      if p.is_Number:
          rv = 0
      else:
          if f.is_Poly:
              gens = gens or f.gens
          rv = Poly(p, gens).total_degree()
  
      return Integer(rv)
  
  
  @public
  def degree_list(f, *gens, **args):
      """
      Return a list of degrees of ``f`` in all variables.
  
      Examples
      ========
  
      >>> from sympy import degree_list
      >>> from sympy.abc import x, y
  
      >>> degree_list(x**2 + y*x + 1)
      (2, 1)
  
      """
      options.allowed_flags(args, ['polys'])
  
      try:
          F, opt = poly_from_expr(f, *gens, **args)
      except PolificationFailed as exc:
          raise ComputationFailed('degree_list', 1, exc)
  
      degrees = F.degree_list()
  
      return tuple(map(Integer, degrees))
  
  
  @public
  def LC(f, *gens, **args):
      """
      Return the leading coefficient of ``f``.
  
      Examples
      ========
  
      >>> from sympy import LC
      >>> from sympy.abc import x, y
  
      >>> LC(4*x**2 + 2*x*y**2 + x*y + 3*y)
      4
  
      """
      options.allowed_flags(args, ['polys'])
  
      try:
          F, opt = poly_from_expr(f, *gens, **args)
      except PolificationFailed as exc:
          raise ComputationFailed('LC', 1, exc)
  
      return F.LC(order=opt.order)
  
  
  @public
  def LM(f, *gens, **args):
      """
      Return the leading monomial of ``f``.
  
      Examples
      ========
  
      >>> from sympy import LM
      >>> from sympy.abc import x, y
  
      >>> LM(4*x**2 + 2*x*y**2 + x*y + 3*y)
      x**2
  
      """
      options.allowed_flags(args, ['polys'])
  
      try:
          F, opt = poly_from_expr(f, *gens, **args)
      except PolificationFailed as exc:
          raise ComputationFailed('LM', 1, exc)
  
      monom = F.LM(order=opt.order)
      return monom.as_expr()
  
  
  @public
  def LT(f, *gens, **args):
      """
      Return the leading term of ``f``.
  
      Examples
      ========
  
      >>> from sympy import LT
      >>> from sympy.abc import x, y
  
      >>> LT(4*x**2 + 2*x*y**2 + x*y + 3*y)
      4*x**2
  
      """
      options.allowed_flags(args, ['polys'])
  
      try:
          F, opt = poly_from_expr(f, *gens, **args)
      except PolificationFailed as exc:
          raise ComputationFailed('LT', 1, exc)
  
      monom, coeff = F.LT(order=opt.order)
      return coeff*monom.as_expr()
  
  
  @public
  def pdiv(f, g, *gens, **args):
      """
      Compute polynomial pseudo-division of ``f`` and ``g``.
  
      Examples
      ========
  
      >>> from sympy import pdiv
      >>> from sympy.abc import x
  
      >>> pdiv(x**2 + 1, 2*x - 4)
      (2*x + 4, 20)
  
      """
      options.allowed_flags(args, ['polys'])
  
      try:
          (F, G), opt = parallel_poly_from_expr((f, g), *gens, **args)
      except PolificationFailed as exc:
          raise ComputationFailed('pdiv', 2, exc)
  
      q, r = F.pdiv(G)
  
      if not opt.polys:
          return q.as_expr(), r.as_expr()
      else:
          return q, r
  
  
  @public
  def prem(f, g, *gens, **args):
      """
      Compute polynomial pseudo-remainder of ``f`` and ``g``.
  
      Examples
      ========
  
      >>> from sympy import prem
      >>> from sympy.abc import x
  
      >>> prem(x**2 + 1, 2*x - 4)
      20
  
      """
      options.allowed_flags(args, ['polys'])
  
      try:
          (F, G), opt = parallel_poly_from_expr((f, g), *gens, **args)
      except PolificationFailed as exc:
          raise ComputationFailed('prem', 2, exc)
  
      r = F.prem(G)
  
      if not opt.polys:
          return r.as_expr()
      else:
          return r
  
  
  @public
  def pquo(f, g, *gens, **args):
      """
      Compute polynomial pseudo-quotient of ``f`` and ``g``.
  
      Examples
      ========
  
      >>> from sympy import pquo
      >>> from sympy.abc import x
  
      >>> pquo(x**2 + 1, 2*x - 4)
      2*x + 4
      >>> pquo(x**2 - 1, 2*x - 1)
      2*x + 1
  
      """
      options.allowed_flags(args, ['polys'])
  
      try:
          (F, G), opt = parallel_poly_from_expr((f, g), *gens, **args)
      except PolificationFailed as exc:
          raise ComputationFailed('pquo', 2, exc)
  
      try:
          q = F.pquo(G)
      except ExactQuotientFailed:
          raise ExactQuotientFailed(f, g)
  
      if not opt.polys:
          return q.as_expr()
      else:
          return q
  
  
  @public
  def pexquo(f, g, *gens, **args):
      """
      Compute polynomial exact pseudo-quotient of ``f`` and ``g``.
  
      Examples
      ========
  
      >>> from sympy import pexquo
      >>> from sympy.abc import x
  
      >>> pexquo(x**2 - 1, 2*x - 2)
      2*x + 2
  
      >>> pexquo(x**2 + 1, 2*x - 4)
      Traceback (most recent call last):
      ...
      ExactQuotientFailed: 2*x - 4 does not divide x**2 + 1
  
      """
      options.allowed_flags(args, ['polys'])
  
      try:
          (F, G), opt = parallel_poly_from_expr((f, g), *gens, **args)
      except PolificationFailed as exc:
          raise ComputationFailed('pexquo', 2, exc)
  
      q = F.pexquo(G)
  
      if not opt.polys:
          return q.as_expr()
      else:
          return q
  
  
  @public
  def div(f, g, *gens, **args):
      """
      Compute polynomial division of ``f`` and ``g``.
  
      Examples
      ========
  
      >>> from sympy import div, ZZ, QQ
      >>> from sympy.abc import x
  
      >>> div(x**2 + 1, 2*x - 4, domain=ZZ)
      (0, x**2 + 1)
      >>> div(x**2 + 1, 2*x - 4, domain=QQ)
      (x/2 + 1, 5)
  
      """
      options.allowed_flags(args, ['auto', 'polys'])
  
      try:
          (F, G), opt = parallel_poly_from_expr((f, g), *gens, **args)
      except PolificationFailed as exc:
          raise ComputationFailed('div', 2, exc)
  
      q, r = F.div(G, auto=opt.auto)
  
      if not opt.polys:
          return q.as_expr(), r.as_expr()
      else:
          return q, r
  
  
  @public
  def rem(f, g, *gens, **args):
      """
      Compute polynomial remainder of ``f`` and ``g``.
  
      Examples
      ========
  
      >>> from sympy import rem, ZZ, QQ
      >>> from sympy.abc import x
  
      >>> rem(x**2 + 1, 2*x - 4, domain=ZZ)
      x**2 + 1
      >>> rem(x**2 + 1, 2*x - 4, domain=QQ)
      5
  
      """
      options.allowed_flags(args, ['auto', 'polys'])
  
      try:
          (F, G), opt = parallel_poly_from_expr((f, g), *gens, **args)
      except PolificationFailed as exc:
          raise ComputationFailed('rem', 2, exc)
  
      r = F.rem(G, auto=opt.auto)
  
      if not opt.polys:
          return r.as_expr()
      else:
          return r
  
  
  @public
  def quo(f, g, *gens, **args):
      """
      Compute polynomial quotient of ``f`` and ``g``.
  
      Examples
      ========
  
      >>> from sympy import quo
      >>> from sympy.abc import x
  
      >>> quo(x**2 + 1, 2*x - 4)
      x/2 + 1
      >>> quo(x**2 - 1, x - 1)
      x + 1
  
      """
      options.allowed_flags(args, ['auto', 'polys'])
  
      try:
          (F, G), opt = parallel_poly_from_expr((f, g), *gens, **args)
      except PolificationFailed as exc:
          raise ComputationFailed('quo', 2, exc)
  
      q = F.quo(G, auto=opt.auto)
  
      if not opt.polys:
          return q.as_expr()
      else:
          return q
  
  
  @public
  def exquo(f, g, *gens, **args):
      """
      Compute polynomial exact quotient of ``f`` and ``g``.
  
      Examples
      ========
  
      >>> from sympy import exquo
      >>> from sympy.abc import x
  
      >>> exquo(x**2 - 1, x - 1)
      x + 1
  
      >>> exquo(x**2 + 1, 2*x - 4)
      Traceback (most recent call last):
      ...
      ExactQuotientFailed: 2*x - 4 does not divide x**2 + 1
  
      """
      options.allowed_flags(args, ['auto', 'polys'])
  
      try:
          (F, G), opt = parallel_poly_from_expr((f, g), *gens, **args)
      except PolificationFailed as exc:
          raise ComputationFailed('exquo', 2, exc)
  
      q = F.exquo(G, auto=opt.auto)
  
      if not opt.polys:
          return q.as_expr()
      else:
          return q
  
  
  @public
  def half_gcdex(f, g, *gens, **args):
      """
      Half extended Euclidean algorithm of ``f`` and ``g``.
  
      Returns ``(s, h)`` such that ``h = gcd(f, g)`` and ``s*f = h (mod g)``.
  
      Examples
      ========
  
      >>> from sympy import half_gcdex
      >>> from sympy.abc import x
  
      >>> half_gcdex(x**4 - 2*x**3 - 6*x**2 + 12*x + 15, x**3 + x**2 - 4*x - 4)
      (3/5 - x/5, x + 1)
  
      """
      options.allowed_flags(args, ['auto', 'polys'])
  
      try:
          (F, G), opt = parallel_poly_from_expr((f, g), *gens, **args)
      except PolificationFailed as exc:
          domain, (a, b) = construct_domain(exc.exprs)
  
          try:
              s, h = domain.half_gcdex(a, b)
          except NotImplementedError:
              raise ComputationFailed('half_gcdex', 2, exc)
          else:
              return domain.to_sympy(s), domain.to_sympy(h)
  
      s, h = F.half_gcdex(G, auto=opt.auto)
  
      if not opt.polys:
          return s.as_expr(), h.as_expr()
      else:
          return s, h
  
  
  @public
  def gcdex(f, g, *gens, **args):
      """
      Extended Euclidean algorithm of ``f`` and ``g``.
  
      Returns ``(s, t, h)`` such that ``h = gcd(f, g)`` and ``s*f + t*g = h``.
  
      Examples
      ========
  
      >>> from sympy import gcdex
      >>> from sympy.abc import x
  
      >>> gcdex(x**4 - 2*x**3 - 6*x**2 + 12*x + 15, x**3 + x**2 - 4*x - 4)
      (3/5 - x/5, x**2/5 - 6*x/5 + 2, x + 1)
  
      """
      options.allowed_flags(args, ['auto', 'polys'])
  
      try:
          (F, G), opt = parallel_poly_from_expr((f, g), *gens, **args)
      except PolificationFailed as exc:
          domain, (a, b) = construct_domain(exc.exprs)
  
          try:
              s, t, h = domain.gcdex(a, b)
          except NotImplementedError:
              raise ComputationFailed('gcdex', 2, exc)
          else:
              return domain.to_sympy(s), domain.to_sympy(t), domain.to_sympy(h)
  
      s, t, h = F.gcdex(G, auto=opt.auto)
  
      if not opt.polys:
          return s.as_expr(), t.as_expr(), h.as_expr()
      else:
          return s, t, h
  
  
  @public
  def invert(f, g, *gens, **args):
      """
      Invert ``f`` modulo ``g`` when possible.
  
      Examples
      ========
  
      >>> from sympy import invert, S
      >>> from sympy.core.numbers import mod_inverse
      >>> from sympy.abc import x
  
      >>> invert(x**2 - 1, 2*x - 1)
      -4/3
  
      >>> invert(x**2 - 1, x - 1)
      Traceback (most recent call last):
      ...
      NotInvertible: zero divisor
  
      For more efficient inversion of Rationals,
      use the :obj:`~.mod_inverse` function:
  
      >>> mod_inverse(3, 5)
      2
      >>> (S(2)/5).invert(S(7)/3)
      5/2
  
      See Also
      ========
  
      sympy.core.numbers.mod_inverse
  
      """
      options.allowed_flags(args, ['auto', 'polys'])
  
      try:
          (F, G), opt = parallel_poly_from_expr((f, g), *gens, **args)
      except PolificationFailed as exc:
          domain, (a, b) = construct_domain(exc.exprs)
  
          try:
              return domain.to_sympy(domain.invert(a, b))
          except NotImplementedError:
              raise ComputationFailed('invert', 2, exc)
  
      h = F.invert(G, auto=opt.auto)
  
      if not opt.polys:
          return h.as_expr()
      else:
          return h
  
  
  @public
  def subresultants(f, g, *gens, **args):
      """
      Compute subresultant PRS of ``f`` and ``g``.
  
      Examples
      ========
  
      >>> from sympy import subresultants
      >>> from sympy.abc import x
  
      >>> subresultants(x**2 + 1, x**2 - 1)
      [x**2 + 1, x**2 - 1, -2]
  
      """
      options.allowed_flags(args, ['polys'])
  
      try:
          (F, G), opt = parallel_poly_from_expr((f, g), *gens, **args)
      except PolificationFailed as exc:
          raise ComputationFailed('subresultants', 2, exc)
  
      result = F.subresultants(G)
  
      if not opt.polys:
          return [r.as_expr() for r in result]
      else:
          return result
  
  
  @public
  def resultant(f, g, *gens, **args):
      """
      Compute resultant of ``f`` and ``g``.
  
      Examples
      ========
  
      >>> from sympy import resultant
      >>> from sympy.abc import x
  
      >>> resultant(x**2 + 1, x**2 - 1)
      4
  
      """
      includePRS = args.pop('includePRS', False)
      options.allowed_flags(args, ['polys'])
  
      try:
          (F, G), opt = parallel_poly_from_expr((f, g), *gens, **args)
      except PolificationFailed as exc:
          raise ComputationFailed('resultant', 2, exc)
  
      if includePRS:
          result, R = F.resultant(G, includePRS=includePRS)
      else:
          result = F.resultant(G)
  
      if not opt.polys:
          if includePRS:
              return result.as_expr(), [r.as_expr() for r in R]
          return result.as_expr()
      else:
          if includePRS:
              return result, R
          return result
  
  
  @public
  def discriminant(f, *gens, **args):
      """
      Compute discriminant of ``f``.
  
      Examples
      ========
  
      >>> from sympy import discriminant
      >>> from sympy.abc import x
  
      >>> discriminant(x**2 + 2*x + 3)
      -8
  
      """
      options.allowed_flags(args, ['polys'])
  
      try:
          F, opt = poly_from_expr(f, *gens, **args)
      except PolificationFailed as exc:
          raise ComputationFailed('discriminant', 1, exc)
  
      result = F.discriminant()
  
      if not opt.polys:
          return result.as_expr()
      else:
          return result
  
  
  @public
  def cofactors(f, g, *gens, **args):
      """
      Compute GCD and cofactors of ``f`` and ``g``.
  
      Returns polynomials ``(h, cff, cfg)`` such that ``h = gcd(f, g)``, and
      ``cff = quo(f, h)`` and ``cfg = quo(g, h)`` are, so called, cofactors
      of ``f`` and ``g``.
  
      Examples
      ========
  
      >>> from sympy import cofactors
      >>> from sympy.abc import x
  
      >>> cofactors(x**2 - 1, x**2 - 3*x + 2)
      (x - 1, x + 1, x - 2)
  
      """
      options.allowed_flags(args, ['polys'])
  
      try:
          (F, G), opt = parallel_poly_from_expr((f, g), *gens, **args)
      except PolificationFailed as exc:
          domain, (a, b) = construct_domain(exc.exprs)
  
          try:
              h, cff, cfg = domain.cofactors(a, b)
          except NotImplementedError:
              raise ComputationFailed('cofactors', 2, exc)
          else:
              return domain.to_sympy(h), domain.to_sympy(cff), domain.to_sympy(cfg)
  
      h, cff, cfg = F.cofactors(G)
  
      if not opt.polys:
          return h.as_expr(), cff.as_expr(), cfg.as_expr()
      else:
          return h, cff, cfg
  
  
  @public
  def gcd_list(seq, *gens, **args):
      """
      Compute GCD of a list of polynomials.
  
      Examples
      ========
  
      >>> from sympy import gcd_list
      >>> from sympy.abc import x
  
      >>> gcd_list([x**3 - 1, x**2 - 1, x**2 - 3*x + 2])
      x - 1
  
      """
      seq = sympify(seq)
  
      def try_non_polynomial_gcd(seq):
          if not gens and not args:
              domain, numbers = construct_domain(seq)
  
              if not numbers:
                  return domain.zero
              elif domain.is_Numerical:
                  result, numbers = numbers[0], numbers[1:]
  
                  for number in numbers:
                      result = domain.gcd(result, number)
  
                      if domain.is_one(result):
                          break
  
                  return domain.to_sympy(result)
  
          return None
  
      result = try_non_polynomial_gcd(seq)
  
      if result is not None:
          return result
  
      options.allowed_flags(args, ['polys'])
  
      try:
          polys, opt = parallel_poly_from_expr(seq, *gens, **args)
  
          # gcd for domain Q[irrational] (purely algebraic irrational)
          if len(seq) > 1 and all(elt.is_algebraic and elt.is_irrational for elt in seq):
              a = seq[-1]
              lst = [ (a/elt).ratsimp() for elt in seq[:-1] ]
              if all(frc.is_rational for frc in lst):
                  lc = 1
                  for frc in lst:
                      lc = lcm(lc, frc.as_numer_denom()[0])
                  return a/lc
  
      except PolificationFailed as exc:
          result = try_non_polynomial_gcd(exc.exprs)
  
          if result is not None:
              return result
          else:
              raise ComputationFailed('gcd_list', len(seq), exc)
  
      if not polys:
          if not opt.polys:
              return S.Zero
          else:
              return Poly(0, opt=opt)
  
      result, polys = polys[0], polys[1:]
  
      for poly in polys:
          result = result.gcd(poly)
  
          if result.is_one:
              break
  
      if not opt.polys:
          return result.as_expr()
      else:
          return result
  
  
  @public
  def gcd(f, g=None, *gens, **args):
      """
      Compute GCD of ``f`` and ``g``.
  
      Examples
      ========
  
      >>> from sympy import gcd
      >>> from sympy.abc import x
  
      >>> gcd(x**2 - 1, x**2 - 3*x + 2)
      x - 1
  
      """
      if hasattr(f, '__iter__'):
          if g is not None:
              gens = (g,) + gens
  
          return gcd_list(f, *gens, **args)
      elif g is None:
          raise TypeError("gcd() takes 2 arguments or a sequence of arguments")
  
      options.allowed_flags(args, ['polys'])
  
      try:
          (F, G), opt = parallel_poly_from_expr((f, g), *gens, **args)
  
          # gcd for domain Q[irrational] (purely algebraic irrational)
          a, b = map(sympify, (f, g))
          if a.is_algebraic and a.is_irrational and b.is_algebraic and b.is_irrational:
              frc = (a/b).ratsimp()
              if frc.is_rational:
                  return a/frc.as_numer_denom()[0]
  
      except PolificationFailed as exc:
          domain, (a, b) = construct_domain(exc.exprs)
  
          try:
              return domain.to_sympy(domain.gcd(a, b))
          except NotImplementedError:
              raise ComputationFailed('gcd', 2, exc)
  
      result = F.gcd(G)
  
      if not opt.polys:
          return result.as_expr()
      else:
          return result
  
  
  @public
  def lcm_list(seq, *gens, **args):
      """
      Compute LCM of a list of polynomials.
  
      Examples
      ========
  
      >>> from sympy import lcm_list
      >>> from sympy.abc import x
  
      >>> lcm_list([x**3 - 1, x**2 - 1, x**2 - 3*x + 2])
      x**5 - x**4 - 2*x**3 - x**2 + x + 2
  
      """
      seq = sympify(seq)
  
      def try_non_polynomial_lcm(seq):
          if not gens and not args:
              domain, numbers = construct_domain(seq)
  
              if not numbers:
                  return domain.one
              elif domain.is_Numerical:
                  result, numbers = numbers[0], numbers[1:]
  
                  for number in numbers:
                      result = domain.lcm(result, number)
  
                  return domain.to_sympy(result)
  
          return None
  
      result = try_non_polynomial_lcm(seq)
  
      if result is not None:
          return result
  
      options.allowed_flags(args, ['polys'])
  
      try:
          polys, opt = parallel_poly_from_expr(seq, *gens, **args)
  
          # lcm for domain Q[irrational] (purely algebraic irrational)
          if len(seq) > 1 and all(elt.is_algebraic and elt.is_irrational for elt in seq):
              a = seq[-1]
              lst = [ (a/elt).ratsimp() for elt in seq[:-1] ]
              if all(frc.is_rational for frc in lst):
                  lc = 1
                  for frc in lst:
                      lc = lcm(lc, frc.as_numer_denom()[1])
                  return a*lc
  
      except PolificationFailed as exc:
          result = try_non_polynomial_lcm(exc.exprs)
  
          if result is not None:
              return result
          else:
              raise ComputationFailed('lcm_list', len(seq), exc)
  
      if not polys:
          if not opt.polys:
              return S.One
          else:
              return Poly(1, opt=opt)
  
      result, polys = polys[0], polys[1:]
  
      for poly in polys:
          result = result.lcm(poly)
  
      if not opt.polys:
          return result.as_expr()
      else:
          return result
  
  
  @public
  def lcm(f, g=None, *gens, **args):
      """
      Compute LCM of ``f`` and ``g``.
  
      Examples
      ========
  
      >>> from sympy import lcm
      >>> from sympy.abc import x
  
      >>> lcm(x**2 - 1, x**2 - 3*x + 2)
      x**3 - 2*x**2 - x + 2
  
      """
      if hasattr(f, '__iter__'):
          if g is not None:
              gens = (g,) + gens
  
          return lcm_list(f, *gens, **args)
      elif g is None:
          raise TypeError("lcm() takes 2 arguments or a sequence of arguments")
  
      options.allowed_flags(args, ['polys'])
  
      try:
          (F, G), opt = parallel_poly_from_expr((f, g), *gens, **args)
  
          # lcm for domain Q[irrational] (purely algebraic irrational)
          a, b = map(sympify, (f, g))
          if a.is_algebraic and a.is_irrational and b.is_algebraic and b.is_irrational:
              frc = (a/b).ratsimp()
              if frc.is_rational:
                  return a*frc.as_numer_denom()[1]
  
      except PolificationFailed as exc:
          domain, (a, b) = construct_domain(exc.exprs)
  
          try:
              return domain.to_sympy(domain.lcm(a, b))
          except NotImplementedError:
              raise ComputationFailed('lcm', 2, exc)
  
      result = F.lcm(G)
  
      if not opt.polys:
          return result.as_expr()
      else:
          return result
  
  
  @public
  def terms_gcd(f, *gens, **args):
      """
      Remove GCD of terms from ``f``.
  
      If the ``deep`` flag is True, then the arguments of ``f`` will have
      terms_gcd applied to them.
  
      If a fraction is factored out of ``f`` and ``f`` is an Add, then
      an unevaluated Mul will be returned so that automatic simplification
      does not redistribute it. The hint ``clear``, when set to False, can be
      used to prevent such factoring when all coefficients are not fractions.
  
      Examples
      ========
  
      >>> from sympy import terms_gcd, cos
      >>> from sympy.abc import x, y
      >>> terms_gcd(x**6*y**2 + x**3*y, x, y)
      x**3*y*(x**3*y + 1)
  
      The default action of polys routines is to expand the expression
      given to them. terms_gcd follows this behavior:
  
      >>> terms_gcd((3+3*x)*(x+x*y))
      3*x*(x*y + x + y + 1)
  
      If this is not desired then the hint ``expand`` can be set to False.
      In this case the expression will be treated as though it were comprised
      of one or more terms:
  
      >>> terms_gcd((3+3*x)*(x+x*y), expand=False)
      (3*x + 3)*(x*y + x)
  
      In order to traverse factors of a Mul or the arguments of other
      functions, the ``deep`` hint can be used:
  
      >>> terms_gcd((3 + 3*x)*(x + x*y), expand=False, deep=True)
      3*x*(x + 1)*(y + 1)
      >>> terms_gcd(cos(x + x*y), deep=True)
      cos(x*(y + 1))
  
      Rationals are factored out by default:
  
      >>> terms_gcd(x + y/2)
      (2*x + y)/2
  
      Only the y-term had a coefficient that was a fraction; if one
      does not want to factor out the 1/2 in cases like this, the
      flag ``clear`` can be set to False:
  
      >>> terms_gcd(x + y/2, clear=False)
      x + y/2
      >>> terms_gcd(x*y/2 + y**2, clear=False)
      y*(x/2 + y)
  
      The ``clear`` flag is ignored if all coefficients are fractions:
  
      >>> terms_gcd(x/3 + y/2, clear=False)
      (2*x + 3*y)/6
  
      See Also
      ========
      sympy.core.exprtools.gcd_terms, sympy.core.exprtools.factor_terms
  
      """
      from sympy.core.relational import Equality
  
      orig = sympify(f)
  
      if isinstance(f, Equality):
          return Equality(*(terms_gcd(s, *gens, **args) for s in [f.lhs, f.rhs]))
      elif isinstance(f, Relational):
          raise TypeError("Inequalities can not be used with terms_gcd. Found: %s" %(f,))
  
      if not isinstance(f, Expr) or f.is_Atom:
          return orig
  
      if args.get('deep', False):
          new = f.func(*[terms_gcd(a, *gens, **args) for a in f.args])
          args.pop('deep')
          args['expand'] = False
          return terms_gcd(new, *gens, **args)
  
      clear = args.pop('clear', True)
      options.allowed_flags(args, ['polys'])
  
      try:
          F, opt = poly_from_expr(f, *gens, **args)
      except PolificationFailed as exc:
          return exc.expr
  
      J, f = F.terms_gcd()
  
      if opt.domain.is_Ring:
          if opt.domain.is_Field:
              denom, f = f.clear_denoms(convert=True)
  
          coeff, f = f.primitive()
  
          if opt.domain.is_Field:
              coeff /= denom
      else:
          coeff = S.One
  
      term = Mul(*[x**j for x, j in zip(f.gens, J)])
      if coeff == 1:
          coeff = S.One
          if term == 1:
              return orig
  
      if clear:
          return _keep_coeff(coeff, term*f.as_expr())
      # base the clearing on the form of the original expression, not
      # the (perhaps) Mul that we have now
      coeff, f = _keep_coeff(coeff, f.as_expr(), clear=False).as_coeff_Mul()
      return _keep_coeff(coeff, term*f, clear=False)
  
  
  @public
  def trunc(f, p, *gens, **args):
      """
      Reduce ``f`` modulo a constant ``p``.
  
      Examples
      ========
  
      >>> from sympy import trunc
      >>> from sympy.abc import x
  
      >>> trunc(2*x**3 + 3*x**2 + 5*x + 7, 3)
      -x**3 - x + 1
  
      """
      options.allowed_flags(args, ['auto', 'polys'])
  
      try:
          F, opt = poly_from_expr(f, *gens, **args)
      except PolificationFailed as exc:
          raise ComputationFailed('trunc', 1, exc)
  
      result = F.trunc(sympify(p))
  
      if not opt.polys:
          return result.as_expr()
      else:
          return result
  
  
  @public
  def monic(f, *gens, **args):
      """
      Divide all coefficients of ``f`` by ``LC(f)``.
  
      Examples
      ========
  
      >>> from sympy import monic
      >>> from sympy.abc import x
  
      >>> monic(3*x**2 + 4*x + 2)
      x**2 + 4*x/3 + 2/3
  
      """
      options.allowed_flags(args, ['auto', 'polys'])
  
      try:
          F, opt = poly_from_expr(f, *gens, **args)
      except PolificationFailed as exc:
          raise ComputationFailed('monic', 1, exc)
  
      result = F.monic(auto=opt.auto)
  
      if not opt.polys:
          return result.as_expr()
      else:
          return result
  
  
  @public
  def content(f, *gens, **args):
      """
      Compute GCD of coefficients of ``f``.
  
      Examples
      ========
  
      >>> from sympy import content
      >>> from sympy.abc import x
  
      >>> content(6*x**2 + 8*x + 12)
      2
  
      """
      options.allowed_flags(args, ['polys'])
  
      try:
          F, opt = poly_from_expr(f, *gens, **args)
      except PolificationFailed as exc:
          raise ComputationFailed('content', 1, exc)
  
      return F.content()
  
  
  @public
  def primitive(f, *gens, **args):
      """
      Compute content and the primitive form of ``f``.
  
      Examples
      ========
  
      >>> from sympy.polys.polytools import primitive
      >>> from sympy.abc import x
  
      >>> primitive(6*x**2 + 8*x + 12)
      (2, 3*x**2 + 4*x + 6)
  
      >>> eq = (2 + 2*x)*x + 2
  
      Expansion is performed by default:
  
      >>> primitive(eq)
      (2, x**2 + x + 1)
  
      Set ``expand`` to False to shut this off. Note that the
      extraction will not be recursive; use the as_content_primitive method
      for recursive, non-destructive Rational extraction.
  
      >>> primitive(eq, expand=False)
      (1, x*(2*x + 2) + 2)
  
      >>> eq.as_content_primitive()
      (2, x*(x + 1) + 1)
  
      """
      options.allowed_flags(args, ['polys'])
  
      try:
          F, opt = poly_from_expr(f, *gens, **args)
      except PolificationFailed as exc:
          raise ComputationFailed('primitive', 1, exc)
  
      cont, result = F.primitive()
      if not opt.polys:
          return cont, result.as_expr()
      else:
          return cont, result
  
  
  @public
  def compose(f, g, *gens, **args):
      """
      Compute functional composition ``f(g)``.
  
      Examples
      ========
  
      >>> from sympy import compose
      >>> from sympy.abc import x
  
      >>> compose(x**2 + x, x - 1)
      x**2 - x
  
      """
      options.allowed_flags(args, ['polys'])
  
      try:
          (F, G), opt = parallel_poly_from_expr((f, g), *gens, **args)
      except PolificationFailed as exc:
          raise ComputationFailed('compose', 2, exc)
  
      result = F.compose(G)
  
      if not opt.polys:
          return result.as_expr()
      else:
          return result
  
  
  @public
  def decompose(f, *gens, **args):
      """
      Compute functional decomposition of ``f``.
  
      Examples
      ========
  
      >>> from sympy import decompose
      >>> from sympy.abc import x
  
      >>> decompose(x**4 + 2*x**3 - x - 1)
      [x**2 - x - 1, x**2 + x]
  
      """
      options.allowed_flags(args, ['polys'])
  
      try:
          F, opt = poly_from_expr(f, *gens, **args)
      except PolificationFailed as exc:
          raise ComputationFailed('decompose', 1, exc)
  
      result = F.decompose()
  
      if not opt.polys:
          return [r.as_expr() for r in result]
      else:
          return result
  
  
  @public
  def sturm(f, *gens, **args):
      """
      Compute Sturm sequence of ``f``.
  
      Examples
      ========
  
      >>> from sympy import sturm
      >>> from sympy.abc import x
  
      >>> sturm(x**3 - 2*x**2 + x - 3)
      [x**3 - 2*x**2 + x - 3, 3*x**2 - 4*x + 1, 2*x/9 + 25/9, -2079/4]
  
      """
      options.allowed_flags(args, ['auto', 'polys'])
  
      try:
          F, opt = poly_from_expr(f, *gens, **args)
      except PolificationFailed as exc:
          raise ComputationFailed('sturm', 1, exc)
  
      result = F.sturm(auto=opt.auto)
  
      if not opt.polys:
          return [r.as_expr() for r in result]
      else:
          return result
  
  
  @public
  def gff_list(f, *gens, **args):
      """
      Compute a list of greatest factorial factors of ``f``.
  
      Note that the input to ff() and rf() should be Poly instances to use the
      definitions here.
  
      Examples
      ========
  
      >>> from sympy import gff_list, ff, Poly
      >>> from sympy.abc import x
  
      >>> f = Poly(x**5 + 2*x**4 - x**3 - 2*x**2, x)
  
      >>> gff_list(f)
      [(Poly(x, x, domain='ZZ'), 1), (Poly(x + 2, x, domain='ZZ'), 4)]
  
      >>> (ff(Poly(x), 1)*ff(Poly(x + 2), 4)) == f
      True
  
      >>> f = Poly(x**12 + 6*x**11 - 11*x**10 - 56*x**9 + 220*x**8 + 208*x**7 - \
          1401*x**6 + 1090*x**5 + 2715*x**4 - 6720*x**3 - 1092*x**2 + 5040*x, x)
  
      >>> gff_list(f)
      [(Poly(x**3 + 7, x, domain='ZZ'), 2), (Poly(x**2 + 5*x, x, domain='ZZ'), 3)]
  
      >>> ff(Poly(x**3 + 7, x), 2)*ff(Poly(x**2 + 5*x, x), 3) == f
      True
  
      """
      options.allowed_flags(args, ['polys'])
  
      try:
          F, opt = poly_from_expr(f, *gens, **args)
      except PolificationFailed as exc:
          raise ComputationFailed('gff_list', 1, exc)
  
      factors = F.gff_list()
  
      if not opt.polys:
          return [(g.as_expr(), k) for g, k in factors]
      else:
          return factors
  
  
  @public
  def gff(f, *gens, **args):
      """Compute greatest factorial factorization of ``f``. """
      raise NotImplementedError('symbolic falling factorial')
  
  
  @public
  def sqf_norm(f, *gens, **args):
      """
      Compute square-free norm of ``f``.
  
      Returns ``s``, ``f``, ``r``, such that ``g(x) = f(x-sa)`` and
      ``r(x) = Norm(g(x))`` is a square-free polynomial over ``K``,
      where ``a`` is the algebraic extension of the ground domain.
  
      Examples
      ========
  
      >>> from sympy import sqf_norm, sqrt
      >>> from sympy.abc import x
  
      >>> sqf_norm(x**2 + 1, extension=[sqrt(3)])
      (1, x**2 - 2*sqrt(3)*x + 4, x**4 - 4*x**2 + 16)
  
      """
      options.allowed_flags(args, ['polys'])
  
      try:
          F, opt = poly_from_expr(f, *gens, **args)
      except PolificationFailed as exc:
          raise ComputationFailed('sqf_norm', 1, exc)
  
      s, g, r = F.sqf_norm()
  
      if not opt.polys:
          return Integer(s), g.as_expr(), r.as_expr()
      else:
          return Integer(s), g, r
  
  
  @public
  def sqf_part(f, *gens, **args):
      """
      Compute square-free part of ``f``.
  
      Examples
      ========
  
      >>> from sympy import sqf_part
      >>> from sympy.abc import x
  
      >>> sqf_part(x**3 - 3*x - 2)
      x**2 - x - 2
  
      """
      options.allowed_flags(args, ['polys'])
  
      try:
          F, opt = poly_from_expr(f, *gens, **args)
      except PolificationFailed as exc:
          raise ComputationFailed('sqf_part', 1, exc)
  
      result = F.sqf_part()
  
      if not opt.polys:
          return result.as_expr()
      else:
          return result
  
  
  def _sorted_factors(factors, method):
      """Sort a list of ``(expr, exp)`` pairs. """
      if method == 'sqf':
          def key(obj):
              poly, exp = obj
              rep = poly.rep.rep
              return (exp, len(rep), len(poly.gens), rep)
      else:
          def key(obj):
              poly, exp = obj
              rep = poly.rep.rep
              return (len(rep), len(poly.gens), exp, rep)
  
      return sorted(factors, key=key)
  
  
  def _factors_product(factors):
      """Multiply a list of ``(expr, exp)`` pairs. """
      return Mul(*[f.as_expr()**k for f, k in factors])
  
  
  def _symbolic_factor_list(expr, opt, method):
      """Helper function for :func:`_symbolic_factor`. """
      coeff, factors = S.One, []
  
      args = [i._eval_factor() if hasattr(i, '_eval_factor') else i
          for i in Mul.make_args(expr)]
      for arg in args:
          if arg.is_Number:
              coeff *= arg
              continue
          elif arg.is_Pow:
              base, exp = arg.args
              if base.is_Number and exp.is_Number:
                  coeff *= arg
                  continue
              if base.is_Number:
                  factors.append((base, exp))
                  continue
          else:
              base, exp = arg, S.One
  
          try:
              poly, _ = _poly_from_expr(base, opt)
          except PolificationFailed as exc:
              factors.append((exc.expr, exp))
          else:
              func = getattr(poly, method + '_list')
  
              _coeff, _factors = func()
              if _coeff is not S.One:
                  if exp.is_Integer:
                      coeff *= _coeff**exp
                  elif _coeff.is_positive:
                      factors.append((_coeff, exp))
                  else:
                      _factors.append((_coeff, S.One))
  
              if exp is S.One:
                  factors.extend(_factors)
              elif exp.is_integer:
                  factors.extend([(f, k*exp) for f, k in _factors])
              else:
                  other = []
  
                  for f, k in _factors:
                      if f.as_expr().is_positive:
                          factors.append((f, k*exp))
                      else:
                          other.append((f, k))
  
                  factors.append((_factors_product(other), exp))
      if method == 'sqf':
          factors = [(reduce(mul, (f for f, _ in factors if _ == k)), k)
                     for k in set(i for _, i in factors)]
  
      return coeff, factors
  
  
  def _symbolic_factor(expr, opt, method):
      """Helper function for :func:`_factor`. """
      if isinstance(expr, Expr):
          if hasattr(expr,'_eval_factor'):
              return expr._eval_factor()
          coeff, factors = _symbolic_factor_list(together(expr, fraction=opt['fraction']), opt, method)
          return _keep_coeff(coeff, _factors_product(factors))
      elif hasattr(expr, 'args'):
          return expr.func(*[_symbolic_factor(arg, opt, method) for arg in expr.args])
      elif hasattr(expr, '__iter__'):
          return expr.__class__([_symbolic_factor(arg, opt, method) for arg in expr])
      else:
          return expr
  
  
  def _generic_factor_list(expr, gens, args, method):
      """Helper function for :func:`sqf_list` and :func:`factor_list`. """
      options.allowed_flags(args, ['frac', 'polys'])
      opt = options.build_options(gens, args)
  
      expr = sympify(expr)
  
      if isinstance(expr, (Expr, Poly)):
          if isinstance(expr, Poly):
              numer, denom = expr, 1
          else:
              numer, denom = together(expr).as_numer_denom()
  
          cp, fp = _symbolic_factor_list(numer, opt, method)
          cq, fq = _symbolic_factor_list(denom, opt, method)
  
          if fq and not opt.frac:
              raise PolynomialError("a polynomial expected, got %s" % expr)
  
          _opt = opt.clone(dict(expand=True))
  
          for factors in (fp, fq):
              for i, (f, k) in enumerate(factors):
                  if not f.is_Poly:
                      f, _ = _poly_from_expr(f, _opt)
                      factors[i] = (f, k)
  
          fp = _sorted_factors(fp, method)
          fq = _sorted_factors(fq, method)
  
          if not opt.polys:
              fp = [(f.as_expr(), k) for f, k in fp]
              fq = [(f.as_expr(), k) for f, k in fq]
  
          coeff = cp/cq
  
          if not opt.frac:
              return coeff, fp
          else:
              return coeff, fp, fq
      else:
          raise PolynomialError("a polynomial expected, got %s" % expr)
  
  
  def _generic_factor(expr, gens, args, method):
      """Helper function for :func:`sqf` and :func:`factor`. """
      fraction = args.pop('fraction', True)
      options.allowed_flags(args, [])
      opt = options.build_options(gens, args)
      opt['fraction'] = fraction
      return _symbolic_factor(sympify(expr), opt, method)
  
  
  def to_rational_coeffs(f):
      """
      try to transform a polynomial to have rational coefficients
  
      try to find a transformation ``x = alpha*y``
  
      ``f(x) = lc*alpha**n * g(y)`` where ``g`` is a polynomial with
      rational coefficients, ``lc`` the leading coefficient.
  
      If this fails, try ``x = y + beta``
      ``f(x) = g(y)``
  
      Returns ``None`` if ``g`` not found;
      ``(lc, alpha, None, g)`` in case of rescaling
      ``(None, None, beta, g)`` in case of translation
  
      Notes
      =====
  
      Currently it transforms only polynomials without roots larger than 2.
  
      Examples
      ========
  
      >>> from sympy import sqrt, Poly, simplify
      >>> from sympy.polys.polytools import to_rational_coeffs
      >>> from sympy.abc import x
      >>> p = Poly(((x**2-1)*(x-2)).subs({x:x*(1 + sqrt(2))}), x, domain='EX')
      >>> lc, r, _, g = to_rational_coeffs(p)
      >>> lc, r
      (7 + 5*sqrt(2), 2 - 2*sqrt(2))
      >>> g
      Poly(x**3 + x**2 - 1/4*x - 1/4, x, domain='QQ')
      >>> r1 = simplify(1/r)
      >>> Poly(lc*r**3*(g.as_expr()).subs({x:x*r1}), x, domain='EX') == p
      True
  
      """
      from sympy.simplify.simplify import simplify
  
      def _try_rescale(f, f1=None):
          """
          try rescaling ``x -> alpha*x`` to convert f to a polynomial
          with rational coefficients.
          Returns ``alpha, f``; if the rescaling is successful,
          ``alpha`` is the rescaling factor, and ``f`` is the rescaled
          polynomial; else ``alpha`` is ``None``.
          """
          from sympy.core.add import Add
          if not len(f.gens) == 1 or not (f.gens[0]).is_Atom:
              return None, f
          n = f.degree()
          lc = f.LC()
          f1 = f1 or f1.monic()
          coeffs = f1.all_coeffs()[1:]
          coeffs = [simplify(coeffx) for coeffx in coeffs]
          if coeffs[-2]:
              rescale1_x = simplify(coeffs[-2]/coeffs[-1])
              coeffs1 = []
              for i in range(len(coeffs)):
                  coeffx = simplify(coeffs[i]*rescale1_x**(i + 1))
                  if not coeffx.is_rational:
                      break
                  coeffs1.append(coeffx)
              else:
                  rescale_x = simplify(1/rescale1_x)
                  x = f.gens[0]
                  v = [x**n]
                  for i in range(1, n + 1):
                      v.append(coeffs1[i - 1]*x**(n - i))
                  f = Add(*v)
                  f = Poly(f)
                  return lc, rescale_x, f
          return None
  
      def _try_translate(f, f1=None):
          """
          try translating ``x -> x + alpha`` to convert f to a polynomial
          with rational coefficients.
          Returns ``alpha, f``; if the translating is successful,
          ``alpha`` is the translating factor, and ``f`` is the shifted
          polynomial; else ``alpha`` is ``None``.
          """
          from sympy.core.add import Add
          if not len(f.gens) == 1 or not (f.gens[0]).is_Atom:
              return None, f
          n = f.degree()
          f1 = f1 or f1.monic()
          coeffs = f1.all_coeffs()[1:]
          c = simplify(coeffs[0])
          if c and not c.is_rational:
              func = Add
              if c.is_Add:
                  args = c.args
                  func = c.func
              else:
                  args = [c]
              c1, c2 = sift(args, lambda z: z.is_rational, binary=True)
              alpha = -func(*c2)/n
              f2 = f1.shift(alpha)
              return alpha, f2
          return None
  
      def _has_square_roots(p):
          """
          Return True if ``f`` is a sum with square roots but no other root
          """
          from sympy.core.exprtools import Factors
          coeffs = p.coeffs()
          has_sq = False
          for y in coeffs:
              for x in Add.make_args(y):
                  f = Factors(x).factors
                  r = [wx.q for b, wx in f.items() if
                      b.is_number and wx.is_Rational and wx.q >= 2]
                  if not r:
                      continue
                  if min(r) == 2:
                      has_sq = True
                  if max(r) > 2:
                      return False
          return has_sq
  
      if f.get_domain().is_EX and _has_square_roots(f):
          f1 = f.monic()
          r = _try_rescale(f, f1)
          if r:
              return r[0], r[1], None, r[2]
          else:
              r = _try_translate(f, f1)
              if r:
                  return None, None, r[0], r[1]
      return None
  
  
  def _torational_factor_list(p, x):
      """
      helper function to factor polynomial using to_rational_coeffs
  
      Examples
      ========
  
      >>> from sympy.polys.polytools import _torational_factor_list
      >>> from sympy.abc import x
      >>> from sympy import sqrt, expand, Mul
      >>> p = expand(((x**2-1)*(x-2)).subs({x:x*(1 + sqrt(2))}))
      >>> factors = _torational_factor_list(p, x); factors
      (-2, [(-x*(1 + sqrt(2))/2 + 1, 1), (-x*(1 + sqrt(2)) - 1, 1), (-x*(1 + sqrt(2)) + 1, 1)])
      >>> expand(factors[0]*Mul(*[z[0] for z in factors[1]])) == p
      True
      >>> p = expand(((x**2-1)*(x-2)).subs({x:x + sqrt(2)}))
      >>> factors = _torational_factor_list(p, x); factors
      (1, [(x - 2 + sqrt(2), 1), (x - 1 + sqrt(2), 1), (x + 1 + sqrt(2), 1)])
      >>> expand(factors[0]*Mul(*[z[0] for z in factors[1]])) == p
      True
  
      """
      from sympy.simplify.simplify import simplify
      p1 = Poly(p, x, domain='EX')
      n = p1.degree()
      res = to_rational_coeffs(p1)
      if not res:
          return None
      lc, r, t, g = res
      factors = factor_list(g.as_expr())
      if lc:
          c = simplify(factors[0]*lc*r**n)
          r1 = simplify(1/r)
          a = []
          for z in factors[1:][0]:
              a.append((simplify(z[0].subs({x: x*r1})), z[1]))
      else:
          c = factors[0]
          a = []
          for z in factors[1:][0]:
              a.append((z[0].subs({x: x - t}), z[1]))
      return (c, a)
  
  
  @public
  def sqf_list(f, *gens, **args):
      """
      Compute a list of square-free factors of ``f``.
  
      Examples
      ========
  
      >>> from sympy import sqf_list
      >>> from sympy.abc import x
  
      >>> sqf_list(2*x**5 + 16*x**4 + 50*x**3 + 76*x**2 + 56*x + 16)
      (2, [(x + 1, 2), (x + 2, 3)])
  
      """
      return _generic_factor_list(f, gens, args, method='sqf')
  
  
  @public
  def sqf(f, *gens, **args):
      """
      Compute square-free factorization of ``f``.
  
      Examples
      ========
  
      >>> from sympy import sqf
      >>> from sympy.abc import x
  
      >>> sqf(2*x**5 + 16*x**4 + 50*x**3 + 76*x**2 + 56*x + 16)
      2*(x + 1)**2*(x + 2)**3
  
      """
      return _generic_factor(f, gens, args, method='sqf')
  
  
  @public
  def factor_list(f, *gens, **args):
      """
      Compute a list of irreducible factors of ``f``.
  
      Examples
      ========
  
      >>> from sympy import factor_list
      >>> from sympy.abc import x, y
  
      >>> factor_list(2*x**5 + 2*x**4*y + 4*x**3 + 4*x**2*y + 2*x + 2*y)
      (2, [(x + y, 1), (x**2 + 1, 2)])
  
      """
      return _generic_factor_list(f, gens, args, method='factor')
  
  
  @public
  def factor(f, *gens, **args):
      """
      Compute the factorization of expression, ``f``, into irreducibles. (To
      factor an integer into primes, use ``factorint``.)
  
      There two modes implemented: symbolic and formal. If ``f`` is not an
      instance of :class:`Poly` and generators are not specified, then the
      former mode is used. Otherwise, the formal mode is used.
  
      In symbolic mode, :func:`factor` will traverse the expression tree and
      factor its components without any prior expansion, unless an instance
      of :class:`~.Add` is encountered (in this case formal factorization is
      used). This way :func:`factor` can handle large or symbolic exponents.
  
      By default, the factorization is computed over the rationals. To factor
      over other domain, e.g. an algebraic or finite field, use appropriate
      options: ``extension``, ``modulus`` or ``domain``.
  
      Examples
      ========
  
      >>> from sympy import factor, sqrt, exp
      >>> from sympy.abc import x, y
  
      >>> factor(2*x**5 + 2*x**4*y + 4*x**3 + 4*x**2*y + 2*x + 2*y)
      2*(x + y)*(x**2 + 1)**2
  
      >>> factor(x**2 + 1)
      x**2 + 1
      >>> factor(x**2 + 1, modulus=2)
      (x + 1)**2
      >>> factor(x**2 + 1, gaussian=True)
      (x - I)*(x + I)
  
      >>> factor(x**2 - 2, extension=sqrt(2))
      (x - sqrt(2))*(x + sqrt(2))
  
      >>> factor((x**2 - 1)/(x**2 + 4*x + 4))
      (x - 1)*(x + 1)/(x + 2)**2
      >>> factor((x**2 + 4*x + 4)**10000000*(x**2 + 1))
      (x + 2)**20000000*(x**2 + 1)
  
      By default, factor deals with an expression as a whole:
  
      >>> eq = 2**(x**2 + 2*x + 1)
      >>> factor(eq)
      2**(x**2 + 2*x + 1)
  
      If the ``deep`` flag is True then subexpressions will
      be factored:
  
      >>> factor(eq, deep=True)
      2**((x + 1)**2)
  
      If the ``fraction`` flag is False then rational expressions
      won't be combined. By default it is True.
  
      >>> factor(5*x + 3*exp(2 - 7*x), deep=True)
      (5*x*exp(7*x) + 3*exp(2))*exp(-7*x)
      >>> factor(5*x + 3*exp(2 - 7*x), deep=True, fraction=False)
      5*x + 3*exp(2)*exp(-7*x)
  
      See Also
      ========
      sympy.ntheory.factor_.factorint
  
      """
      f = sympify(f)
      if args.pop('deep', False):
          from sympy.simplify.simplify import bottom_up
          def _try_factor(expr):
              """
              Factor, but avoid changing the expression when unable to.
              """
              fac = factor(expr, *gens, **args)
              if fac.is_Mul or fac.is_Pow:
                  return fac
              return expr
  
          f = bottom_up(f, _try_factor)
          # clean up any subexpressions that may have been expanded
          # while factoring out a larger expression
          partials = {}
          muladd = f.atoms(Mul, Add)
          for p in muladd:
              fac = factor(p, *gens, **args)
              if (fac.is_Mul or fac.is_Pow) and fac != p:
                  partials[p] = fac
          return f.xreplace(partials)
  
      try:
          return _generic_factor(f, gens, args, method='factor')
      except PolynomialError as msg:
          if not f.is_commutative:
              from sympy.core.exprtools import factor_nc
              return factor_nc(f)
          else:
              raise PolynomialError(msg)
  
  
  @public
  def intervals(F, all=False, eps=None, inf=None, sup=None, strict=False, fast=False, sqf=False):
      """
      Compute isolating intervals for roots of ``f``.
  
      Examples
      ========
  
      >>> from sympy import intervals
      >>> from sympy.abc import x
  
      >>> intervals(x**2 - 3)
      [((-2, -1), 1), ((1, 2), 1)]
      >>> intervals(x**2 - 3, eps=1e-2)
      [((-26/15, -19/11), 1), ((19/11, 26/15), 1)]
  
      """
      if not hasattr(F, '__iter__'):
          try:
              F = Poly(F)
          except GeneratorsNeeded:
              return []
  
          return F.intervals(all=all, eps=eps, inf=inf, sup=sup, fast=fast, sqf=sqf)
      else:
          polys, opt = parallel_poly_from_expr(F, domain='QQ')
  
          if len(opt.gens) > 1:
              raise MultivariatePolynomialError
  
          for i, poly in enumerate(polys):
              polys[i] = poly.rep.rep
  
          if eps is not None:
              eps = opt.domain.convert(eps)
  
              if eps <= 0:
                  raise ValueError("'eps' must be a positive rational")
  
          if inf is not None:
              inf = opt.domain.convert(inf)
          if sup is not None:
              sup = opt.domain.convert(sup)
  
          intervals = dup_isolate_real_roots_list(polys, opt.domain,
              eps=eps, inf=inf, sup=sup, strict=strict, fast=fast)
  
          result = []
  
          for (s, t), indices in intervals:
              s, t = opt.domain.to_sympy(s), opt.domain.to_sympy(t)
              result.append(((s, t), indices))
  
          return result
  
  
  @public
  def refine_root(f, s, t, eps=None, steps=None, fast=False, check_sqf=False):
      """
      Refine an isolating interval of a root to the given precision.
  
      Examples
      ========
  
      >>> from sympy import refine_root
      >>> from sympy.abc import x
  
      >>> refine_root(x**2 - 3, 1, 2, eps=1e-2)
      (19/11, 26/15)
  
      """
      try:
          F = Poly(f)
          if not isinstance(f, Poly) and not F.gen.is_Symbol:
              # root of sin(x) + 1 is -1 but when someone
              # passes an Expr instead of Poly they may not expect
              # that the generator will be sin(x), not x
              raise PolynomialError("generator must be a Symbol")
      except GeneratorsNeeded:
          raise PolynomialError(
              "can't refine a root of %s, not a polynomial" % f)
  
      return F.refine_root(s, t, eps=eps, steps=steps, fast=fast, check_sqf=check_sqf)
  
  
  @public
  def count_roots(f, inf=None, sup=None):
      """
      Return the number of roots of ``f`` in ``[inf, sup]`` interval.
  
      If one of ``inf`` or ``sup`` is complex, it will return the number of roots
      in the complex rectangle with corners at ``inf`` and ``sup``.
  
      Examples
      ========
  
      >>> from sympy import count_roots, I
      >>> from sympy.abc import x
  
      >>> count_roots(x**4 - 4, -3, 3)
      2
      >>> count_roots(x**4 - 4, 0, 1 + 3*I)
      1
  
      """
      try:
          F = Poly(f, greedy=False)
          if not isinstance(f, Poly) and not F.gen.is_Symbol:
              # root of sin(x) + 1 is -1 but when someone
              # passes an Expr instead of Poly they may not expect
              # that the generator will be sin(x), not x
              raise PolynomialError("generator must be a Symbol")
      except GeneratorsNeeded:
          raise PolynomialError("can't count roots of %s, not a polynomial" % f)
  
      return F.count_roots(inf=inf, sup=sup)
  
  
  @public
  def real_roots(f, multiple=True):
      """
      Return a list of real roots with multiplicities of ``f``.
  
      Examples
      ========
  
      >>> from sympy import real_roots
      >>> from sympy.abc import x
  
      >>> real_roots(2*x**3 - 7*x**2 + 4*x + 4)
      [-1/2, 2, 2]
      """
      try:
          F = Poly(f, greedy=False)
          if not isinstance(f, Poly) and not F.gen.is_Symbol:
              # root of sin(x) + 1 is -1 but when someone
              # passes an Expr instead of Poly they may not expect
              # that the generator will be sin(x), not x
              raise PolynomialError("generator must be a Symbol")
      except GeneratorsNeeded:
          raise PolynomialError(
              "can't compute real roots of %s, not a polynomial" % f)
  
      return F.real_roots(multiple=multiple)
  
  
  @public
  def nroots(f, n=15, maxsteps=50, cleanup=True):
      """
      Compute numerical approximations of roots of ``f``.
  
      Examples
      ========
  
      >>> from sympy import nroots
      >>> from sympy.abc import x
  
      >>> nroots(x**2 - 3, n=15)
      [-1.73205080756888, 1.73205080756888]
      >>> nroots(x**2 - 3, n=30)
      [-1.73205080756887729352744634151, 1.73205080756887729352744634151]
  
      """
      try:
          F = Poly(f, greedy=False)
          if not isinstance(f, Poly) and not F.gen.is_Symbol:
              # root of sin(x) + 1 is -1 but when someone
              # passes an Expr instead of Poly they may not expect
              # that the generator will be sin(x), not x
              raise PolynomialError("generator must be a Symbol")
      except GeneratorsNeeded:
          raise PolynomialError(
              "can't compute numerical roots of %s, not a polynomial" % f)
  
      return F.nroots(n=n, maxsteps=maxsteps, cleanup=cleanup)
  
  
  @public
  def ground_roots(f, *gens, **args):
      """
      Compute roots of ``f`` by factorization in the ground domain.
  
      Examples
      ========
  
      >>> from sympy import ground_roots
      >>> from sympy.abc import x
  
      >>> ground_roots(x**6 - 4*x**4 + 4*x**3 - x**2)
      {0: 2, 1: 2}
  
      """
      options.allowed_flags(args, [])
  
      try:
          F, opt = poly_from_expr(f, *gens, **args)
          if not isinstance(f, Poly) and not F.gen.is_Symbol:
              # root of sin(x) + 1 is -1 but when someone
              # passes an Expr instead of Poly they may not expect
              # that the generator will be sin(x), not x
              raise PolynomialError("generator must be a Symbol")
      except PolificationFailed as exc:
          raise ComputationFailed('ground_roots', 1, exc)
  
      return F.ground_roots()
  
  
  @public
  def nth_power_roots_poly(f, n, *gens, **args):
      """
      Construct a polynomial with n-th powers of roots of ``f``.
  
      Examples
      ========
  
      >>> from sympy import nth_power_roots_poly, factor, roots
      >>> from sympy.abc import x
  
      >>> f = x**4 - x**2 + 1
      >>> g = factor(nth_power_roots_poly(f, 2))
  
      >>> g
      (x**2 - x + 1)**2
  
      >>> R_f = [ (r**2).expand() for r in roots(f) ]
      >>> R_g = roots(g).keys()
  
      >>> set(R_f) == set(R_g)
      True
  
      """
      options.allowed_flags(args, [])
  
      try:
          F, opt = poly_from_expr(f, *gens, **args)
          if not isinstance(f, Poly) and not F.gen.is_Symbol:
              # root of sin(x) + 1 is -1 but when someone
              # passes an Expr instead of Poly they may not expect
              # that the generator will be sin(x), not x
              raise PolynomialError("generator must be a Symbol")
      except PolificationFailed as exc:
          raise ComputationFailed('nth_power_roots_poly', 1, exc)
  
      result = F.nth_power_roots_poly(n)
  
      if not opt.polys:
          return result.as_expr()
      else:
          return result
  
  
  @public
  def cancel(f, *gens, **args):
      """
      Cancel common factors in a rational function ``f``.
  
      Examples
      ========
  
      >>> from sympy import cancel, sqrt, Symbol, together
      >>> from sympy.abc import x
      >>> A = Symbol('A', commutative=False)
  
      >>> cancel((2*x**2 - 2)/(x**2 - 2*x + 1))
      (2*x + 2)/(x - 1)
      >>> cancel((sqrt(3) + sqrt(15)*A)/(sqrt(2) + sqrt(10)*A))
      sqrt(6)/2
  
      Note: due to automatic distribution of Rationals, a sum divided by an integer
      will appear as a sum. To recover a rational form use `together` on the result:
  
      >>> cancel(x/2 + 1)
      x/2 + 1
      >>> together(_)
      (x + 2)/2
      """
      from sympy.core.exprtools import factor_terms
      from sympy.functions.elementary.piecewise import Piecewise
      options.allowed_flags(args, ['polys'])
  
      f = sympify(f)
  
      if not isinstance(f, (tuple, Tuple)):
          if f.is_Number or isinstance(f, Relational) or not isinstance(f, Expr):
              return f
          f = factor_terms(f, radical=True)
          p, q = f.as_numer_denom()
  
      elif len(f) == 2:
          p, q = f
      elif isinstance(f, Tuple):
          return factor_terms(f)
      else:
          raise ValueError('unexpected argument: %s' % f)
  
      try:
          (F, G), opt = parallel_poly_from_expr((p, q), *gens, **args)
      except PolificationFailed:
          if not isinstance(f, (tuple, Tuple)):
              return f.expand()
          else:
              return S.One, p, q
      except PolynomialError as msg:
          if f.is_commutative and not f.has(Piecewise):
              raise PolynomialError(msg)
          # Handling of noncommutative and/or piecewise expressions
          if f.is_Add or f.is_Mul:
              c, nc = sift(f.args, lambda x:
                  x.is_commutative is True and not x.has(Piecewise),
                  binary=True)
              nc = [cancel(i) for i in nc]
              return f.func(cancel(f.func(*c)), *nc)
          else:
              reps = []
              pot = preorder_traversal(f)
              next(pot)
              for e in pot:
                  # XXX: This should really skip anything that's not Expr.
                  if isinstance(e, (tuple, Tuple, BooleanAtom)):
                      continue
                  try:
                      reps.append((e, cancel(e)))
                      pot.skip()  # this was handled successfully
                  except NotImplementedError:
                      pass
              return f.xreplace(dict(reps))
  
      c, P, Q = F.cancel(G)
  
      if not isinstance(f, (tuple, Tuple)):
          return c*(P.as_expr()/Q.as_expr())
      else:
          if not opt.polys:
              return c, P.as_expr(), Q.as_expr()
          else:
              return c, P, Q
  
  
  @public
  def reduced(f, G, *gens, **args):
      """
      Reduces a polynomial ``f`` modulo a set of polynomials ``G``.
  
      Given a polynomial ``f`` and a set of polynomials ``G = (g_1, ..., g_n)``,
      computes a set of quotients ``q = (q_1, ..., q_n)`` and the remainder ``r``
      such that ``f = q_1*g_1 + ... + q_n*g_n + r``, where ``r`` vanishes or ``r``
      is a completely reduced polynomial with respect to ``G``.
  
      Examples
      ========
  
      >>> from sympy import reduced
      >>> from sympy.abc import x, y
  
      >>> reduced(2*x**4 + y**2 - x**2 + y**3, [x**3 - x, y**3 - y])
      ([2*x, 1], x**2 + y**2 + y)
  
      """
      options.allowed_flags(args, ['polys', 'auto'])
  
      try:
          polys, opt = parallel_poly_from_expr([f] + list(G), *gens, **args)
      except PolificationFailed as exc:
          raise ComputationFailed('reduced', 0, exc)
  
      domain = opt.domain
      retract = False
  
      if opt.auto and domain.is_Ring and not domain.is_Field:
          opt = opt.clone(dict(domain=domain.get_field()))
          retract = True
  
      from sympy.polys.rings import xring
      _ring, _ = xring(opt.gens, opt.domain, opt.order)
  
      for i, poly in enumerate(polys):
          poly = poly.set_domain(opt.domain).rep.to_dict()
          polys[i] = _ring.from_dict(poly)
  
      Q, r = polys[0].div(polys[1:])
  
      Q = [Poly._from_dict(dict(q), opt) for q in Q]
      r = Poly._from_dict(dict(r), opt)
  
      if retract:
          try:
              _Q, _r = [q.to_ring() for q in Q], r.to_ring()
          except CoercionFailed:
              pass
          else:
              Q, r = _Q, _r
  
      if not opt.polys:
          return [q.as_expr() for q in Q], r.as_expr()
      else:
          return Q, r
  
  
  @public
  def groebner(F, *gens, **args):
      """
      Computes the reduced Groebner basis for a set of polynomials.
  
      Use the ``order`` argument to set the monomial ordering that will be
      used to compute the basis. Allowed orders are ``lex``, ``grlex`` and
      ``grevlex``. If no order is specified, it defaults to ``lex``.
  
      For more information on Groebner bases, see the references and the docstring
      of :func:`~.solve_poly_system`.
  
      Examples
      ========
  
      Example taken from [1].
  
      >>> from sympy import groebner
      >>> from sympy.abc import x, y
  
      >>> F = [x*y - 2*y, 2*y**2 - x**2]
  
      >>> groebner(F, x, y, order='lex')
      GroebnerBasis([x**2 - 2*y**2, x*y - 2*y, y**3 - 2*y], x, y,
                    domain='ZZ', order='lex')
      >>> groebner(F, x, y, order='grlex')
      GroebnerBasis([y**3 - 2*y, x**2 - 2*y**2, x*y - 2*y], x, y,
                    domain='ZZ', order='grlex')
      >>> groebner(F, x, y, order='grevlex')
      GroebnerBasis([y**3 - 2*y, x**2 - 2*y**2, x*y - 2*y], x, y,
                    domain='ZZ', order='grevlex')
  
      By default, an improved implementation of the Buchberger algorithm is
      used. Optionally, an implementation of the F5B algorithm can be used. The
      algorithm can be set using the ``method`` flag or with the
      :func:`sympy.polys.polyconfig.setup` function.
  
      >>> F = [x**2 - x - 1, (2*x - 1) * y - (x**10 - (1 - x)**10)]
  
      >>> groebner(F, x, y, method='buchberger')
      GroebnerBasis([x**2 - x - 1, y - 55], x, y, domain='ZZ', order='lex')
      >>> groebner(F, x, y, method='f5b')
      GroebnerBasis([x**2 - x - 1, y - 55], x, y, domain='ZZ', order='lex')
  
      References
      ==========
  
      1. [Buchberger01]_
      2. [Cox97]_
  
      """
      return GroebnerBasis(F, *gens, **args)
  
  
  @public
  def is_zero_dimensional(F, *gens, **args):
      """
      Checks if the ideal generated by a Groebner basis is zero-dimensional.
  
      The algorithm checks if the set of monomials not divisible by the
      leading monomial of any element of ``F`` is bounded.
  
      References
      ==========
  
      David A. Cox, John B. Little, Donal O'Shea. Ideals, Varieties and
      Algorithms, 3rd edition, p. 230
  
      """
      return GroebnerBasis(F, *gens, **args).is_zero_dimensional
  
  
  @public
  class GroebnerBasis(Basic):
      """Represents a reduced Groebner basis. """
  
      def __new__(cls, F, *gens, **args):
          """Compute a reduced Groebner basis for a system of polynomials. """
          options.allowed_flags(args, ['polys', 'method'])
  
          try:
              polys, opt = parallel_poly_from_expr(F, *gens, **args)
          except PolificationFailed as exc:
              raise ComputationFailed('groebner', len(F), exc)
  
          from sympy.polys.rings import PolyRing
          ring = PolyRing(opt.gens, opt.domain, opt.order)
  
          polys = [ring.from_dict(poly.rep.to_dict()) for poly in polys if poly]
  
          G = _groebner(polys, ring, method=opt.method)
          G = [Poly._from_dict(g, opt) for g in G]
  
          return cls._new(G, opt)
  
      @classmethod
      def _new(cls, basis, options):
          obj = Basic.__new__(cls)
  
          obj._basis = tuple(basis)
          obj._options = options
  
          return obj
  
      @property
      def args(self):
          basis = (p.as_expr() for p in self._basis)
          return (Tuple(*basis), Tuple(*self._options.gens))
  
      @property
      def exprs(self):
          return [poly.as_expr() for poly in self._basis]
  
      @property
      def polys(self):
          return list(self._basis)
  
      @property
      def gens(self):
          return self._options.gens
  
      @property
      def domain(self):
          return self._options.domain
  
      @property
      def order(self):
          return self._options.order
  
      def __len__(self):
          return len(self._basis)
  
      def __iter__(self):
          if self._options.polys:
              return iter(self.polys)
          else:
              return iter(self.exprs)
  
      def __getitem__(self, item):
          if self._options.polys:
              basis = self.polys
          else:
              basis = self.exprs
  
          return basis[item]
  
      def __hash__(self):
          return hash((self._basis, tuple(self._options.items())))
  
      def __eq__(self, other):
          if isinstance(other, self.__class__):
              return self._basis == other._basis and self._options == other._options
          elif iterable(other):
              return self.polys == list(other) or self.exprs == list(other)
          else:
              return False
  
      def __ne__(self, other):
          return not self == other
  
      @property
      def is_zero_dimensional(self):
          """
          Checks if the ideal generated by a Groebner basis is zero-dimensional.
  
          The algorithm checks if the set of monomials not divisible by the
          leading monomial of any element of ``F`` is bounded.
  
          References
          ==========
  
          David A. Cox, John B. Little, Donal O'Shea. Ideals, Varieties and
          Algorithms, 3rd edition, p. 230
  
          """
          def single_var(monomial):
              return sum(map(bool, monomial)) == 1
  
          exponents = Monomial([0]*len(self.gens))
          order = self._options.order
  
          for poly in self.polys:
              monomial = poly.LM(order=order)
  
              if single_var(monomial):
                  exponents *= monomial
  
          # If any element of the exponents vector is zero, then there's
          # a variable for which there's no degree bound and the ideal
          # generated by this Groebner basis isn't zero-dimensional.
          return all(exponents)
  
      def fglm(self, order):
          """
          Convert a Groebner basis from one ordering to another.
  
          The FGLM algorithm converts reduced Groebner bases of zero-dimensional
          ideals from one ordering to another. This method is often used when it
          is infeasible to compute a Groebner basis with respect to a particular
          ordering directly.
  
          Examples
          ========
  
          >>> from sympy.abc import x, y
          >>> from sympy import groebner
  
          >>> F = [x**2 - 3*y - x + 1, y**2 - 2*x + y - 1]
          >>> G = groebner(F, x, y, order='grlex')
  
          >>> list(G.fglm('lex'))
          [2*x - y**2 - y + 1, y**4 + 2*y**3 - 3*y**2 - 16*y + 7]
          >>> list(groebner(F, x, y, order='lex'))
          [2*x - y**2 - y + 1, y**4 + 2*y**3 - 3*y**2 - 16*y + 7]
  
          References
          ==========
  
          .. [1] J.C. Faugere, P. Gianni, D. Lazard, T. Mora (1994). Efficient
                 Computation of Zero-dimensional Groebner Bases by Change of
                 Ordering
  
          """
          opt = self._options
  
          src_order = opt.order
          dst_order = monomial_key(order)
  
          if src_order == dst_order:
              return self
  
          if not self.is_zero_dimensional:
              raise NotImplementedError("can't convert Groebner bases of ideals with positive dimension")
  
          polys = list(self._basis)
          domain = opt.domain
  
          opt = opt.clone(dict(
              domain=domain.get_field(),
              order=dst_order,
          ))
  
          from sympy.polys.rings import xring
          _ring, _ = xring(opt.gens, opt.domain, src_order)
  
          for i, poly in enumerate(polys):
              poly = poly.set_domain(opt.domain).rep.to_dict()
              polys[i] = _ring.from_dict(poly)
  
          G = matrix_fglm(polys, _ring, dst_order)
          G = [Poly._from_dict(dict(g), opt) for g in G]
  
          if not domain.is_Field:
              G = [g.clear_denoms(convert=True)[1] for g in G]
              opt.domain = domain
  
          return self._new(G, opt)
  
      def reduce(self, expr, auto=True):
          """
          Reduces a polynomial modulo a Groebner basis.
  
          Given a polynomial ``f`` and a set of polynomials ``G = (g_1, ..., g_n)``,
          computes a set of quotients ``q = (q_1, ..., q_n)`` and the remainder ``r``
          such that ``f = q_1*f_1 + ... + q_n*f_n + r``, where ``r`` vanishes or ``r``
          is a completely reduced polynomial with respect to ``G``.
  
          Examples
          ========
  
          >>> from sympy import groebner, expand
          >>> from sympy.abc import x, y
  
          >>> f = 2*x**4 - x**2 + y**3 + y**2
          >>> G = groebner([x**3 - x, y**3 - y])
  
          >>> G.reduce(f)
          ([2*x, 1], x**2 + y**2 + y)
          >>> Q, r = _
  
          >>> expand(sum(q*g for q, g in zip(Q, G)) + r)
          2*x**4 - x**2 + y**3 + y**2
          >>> _ == f
          True
  
          """
          poly = Poly._from_expr(expr, self._options)
          polys = [poly] + list(self._basis)
  
          opt = self._options
          domain = opt.domain
  
          retract = False
  
          if auto and domain.is_Ring and not domain.is_Field:
              opt = opt.clone(dict(domain=domain.get_field()))
              retract = True
  
          from sympy.polys.rings import xring
          _ring, _ = xring(opt.gens, opt.domain, opt.order)
  
          for i, poly in enumerate(polys):
              poly = poly.set_domain(opt.domain).rep.to_dict()
              polys[i] = _ring.from_dict(poly)
  
          Q, r = polys[0].div(polys[1:])
  
          Q = [Poly._from_dict(dict(q), opt) for q in Q]
          r = Poly._from_dict(dict(r), opt)
  
          if retract:
              try:
                  _Q, _r = [q.to_ring() for q in Q], r.to_ring()
              except CoercionFailed:
                  pass
              else:
                  Q, r = _Q, _r
  
          if not opt.polys:
              return [q.as_expr() for q in Q], r.as_expr()
          else:
              return Q, r
  
      def contains(self, poly):
          """
          Check if ``poly`` belongs the ideal generated by ``self``.
  
          Examples
          ========
  
          >>> from sympy import groebner
          >>> from sympy.abc import x, y
  
          >>> f = 2*x**3 + y**3 + 3*y
          >>> G = groebner([x**2 + y**2 - 1, x*y - 2])
  
          >>> G.contains(f)
          True
          >>> G.contains(f + 1)
          False
  
          """
          return self.reduce(poly)[1] == 0
  
  
  @public
  def poly(expr, *gens, **args):
      """
      Efficiently transform an expression into a polynomial.
  
      Examples
      ========
  
      >>> from sympy import poly
      >>> from sympy.abc import x
  
      >>> poly(x*(x**2 + x - 1)**2)
      Poly(x**5 + 2*x**4 - x**3 - 2*x**2 + x, x, domain='ZZ')
  
      """
      options.allowed_flags(args, [])
  
      def _poly(expr, opt):
          terms, poly_terms = [], []
  
          for term in Add.make_args(expr):
              factors, poly_factors = [], []
  
              for factor in Mul.make_args(term):
                  if factor.is_Add:
                      poly_factors.append(_poly(factor, opt))
                  elif factor.is_Pow and factor.base.is_Add and \
                          factor.exp.is_Integer and factor.exp >= 0:
                      poly_factors.append(
                          _poly(factor.base, opt).pow(factor.exp))
                  else:
                      factors.append(factor)
  
              if not poly_factors:
                  terms.append(term)
              else:
                  product = poly_factors[0]
  
                  for factor in poly_factors[1:]:
                      product = product.mul(factor)
  
                  if factors:
                      factor = Mul(*factors)
  
                      if factor.is_Number:
                          product = product.mul(factor)
                      else:
                          product = product.mul(Poly._from_expr(factor, opt))
  
                  poly_terms.append(product)
  
          if not poly_terms:
              result = Poly._from_expr(expr, opt)
          else:
              result = poly_terms[0]
  
              for term in poly_terms[1:]:
                  result = result.add(term)
  
              if terms:
                  term = Add(*terms)
  
                  if term.is_Number:
                      result = result.add(term)
                  else:
                      result = result.add(Poly._from_expr(term, opt))
  
          return result.reorder(*opt.get('gens', ()), **args)
  
      expr = sympify(expr)
  
      if expr.is_Poly:
          return Poly(expr, *gens, **args)
  
      if 'expand' not in args:
          args['expand'] = False
  
      opt = options.build_options(gens, args)
  
      return _poly(expr, opt)
  