xhs-note-crawling/python/py/Lib/site-packages/sympy/polys/polytools.py

"""User-friendly public interface to polynomial functions. """

from __future__ import annotations

from functools import wraps, reduce
from operator import mul
from typing import Optional
from collections import Counter, defaultdict

from sympy.core import (
    S, Expr, Add, Tuple
)
from sympy.core.basic import Basic
from sympy.core.decorators import _sympifyit
from sympy.core.exprtools import Factors, factor_nc, factor_terms
from sympy.core.evalf import (
    pure_complex, evalf, fastlog, _evalf_with_bounded_error, quad_to_mpmath)
from sympy.core.function import Derivative
from sympy.core.mul import Mul, _keep_coeff
from sympy.core.intfunc import ilcm
from sympy.core.numbers import I, Integer, equal_valued
from sympy.core.relational import Relational, Equality
from sympy.core.sorting import ordered
from sympy.core.symbol import Dummy, Symbol
from sympy.core.sympify import sympify, _sympify
from sympy.core.traversal import preorder_traversal, bottom_up
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.domains.domainelement import DomainElement
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, DMF, ANP
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, public, filldedent
from sympy.utilities.exceptions import sympy_deprecation_warning
from sympy.utilities.iterables import iterable, sift

# Required to avoid errors
import sympy.polys

import mpmath
from mpmath.libmp.libhyper import NoConvergence



def _polifyit(func):
    @wraps(func)
    def wrapper(f, g):
        g = _sympify(g)
        if isinstance(g, Poly):
            return func(f, g)
        elif isinstance(g, Integer):
            g = f.from_expr(g, *f.gens, domain=f.domain)
            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:
                    sympy_deprecation_warning(
                        """
                        Mixing Poly with non-polynomial expressions in binary
                        operations is deprecated. Either explicitly convert
                        the non-Poly operand to a Poly with as_poly() or
                        convert the Poly to an Expr with as_expr().
                        """,
                        deprecated_since_version="1.6",
                        active_deprecations_target="deprecated-poly-nonpoly-binary-operations",
                    )
                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.

    See :ref:`polys-docs` for general documentation.

    Poly is a subclass of Basic rather than Expr but instances can be
    converted to Expr with the :py:meth:`~.Poly.as_expr` method.

    .. deprecated:: 1.6

       Combining Poly with non-Poly objects in binary operations is
       deprecated. Explicitly convert both objects to either Poly or Expr
       first. See :ref:`deprecated-poly-nonpoly-binary-operations`.

    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

    rep: DMP
    gens: tuple[Expr, ...]

    def __new__(cls, rep, *gens, **args) -> Poly:
        """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 isinstance(rep, (DMP, DMF, ANP, DomainElement)):
            return cls._from_domain_element(rep, opt)
        elif 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, evaluate=type(rep) is not str) # type: ignore

            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(
                "Cannot 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(
                "Cannot 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)

    @classmethod
    def _from_domain_element(cls, rep, opt):
        gens = opt.gens
        domain = opt.domain

        level = len(gens) - 1
        rep = [domain.convert(rep)]

        return cls.new(DMP.from_list(rep, level, domain), *gens)

    def __hash__(self):
        return super().__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 a :py:class:`~.Poly`

        Returns
        =======

        :py:class:`~.Domain`:
            Ground domain of the :py:class:`~.Poly`.

        Examples
        ========

        >>> from sympy import Poly, Symbol
        >>> x = Symbol('x')
        >>> p = Poly(x**2 + x)
        >>> p
        Poly(x**2 + x, x, domain='ZZ')
        >>> p.domain
        ZZ
        """
        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)

    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:
                g_coeff = f.rep.dom.from_sympy(g)
            except CoercionFailed:
                raise UnificationFailed("Cannot unify %s with %s" % (f, g))
            else:
                return f.rep.dom, f.per, f.rep, f.rep.ground_new(g_coeff)

        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.from_dict(dict(list(zip(f_monoms, f_coeffs))), lev, dom)
            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.from_dict(dict(list(zip(g_monoms, g_coeffs))), lev, dom)
            else:
                G = g.rep.convert(dom)
        else:
            raise UnificationFailed("Cannot 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 = [gen for j, gen in enumerate(f.gens) if j not in 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("Cannot 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.from_dict(rep, len(gens) - 1, f.rep.dom), 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("Cannot 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("Cannot 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("Cannot 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'):
            d = f.rep.degree(j)
            if d < 0:
                d = S.NegativeInfinity
            return d
        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 isinstance(spec, 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 isinstance(spec, 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("Cannot 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 - 2, x).revert(2)
        Traceback (most recent call last):
        ...
        NotReversible: only units are 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')

        See Also
        ========

        shift_list: Analogous method for multivariate polynomials.
        """
        return f.per(f.rep.shift(a))

    def shift_list(f, a):
        """
        Efficiently compute Taylor shift ``f(X + A)``.

        Examples
        ========

        >>> from sympy import Poly
        >>> from sympy.abc import x, y

        >>> Poly(x*y, [x,y]).shift_list([1, 2]) == Poly((x+1)*(y+2), [x,y])
        True

        See Also
        ========

        shift: Analogous method for univariate polynomials.
        """
        return f.per(f.rep.shift_list(a))

    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) -> tuple[Expr, list[tuple[Poly, int]]]:
        """
        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:
                if f.degree() == 0:
                    return f.as_expr(), []
                else:
                    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.

        See :func:`real_roots` for more explanation.

        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.

        See :func:`all_roots` for more explanation.

        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]

        """
        if f.is_multivariate:
            raise MultivariatePolynomialError(
                "Cannot 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()]
            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()]
            with mpmath.workdps(n):
                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

        from sympy.functions.elementary.complexes import sign
        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:
            try:
                # If roots did not converge try again with more extra precision.
                roots = mpmath.polyroots(coeffs, maxsteps=maxsteps,
                    cleanup=cleanup, error=False, extraprec=f.degree()*15)
                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(
                "Cannot 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 which_real_roots(f, candidates):
        """
        Find roots of a square-free polynomial ``f`` from ``candidates``.

        Explanation
        ===========

        If ``f`` is a square-free polynomial and ``candidates`` is a superset
        of the roots of ``f``, then ``f.which_real_roots(candidates)`` returns a
        list containing exactly the set of roots of ``f``. The domain must be
        :ref:`ZZ`, :ref:`QQ`, or :ref:`QQ(a)` and``f`` must be univariate and
        square-free.

        The list ``candidates`` must be a superset of the real roots of ``f``
        and ``f.which_real_roots(candidates)`` returns the set of real roots
        of ``f``. The output preserves the order of the order of ``candidates``.

        Examples
        ========

        >>> from sympy import Poly, sqrt
        >>> from sympy.abc import x

        >>> f = Poly(x**4 - 1)
        >>> f.which_real_roots([-1, 1, 0, -2, 2])
        [-1, 1]
        >>> f.which_real_roots([-1, 1, 1, 1, 1])
        [-1, 1]

        This method is useful as lifting to rational coefficients
        produced extraneous roots, which we can filter out with
        this method.

        >>> f = Poly(sqrt(2)*x**3 + x**2 - 1, x, extension=True)
        >>> f.lift()
        Poly(-2*x**6 + x**4 - 2*x**2 + 1, x, domain='QQ')
        >>> f.lift().real_roots()
        [-sqrt(2)/2, sqrt(2)/2]
        >>> f.which_real_roots(f.lift().real_roots())
        [sqrt(2)/2]

        This procedure is already done internally when calling
        `.real_roots()` on a polynomial with algebraic coefficients.

        >>> f.real_roots()
        [sqrt(2)/2]

        See Also
        ========

        same_root
        which_all_roots
        """
        if f.is_multivariate:
            raise MultivariatePolynomialError(
                "Must be a univariate polynomial")

        dom = f.get_domain()

        if not (dom.is_ZZ or dom.is_QQ or dom.is_AlgebraicField):
            raise NotImplementedError(
                "root counting not supported over %s" % dom)

        return f._which_roots(candidates, f.count_roots())

    def which_all_roots(f, candidates):
        """
        Find roots of a square-free polynomial ``f`` from ``candidates``.

        Explanation
        ===========

        If ``f`` is a square-free polynomial and ``candidates`` is a superset
        of the roots of ``f``, then ``f.which_all_roots(candidates)`` returns a
        list containing exactly the set of roots of ``f``. The polynomial``f``
        must be univariate and square-free.

        The list ``candidates`` must be a superset of the complex roots of
        ``f`` and ``f.which_all_roots(candidates)`` returns exactly the
        set of all complex roots of ``f``. The output preserves the order of
        the order of ``candidates``.

        Examples
        ========

        >>> from sympy import Poly, I
        >>> from sympy.abc import x

        >>> f = Poly(x**4 - 1)
        >>> f.which_all_roots([-1, 1, -I, I, 0])
        [-1, 1, -I, I]
        >>> f.which_all_roots([-1, 1, -I, I, I, I])
        [-1, 1, -I, I]

        This method is useful as lifting to rational coefficients
        produced extraneous roots, which we can filter out with
        this method.

        >>> f = Poly(x**2 + I*x - 1, x, extension=True)
        >>> f.lift()
        Poly(x**4 - x**2 + 1, x, domain='ZZ')
        >>> f.lift().all_roots()
        [CRootOf(x**4 - x**2 + 1, 0),
        CRootOf(x**4 - x**2 + 1, 1),
        CRootOf(x**4 - x**2 + 1, 2),
        CRootOf(x**4 - x**2 + 1, 3)]
        >>> f.which_all_roots(f.lift().all_roots())
        [CRootOf(x**4 - x**2 + 1, 0), CRootOf(x**4 - x**2 + 1, 2)]

        This procedure is already done internally when calling
        `.all_roots()` on a polynomial with algebraic coefficients,
        or polynomials with Gaussian domains.

        >>> f.all_roots()
        [CRootOf(x**4 - x**2 + 1, 0), CRootOf(x**4 - x**2 + 1, 2)]

        See Also
        ========

        same_root
        which_real_roots
        """
        if f.is_multivariate:
            raise MultivariatePolynomialError(
                "Must be a univariate polynomial")

        return f._which_roots(candidates, f.degree())

    def _which_roots(f, candidates, num_roots):
        prec = 10
        # using Counter bc its like an ordered set
        root_counts = Counter(candidates)

        while len(root_counts) > num_roots:
            for r in list(root_counts.keys()):
                # If f(r) != 0 then f(r).evalf() gives a float/complex with precision.
                f_r = f(r).evalf(prec, maxn=2*prec)
                if abs(f_r)._prec >= 2:
                    root_counts.pop(r)

            prec *= 2

        return list(root_counts.keys())

    def same_root(f, a, b):
        """
        Decide whether two roots of this polynomial are equal.

        Examples
        ========

        >>> from sympy import Poly, cyclotomic_poly, exp, I, pi
        >>> f = Poly(cyclotomic_poly(5))
        >>> r0 = exp(2*I*pi/5)
        >>> indices = [i for i, r in enumerate(f.all_roots()) if f.same_root(r, r0)]
        >>> print(indices)
        [3]

        Raises
        ======

        DomainError
            If the domain of the polynomial is not :ref:`ZZ`, :ref:`QQ`,
            :ref:`RR`, or :ref:`CC`.
        MultivariatePolynomialError
            If the polynomial is not univariate.
        PolynomialError
            If the polynomial is of degree < 2.

        See Also
        ========

        which_real_roots
        which_all_roots
        """
        if f.is_multivariate:
            raise MultivariatePolynomialError(
                "Must be a univariate polynomial")

        dom_delta_sq = f.rep.mignotte_sep_bound_squared()
        delta_sq = f.domain.get_field().to_sympy(dom_delta_sq)
        # We have delta_sq = delta**2, where delta is a lower bound on the
        # minimum separation between any two roots of this polynomial.
        # Let eps = delta/3, and define eps_sq = eps**2 = delta**2/9.
        eps_sq = delta_sq / 9

        r, _, _, _ = evalf(1/eps_sq, 1, {})
        n = fastlog(r)
        # Then 2^n > 1/eps**2.
        m = (n // 2) + (n % 2)
        # Then 2^(-m) < eps.
        ev = lambda x: quad_to_mpmath(_evalf_with_bounded_error(x, m=m))

        # Then for any complex numbers a, b we will have
        # |a - ev(a)| < eps and |b - ev(b)| < eps.
        # So if |ev(a) - ev(b)|**2 < eps**2, then
        # |ev(a) - ev(b)| < eps, hence |a - b| < 3*eps = delta.
        A, B = ev(a), ev(b)
        return (A.real - B.real)**2 + (A.imag - B.imag)**2 < eps_sq

    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))

    def make_monic_over_integers_by_scaling_roots(f):
        """
        Turn any univariate polynomial over :ref:`QQ` or :ref:`ZZ` into a monic
        polynomial over :ref:`ZZ`, by scaling the roots as necessary.

        Explanation
        ===========

        This operation can be performed whether or not *f* is irreducible; when
        it is, this can be understood as determining an algebraic integer
        generating the same field as a root of *f*.

        Examples
        ========

        >>> from sympy import Poly, S
        >>> from sympy.abc import x
        >>> f = Poly(x**2/2 + S(1)/4 * x + S(1)/8, x, domain='QQ')
        >>> f.make_monic_over_integers_by_scaling_roots()
        (Poly(x**2 + 2*x + 4, x, domain='ZZ'), 4)

        Returns
        =======

        Pair ``(g, c)``
            g is the polynomial

            c is the integer by which the roots had to be scaled

        """
        if not f.is_univariate or f.domain not in [ZZ, QQ]:
            raise ValueError('Polynomial must be univariate over ZZ or QQ.')
        if f.is_monic and f.domain == ZZ:
            return f, ZZ.one
        else:
            fm = f.monic()
            c, _ = fm.clear_denoms()
            return fm.transform(Poly(fm.gen), c).to_ring(), c

    def galois_group(f, by_name=False, max_tries=30, randomize=False):
        """
        Compute the Galois group of this polynomial.

        Examples
        ========

        >>> from sympy import Poly
        >>> from sympy.abc import x
        >>> f = Poly(x**4 - 2)
        >>> G, _ = f.galois_group(by_name=True)
        >>> print(G)
        S4TransitiveSubgroups.D4

        See Also
        ========

        sympy.polys.numberfields.galoisgroups.galois_group

        """
        from sympy.polys.numberfields.galoisgroups import (
            _galois_group_degree_3, _galois_group_degree_4_lookup,
            _galois_group_degree_5_lookup_ext_factor,
            _galois_group_degree_6_lookup,
        )
        if (not f.is_univariate
            or not f.is_irreducible
            or f.domain not in [ZZ, QQ]
        ):
            raise ValueError('Polynomial must be irreducible and univariate over ZZ or QQ.')
        gg = {
            3: _galois_group_degree_3,
            4: _galois_group_degree_4_lookup,
            5: _galois_group_degree_5_lookup_ext_factor,
            6: _galois_group_degree_6_lookup,
        }
        max_supported = max(gg.keys())
        n = f.degree()
        if n > max_supported:
            raise ValueError(f"Only polynomials up to degree {max_supported} are supported.")
        elif n < 1:
            raise ValueError("Constant polynomial has no Galois group.")
        elif n == 1:
            from sympy.combinatorics.galois import S1TransitiveSubgroups
            name, alt = S1TransitiveSubgroups.S1, True
        elif n == 2:
            from sympy.combinatorics.galois import S2TransitiveSubgroups
            name, alt = S2TransitiveSubgroups.S2, False
        else:
            g, _ = f.make_monic_over_integers_by_scaling_roots()
            name, alt = gg[n](g, max_tries=max_tries, randomize=randomize)
        G = name if by_name else name.get_perm_group()
        return G, alt

    @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 __truediv__(f, g):
        return f.as_expr()/g.as_expr()

    @_sympifyit('g', NotImplemented)
    def __rtruediv__(f, g):
        return g.as_expr()/f.as_expr()

    @_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 __bool__(f):
        return not f.is_zero

    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().__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("Cannot unify %s with %s" % (f, g))

        if len(f.gens) != len(g.gens):
            raise UnificationFailed("Cannot unify %s with %s" % (f, g))

        if not (isinstance(f.rep, DMP) and isinstance(g.rep, DMP)):
            raise UnificationFailed("Cannot 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. """
    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)

    from sympy.functions.elementary.piecewise import Piecewise
    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
    """

    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)))
    result = p.degree(gen)
    return Integer(result) if isinstance(result, int) else S.NegativeInfinity


@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

    >>> 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, 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:`sympy.core.intfunc.mod_inverse` function:

    >>> mod_inverse(3, 5)
    2
    >>> (S(2)/5).invert(S(7)/3)
    5/2

    See Also
    ========
    sympy.core.intfunc.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, includePRS=False, **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

    """
    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])
                # abs ensures that the gcd is always non-negative
                return abs(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:
                # abs ensures that the returned gcd is always non-negative
                return abs(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) -> Optional[Expr]:
        if not gens and not args:
            domain, numbers = construct_domain(seq)

            if not numbers:
                return domain.to_sympy(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

    """

    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 cannot 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 equal_valued(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()

    s_expr = [Integer(si) for si in s]

    if not opt.polys:
        return s_expr, g.as_expr(), r.as_expr()
    else:
        return s_expr, 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 _poly_sort_key(poly):
    """Sort a list of polys."""
    rep = poly.rep.to_list()
    return (len(rep), len(poly.gens), str(poly.domain), rep)


def _sorted_factors(factors, method):
    """Sort a list of ``(expr, exp)`` pairs. """
    if method == 'sqf':
        def key(obj):
            poly, exp = obj
            rep = poly.rep.to_list()
            return (exp, len(rep), len(poly.gens), str(poly.domain), rep)
    else:
        def key(obj):
            poly, exp = obj
            rep = poly.rep.to_list()
            return (len(rep), len(poly.gens), exp, str(poly.domain), 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 or (isinstance(arg, Expr) and pure_complex(arg)):
            coeff *= arg
            continue
        elif arg.is_Pow and arg.base != S.Exp1:
            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 {i for _, i in factors}]
    #collect duplicates
    rv = defaultdict(int)
    for k, v in factors:
        rv[k] += v
    return coeff, list(rv.items())


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({"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``.
        """
        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 len(coeffs) > 1 and 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``.
        """
        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.is_Add and not c.is_rational:
            rat, nonrat = sift(c.args,
                lambda z: z.is_rational is True, binary=True)
            alpha = -c.func(*nonrat)/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
        """
        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, deep=False, **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
    will not 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 deep:
        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:
        if not f.is_commutative:
            return factor_nc(f)
        else:
            raise


@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.to_list()

        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(
            "Cannot 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("Cannot count roots of %s, not a polynomial" % f)

    return F.count_roots(inf=inf, sup=sup)


@public
def all_roots(f, multiple=True, radicals=True, extension=False):
    """
    Returns the real and complex roots of ``f`` with multiplicities.

    Explanation
    ===========

    Finds all real and complex roots of a univariate polynomial with rational
    coefficients of any degree exactly. The roots are represented in the form
    given by :func:`~.rootof`. This is equivalent to using :func:`~.rootof` to
    find each of the indexed roots.

    Examples
    ========

    >>> from sympy import all_roots
    >>> from sympy.abc import x, y

    >>> print(all_roots(x**3 + 1))
    [-1, 1/2 - sqrt(3)*I/2, 1/2 + sqrt(3)*I/2]

    Simple radical formulae are used in some cases but the cubic and quartic
    formulae are avoided. Instead most non-rational roots will be represented
    as :class:`~.ComplexRootOf`:

    >>> print(all_roots(x**3 + x + 1))
    [CRootOf(x**3 + x + 1, 0), CRootOf(x**3 + x + 1, 1), CRootOf(x**3 + x + 1, 2)]

    All roots of any polynomial with rational coefficients of any degree can be
    represented using :py:class:`~.ComplexRootOf`. The use of
    :py:class:`~.ComplexRootOf` bypasses limitations on the availability of
    radical formulae for quintic and higher degree polynomials _[1]:

    >>> p = x**5 - x - 1
    >>> for r in all_roots(p): print(r)
    CRootOf(x**5 - x - 1, 0)
    CRootOf(x**5 - x - 1, 1)
    CRootOf(x**5 - x - 1, 2)
    CRootOf(x**5 - x - 1, 3)
    CRootOf(x**5 - x - 1, 4)
    >>> [r.evalf(3) for r in all_roots(p)]
    [1.17, -0.765 - 0.352*I, -0.765 + 0.352*I, 0.181 - 1.08*I, 0.181 + 1.08*I]

    Irrational algebraic coefficients are handled by :func:`all_roots`
    if `extension=True` is set.

    >>> from sympy import sqrt, expand
    >>> p = expand((x - sqrt(2))*(x - sqrt(3)))
    >>> print(p)
    x**2 - sqrt(3)*x - sqrt(2)*x + sqrt(6)
    >>> all_roots(p)
    Traceback (most recent call last):
    ...
    NotImplementedError: sorted roots not supported over EX
    >>> all_roots(p, extension=True)
    [sqrt(2), sqrt(3)]

    Algebraic coefficients can be complex as well.

    >>> from sympy import I
    >>> all_roots(x**2 - I, extension=True)
    [-sqrt(2)/2 - sqrt(2)*I/2, sqrt(2)/2 + sqrt(2)*I/2]
    >>> all_roots(x**2 - sqrt(2)*I, extension=True)
    [-2**(3/4)/2 - 2**(3/4)*I/2, 2**(3/4)/2 + 2**(3/4)*I/2]

    Transcendental coefficients cannot currently be handled by
    :func:`all_roots`. In the case of algebraic or transcendental coefficients
    :func:`~.ground_roots` might be able to find some roots by factorisation:

    >>> from sympy import ground_roots
    >>> ground_roots(p, x, extension=True)
    {sqrt(2): 1, sqrt(3): 1}

    If the coefficients are numeric then :func:`~.nroots` can be used to find
    all roots approximately:

    >>> from sympy import nroots
    >>> nroots(p, 5)
    [1.4142, 1.732]

    If the coefficients are symbolic then :func:`sympy.polys.polyroots.roots`
    or :func:`~.ground_roots` should be used instead:

    >>> from sympy import roots, ground_roots
    >>> p = x**2 - 3*x*y + 2*y**2
    >>> roots(p, x)
    {y: 1, 2*y: 1}
    >>> ground_roots(p, x)
    {y: 1, 2*y: 1}

    Parameters
    ==========

    f : :class:`~.Expr` or :class:`~.Poly`
        A univariate polynomial with rational (or ``Float``) coefficients.
    multiple : ``bool`` (default ``True``).
        Whether to return a ``list`` of roots or a list of root/multiplicity
        pairs.
    radicals : ``bool`` (default ``True``)
        Use simple radical formulae rather than :py:class:`~.ComplexRootOf` for
        some irrational roots.
    extension: ``bool`` (default ``False``)
        Whether to construct an algebraic extension domain before computing
        the roots. Setting to ``True`` is necessary for finding roots of a
        polynomial with (irrational) algebraic coefficients but can be slow.

    Returns
    =======

    A list of :class:`~.Expr` (usually :class:`~.ComplexRootOf`) representing
    the roots is returned with each root repeated according to its multiplicity
    as a root of ``f``. The roots are always uniquely ordered with real roots
    coming before complex roots. The real roots are in increasing order.
    Complex roots are ordered by increasing real part and then increasing
    imaginary part.

    If ``multiple=False`` is passed then a list of root/multiplicity pairs is
    returned instead.

    If ``radicals=False`` is passed then all roots will be represented as
    either rational numbers or :class:`~.ComplexRootOf`.

    See also
    ========

    Poly.all_roots:
        The underlying :class:`Poly` method used by :func:`~.all_roots`.
    rootof:
        Compute a single numbered root of a univariate polynomial.
    real_roots:
        Compute all the real roots using :func:`~.rootof`.
    ground_roots:
        Compute some roots in the ground domain by factorisation.
    nroots:
        Compute all roots using approximate numerical techniques.
    sympy.polys.polyroots.roots:
        Compute symbolic expressions for roots using radical formulae.

    References
    ==========

    .. [1] https://en.wikipedia.org/wiki/Abel%E2%80%93Ruffini_theorem
    """
    try:
        if isinstance(f, Poly):
            if extension and not f.domain.is_AlgebraicField:
                F = Poly(f.expr, extension=True)
            else:
                F = f
        else:
            if extension:
                F = Poly(f, extension=True)
            else:
                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(
            "Cannot compute real roots of %s, not a polynomial" % f)

    return F.all_roots(multiple=multiple, radicals=radicals)


@public
def real_roots(f, multiple=True, radicals=True, extension=False):
    """
    Returns the real roots of ``f`` with multiplicities.

    Explanation
    ===========

    Finds all real roots of a univariate polynomial with rational coefficients
    of any degree exactly. The roots are represented in the form given by
    :func:`~.rootof`. This is equivalent to using :func:`~.rootof` or
    :func:`~.all_roots` and filtering out only the real roots. However if only
    the real roots are needed then :func:`real_roots` is more efficient than
    :func:`~.all_roots` because it computes only the real roots and avoids
    costly complex root isolation routines.

    Examples
    ========

    >>> from sympy import real_roots
    >>> from sympy.abc import x, y

    >>> real_roots(2*x**3 - 7*x**2 + 4*x + 4)
    [-1/2, 2, 2]
    >>> real_roots(2*x**3 - 7*x**2 + 4*x + 4, multiple=False)
    [(-1/2, 1), (2, 2)]

    Real roots of any polynomial with rational coefficients of any degree can
    be represented using :py:class:`~.ComplexRootOf`:

    >>> p = x**9 + 2*x + 2
    >>> print(real_roots(p))
    [CRootOf(x**9 + 2*x + 2, 0)]
    >>> [r.evalf(3) for r in real_roots(p)]
    [-0.865]

    All rational roots will be returned as rational numbers. Roots of some
    simple factors will be expressed using radical or other formulae (unless
    ``radicals=False`` is passed). All other roots will be expressed as
    :class:`~.ComplexRootOf`.

    >>> p = (x + 7)*(x**2 - 2)*(x**3 + x + 1)
    >>> print(real_roots(p))
    [-7, -sqrt(2), CRootOf(x**3 + x + 1, 0), sqrt(2)]
    >>> print(real_roots(p, radicals=False))
    [-7, CRootOf(x**2 - 2, 0), CRootOf(x**3 + x + 1, 0), CRootOf(x**2 - 2, 1)]

    All returned root expressions will numerically evaluate to real numbers
    with no imaginary part. This is in contrast to the expressions generated by
    the cubic or quartic formulae as used by :func:`~.roots` which suffer from
    casus irreducibilis [1]_:

    >>> from sympy import roots
    >>> p = 2*x**3 - 9*x**2 - 6*x + 3
    >>> [r.evalf(5) for r in roots(p, multiple=True)]
    [5.0365 - 0.e-11*I, 0.33984 + 0.e-13*I, -0.87636 + 0.e-10*I]
    >>> [r.evalf(5) for r in real_roots(p, x)]
    [-0.87636, 0.33984, 5.0365]
    >>> [r.is_real for r in roots(p, multiple=True)]
    [None, None, None]
    >>> [r.is_real for r in real_roots(p)]
    [True, True, True]

    Using :func:`real_roots` is equivalent to using :func:`~.all_roots` (or
    :func:`~.rootof`) and filtering out only the real roots:

    >>> from sympy import all_roots
    >>> r = [r for r in all_roots(p) if r.is_real]
    >>> real_roots(p) == r
    True

    If only the real roots are wanted then using :func:`real_roots` is faster
    than using :func:`~.all_roots`. Using :func:`real_roots` avoids complex root
    isolation which can be a lot slower than real root isolation especially for
    polynomials of high degree which typically have many more complex roots
    than real roots.

    Irrational algebraic coefficients are handled by :func:`real_roots`
    if `extension=True` is set.

    >>> from sympy import sqrt, expand
    >>> p = expand((x - sqrt(2))*(x - sqrt(3)))
    >>> print(p)
    x**2 - sqrt(3)*x - sqrt(2)*x + sqrt(6)
    >>> real_roots(p)
    Traceback (most recent call last):
    ...
    NotImplementedError: sorted roots not supported over EX
    >>> real_roots(p, extension=True)
    [sqrt(2), sqrt(3)]

    Transcendental coefficients cannot currently be handled by
    :func:`real_roots`. In the case of algebraic or transcendental coefficients
    :func:`~.ground_roots` might be able to find some roots by factorisation:

    >>> from sympy import ground_roots
    >>> ground_roots(p, x, extension=True)
    {sqrt(2): 1, sqrt(3): 1}

    If the coefficients are numeric then :func:`~.nroots` can be used to find
    all roots approximately:

    >>> from sympy import nroots
    >>> nroots(p, 5)
    [1.4142, 1.732]

    If the coefficients are symbolic then :func:`sympy.polys.polyroots.roots`
    or :func:`~.ground_roots` should be used instead.

    >>> from sympy import roots, ground_roots
    >>> p = x**2 - 3*x*y + 2*y**2
    >>> roots(p, x)
    {y: 1, 2*y: 1}
    >>> ground_roots(p, x)
    {y: 1, 2*y: 1}

    Parameters
    ==========

    f : :class:`~.Expr` or :class:`~.Poly`
        A univariate polynomial with rational (or ``Float``) coefficients.
    multiple : ``bool`` (default ``True``).
        Whether to return a ``list`` of roots or a list of root/multiplicity
        pairs.
    radicals : ``bool`` (default ``True``)
        Use simple radical formulae rather than :py:class:`~.ComplexRootOf` for
        some irrational roots.
    extension: ``bool`` (default ``False``)
        Whether to construct an algebraic extension domain before computing
        the roots. Setting to ``True`` is necessary for finding roots of a
        polynomial with (irrational) algebraic coefficients but can be slow.

    Returns
    =======

    A list of :class:`~.Expr` (usually :class:`~.ComplexRootOf`) representing
    the real roots is returned. The roots are arranged in increasing order and
    are repeated according to their multiplicities as roots of ``f``.

    If ``multiple=False`` is passed then a list of root/multiplicity pairs is
    returned instead.

    If ``radicals=False`` is passed then all roots will be represented as
    either rational numbers or :class:`~.ComplexRootOf`.

    See also
    ========

    Poly.real_roots:
        The underlying :class:`Poly` method used by :func:`real_roots`.
    rootof:
        Compute a single numbered root of a univariate polynomial.
    all_roots:
        Compute all real and non-real roots using :func:`~.rootof`.
    ground_roots:
        Compute some roots in the ground domain by factorisation.
    nroots:
        Compute all roots using approximate numerical techniques.
    sympy.polys.polyroots.roots:
        Compute symbolic expressions for roots using radical formulae.

    References
    ==========

    .. [1] https://en.wikipedia.org/wiki/Casus_irreducibilis
    """
    try:
        if isinstance(f, Poly):
            if extension and not f.domain.is_AlgebraicField:
                F = Poly(f.expr, extension=True)
            else:
                F = f
        else:
            if extension:
                F = Poly(f, extension=True)
            else:
                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(
            "Cannot compute real roots of %s, not a polynomial" % f)

    return F.real_roots(multiple=multiple, radicals=radicals)


@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(
            "Cannot 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, _signsimp=True, **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.simplify.simplify import signsimp
    from sympy.polys.rings import sring
    options.allowed_flags(args, ['polys'])

    f = sympify(f)
    if _signsimp:
        f = signsimp(f)
    opt = {}
    if 'polys' in args:
        opt['polys'] = args['polys']

    if not isinstance(f, 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
        if isinstance(p, Poly) and isinstance(q, Poly):
            opt['gens'] = p.gens
            opt['domain'] = p.domain
            opt['polys'] = opt.get('polys', True)
        p, q = p.as_expr(), q.as_expr()
    else:
        raise ValueError('unexpected argument: %s' % f)

    from sympy.functions.elementary.piecewise import Piecewise
    try:
        if f.has(Piecewise):
            raise PolynomialError()
        R, (F, G) = sring((p, q), *gens, **args)
        if not R.ngens:
            if not isinstance(f, 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:
                if isinstance(e, BooleanAtom) or not isinstance(e, Expr):
                    continue
                try:
                    reps.append((e, cancel(e)))
                    pot.skip()  # this was handled successfully
                except NotImplementedError:
                    pass
            return f.xreplace(dict(reps))

    c, (P, Q) = 1, F.cancel(G)
    if opt.get('polys', False) and 'gens' not in opt:
        opt['gens'] = R.symbols

    if not isinstance(f, Tuple):
        return c*(P.as_expr()/Q.as_expr())
    else:
        P, Q = P.as_expr(), Q.as_expr()
        if not opt.get('polys', False):
            return c, P, Q
        else:
            return c, Poly(P, *gens, **opt), Poly(Q, *gens, **opt)


@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({"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("Cannot convert Groebner bases of ideals with positive dimension")

        polys = list(self._basis)
        domain = opt.domain

        opt = opt.clone({
            "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, Poly
        >>> 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

        # Using Poly input
        >>> f_poly = Poly(f, x, y)
        >>> G = groebner([Poly(x**3 - x), Poly(y**3 - y)])

        >>> G.reduce(f_poly)
        ([Poly(2*x, x, y, domain='ZZ'), Poly(1, x, y, domain='ZZ')], Poly(x**2 + y**2 + y, x, y, domain='ZZ'))

        """
        if isinstance(expr, Poly):

            if expr.gens != self._options.gens:
                raise ValueError("Polynomial generators don't match Groebner basis generators")
            poly = expr.set_domain(self._options.domain)
        else:

            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({"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 *= 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 += 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)


def named_poly(n, f, K, name, x, polys):
    r"""Common interface to the low-level polynomial generating functions
    in orthopolys and appellseqs.

    Parameters
    ==========

    n : int
        Index of the polynomial, which may or may not equal its degree.
    f : callable
        Low-level generating function to use.
    K : Domain or None
        Domain in which to perform the computations. If None, use the smallest
        field containing the rationals and the extra parameters of x (see below).
    name : str
        Name of an arbitrary individual polynomial in the sequence generated
        by f, only used in the error message for invalid n.
    x : seq
        The first element of this argument is the main variable of all
        polynomials in this sequence. Any further elements are extra
        parameters required by f.
    polys : bool, optional
        If True, return a Poly, otherwise (default) return an expression.
    """
    if n < 0:
        raise ValueError("Cannot generate %s of index %s" % (name, n))
    head, tail = x[0], x[1:]
    if K is None:
        K, tail = construct_domain(tail, field=True)
    poly = DMP(f(int(n), *tail, K), K)
    if head is None:
        poly = PurePoly.new(poly, Dummy('x'))
    else:
        poly = Poly.new(poly, head)
    return poly if polys else poly.as_expr()