Add files using upload-large-folder tool

lib/python3.12/site-packages/sympy/polys/domains/__init__.py

@@ -0,0 +1,57 @@
+"""Implementation of mathematical domains. """
+
+__all__ = [
+    'Domain', 'FiniteField', 'IntegerRing', 'RationalField', 'RealField',
+    'ComplexField', 'AlgebraicField', 'PolynomialRing', 'FractionField',
+    'ExpressionDomain', 'PythonRational',
+
+    'GF', 'FF', 'ZZ', 'QQ', 'ZZ_I', 'QQ_I', 'RR', 'CC', 'EX', 'EXRAW',
+]
+
+from .domain import Domain
+from .finitefield import FiniteField, FF, GF
+from .integerring import IntegerRing, ZZ
+from .rationalfield import RationalField, QQ
+from .algebraicfield import AlgebraicField
+from .gaussiandomains import ZZ_I, QQ_I
+from .realfield import RealField, RR
+from .complexfield import ComplexField, CC
+from .polynomialring import PolynomialRing
+from .fractionfield import FractionField
+from .expressiondomain import ExpressionDomain, EX
+from .expressionrawdomain import EXRAW
+from .pythonrational import PythonRational
+
+
+# This is imported purely for backwards compatibility because some parts of
+# the codebase used to import this from here and it's possible that downstream
+# does as well:
+from sympy.external.gmpy import GROUND_TYPES  # noqa: F401
+
+#
+# The rest of these are obsolete and provided only for backwards
+# compatibility:
+#
+
+from .pythonfinitefield import PythonFiniteField
+from .gmpyfinitefield import GMPYFiniteField
+from .pythonintegerring import PythonIntegerRing
+from .gmpyintegerring import GMPYIntegerRing
+from .pythonrationalfield import PythonRationalField
+from .gmpyrationalfield import GMPYRationalField
+
+FF_python = PythonFiniteField
+FF_gmpy = GMPYFiniteField
+
+ZZ_python = PythonIntegerRing
+ZZ_gmpy = GMPYIntegerRing
+
+QQ_python = PythonRationalField
+QQ_gmpy = GMPYRationalField
+
+__all__.extend((
+    'PythonFiniteField', 'GMPYFiniteField', 'PythonIntegerRing',
+    'GMPYIntegerRing', 'PythonRational', 'GMPYRationalField',
+
+    'FF_python', 'FF_gmpy', 'ZZ_python', 'ZZ_gmpy', 'QQ_python', 'QQ_gmpy',
+))

lib/python3.12/site-packages/sympy/polys/domains/algebraicfield.py

@@ -0,0 +1,607 @@
+"""Implementation of :class:`AlgebraicField` class. """
+
+
+from sympy.core.add import Add
+from sympy.core.mul import Mul
+from sympy.core.singleton import S
+from sympy.polys.domains.characteristiczero import CharacteristicZero
+from sympy.polys.domains.field import Field
+from sympy.polys.domains.simpledomain import SimpleDomain
+from sympy.polys.polyclasses import ANP
+from sympy.polys.polyerrors import CoercionFailed, DomainError, NotAlgebraic, IsomorphismFailed
+from sympy.utilities import public
+
+@public
+class AlgebraicField(Field, CharacteristicZero, SimpleDomain):
+    r"""Algebraic number field :ref:`QQ(a)`
+
+    A :ref:`QQ(a)` domain represents an `algebraic number field`_
+    `\mathbb{Q}(a)` as a :py:class:`~.Domain` in the domain system (see
+    :ref:`polys-domainsintro`).
+
+    A :py:class:`~.Poly` created from an expression involving `algebraic
+    numbers`_ will treat the algebraic numbers as generators if the generators
+    argument is not specified.
+
+    >>> from sympy import Poly, Symbol, sqrt
+    >>> x = Symbol('x')
+    >>> Poly(x**2 + sqrt(2))
+    Poly(x**2 + (sqrt(2)), x, sqrt(2), domain='ZZ')
+
+    That is a multivariate polynomial with ``sqrt(2)`` treated as one of the
+    generators (variables). If the generators are explicitly specified then
+    ``sqrt(2)`` will be considered to be a coefficient but by default the
+    :ref:`EX` domain is used. To make a :py:class:`~.Poly` with a :ref:`QQ(a)`
+    domain the argument ``extension=True`` can be given.
+
+    >>> Poly(x**2 + sqrt(2), x)
+    Poly(x**2 + sqrt(2), x, domain='EX')
+    >>> Poly(x**2 + sqrt(2), x, extension=True)
+    Poly(x**2 + sqrt(2), x, domain='QQ<sqrt(2)>')
+
+    A generator of the algebraic field extension can also be specified
+    explicitly which is particularly useful if the coefficients are all
+    rational but an extension field is needed (e.g. to factor the
+    polynomial).
+
+    >>> Poly(x**2 + 1)
+    Poly(x**2 + 1, x, domain='ZZ')
+    >>> Poly(x**2 + 1, extension=sqrt(2))
+    Poly(x**2 + 1, x, domain='QQ<sqrt(2)>')
+
+    It is possible to factorise a polynomial over a :ref:`QQ(a)` domain using
+    the ``extension`` argument to :py:func:`~.factor` or by specifying the domain
+    explicitly.
+
+    >>> from sympy import factor, QQ
+    >>> factor(x**2 - 2)
+    x**2 - 2
+    >>> factor(x**2 - 2, extension=sqrt(2))
+    (x - sqrt(2))*(x + sqrt(2))
+    >>> factor(x**2 - 2, domain='QQ<sqrt(2)>')
+    (x - sqrt(2))*(x + sqrt(2))
+    >>> factor(x**2 - 2, domain=QQ.algebraic_field(sqrt(2)))
+    (x - sqrt(2))*(x + sqrt(2))
+
+    The ``extension=True`` argument can be used but will only create an
+    extension that contains the coefficients which is usually not enough to
+    factorise the polynomial.
+
+    >>> p = x**3 + sqrt(2)*x**2 - 2*x - 2*sqrt(2)
+    >>> factor(p)                         # treats sqrt(2) as a symbol
+    (x + sqrt(2))*(x**2 - 2)
+    >>> factor(p, extension=True)
+    (x - sqrt(2))*(x + sqrt(2))**2
+    >>> factor(x**2 - 2, extension=True)  # all rational coefficients
+    x**2 - 2
+
+    It is also possible to use :ref:`QQ(a)` with the :py:func:`~.cancel`
+    and :py:func:`~.gcd` functions.
+
+    >>> from sympy import cancel, gcd
+    >>> cancel((x**2 - 2)/(x - sqrt(2)))
+    (x**2 - 2)/(x - sqrt(2))
+    >>> cancel((x**2 - 2)/(x - sqrt(2)), extension=sqrt(2))
+    x + sqrt(2)
+    >>> gcd(x**2 - 2, x - sqrt(2))
+    1
+    >>> gcd(x**2 - 2, x - sqrt(2), extension=sqrt(2))
+    x - sqrt(2)
+
+    When using the domain directly :ref:`QQ(a)` can be used as a constructor
+    to create instances which then support the operations ``+,-,*,**,/``. The
+    :py:meth:`~.Domain.algebraic_field` method is used to construct a
+    particular :ref:`QQ(a)` domain. The :py:meth:`~.Domain.from_sympy` method
+    can be used to create domain elements from normal SymPy expressions.
+
+    >>> K = QQ.algebraic_field(sqrt(2))
+    >>> K
+    QQ<sqrt(2)>
+    >>> xk = K.from_sympy(3 + 4*sqrt(2))
+    >>> xk  # doctest: +SKIP
+    ANP([4, 3], [1, 0, -2], QQ)
+
+    Elements of :ref:`QQ(a)` are instances of :py:class:`~.ANP` which have
+    limited printing support. The raw display shows the internal
+    representation of the element as the list ``[4, 3]`` representing the
+    coefficients of ``1`` and ``sqrt(2)`` for this element in the form
+    ``a * sqrt(2) + b * 1`` where ``a`` and ``b`` are elements of :ref:`QQ`.
+    The minimal polynomial for the generator ``(x**2 - 2)`` is also shown in
+    the :ref:`dup-representation` as the list ``[1, 0, -2]``. We can use
+    :py:meth:`~.Domain.to_sympy` to get a better printed form for the
+    elements and to see the results of operations.
+
+    >>> xk = K.from_sympy(3 + 4*sqrt(2))
+    >>> yk = K.from_sympy(2 + 3*sqrt(2))
+    >>> xk * yk  # doctest: +SKIP
+    ANP([17, 30], [1, 0, -2], QQ)
+    >>> K.to_sympy(xk * yk)
+    17*sqrt(2) + 30
+    >>> K.to_sympy(xk + yk)
+    5 + 7*sqrt(2)
+    >>> K.to_sympy(xk ** 2)
+    24*sqrt(2) + 41
+    >>> K.to_sympy(xk / yk)
+    sqrt(2)/14 + 9/7
+
+    Any expression representing an algebraic number can be used to generate
+    a :ref:`QQ(a)` domain provided its `minimal polynomial`_ can be computed.
+    The function :py:func:`~.minpoly` function is used for this.
+
+    >>> from sympy import exp, I, pi, minpoly
+    >>> g = exp(2*I*pi/3)
+    >>> g
+    exp(2*I*pi/3)
+    >>> g.is_algebraic
+    True
+    >>> minpoly(g, x)
+    x**2 + x + 1
+    >>> factor(x**3 - 1, extension=g)
+    (x - 1)*(x - exp(2*I*pi/3))*(x + 1 + exp(2*I*pi/3))
+
+    It is also possible to make an algebraic field from multiple extension
+    elements.
+
+    >>> K = QQ.algebraic_field(sqrt(2), sqrt(3))
+    >>> K
+    QQ<sqrt(2) + sqrt(3)>
+    >>> p = x**4 - 5*x**2 + 6
+    >>> factor(p)
+    (x**2 - 3)*(x**2 - 2)
+    >>> factor(p, domain=K)
+    (x - sqrt(2))*(x + sqrt(2))*(x - sqrt(3))*(x + sqrt(3))
+    >>> factor(p, extension=[sqrt(2), sqrt(3)])
+    (x - sqrt(2))*(x + sqrt(2))*(x - sqrt(3))*(x + sqrt(3))
+
+    Multiple extension elements are always combined together to make a single
+    `primitive element`_. In the case of ``[sqrt(2), sqrt(3)]`` the primitive
+    element chosen is ``sqrt(2) + sqrt(3)`` which is why the domain displays
+    as ``QQ<sqrt(2) + sqrt(3)>``. The minimal polynomial for the primitive
+    element is computed using the :py:func:`~.primitive_element` function.
+
+    >>> from sympy import primitive_element
+    >>> primitive_element([sqrt(2), sqrt(3)], x)
+    (x**4 - 10*x**2 + 1, [1, 1])
+    >>> minpoly(sqrt(2) + sqrt(3), x)
+    x**4 - 10*x**2 + 1
+
+    The extension elements that generate the domain can be accessed from the
+    domain using the :py:attr:`~.ext` and :py:attr:`~.orig_ext` attributes as
+    instances of :py:class:`~.AlgebraicNumber`. The minimal polynomial for
+    the primitive element as a :py:class:`~.DMP` instance is available as
+    :py:attr:`~.mod`.
+
+    >>> K = QQ.algebraic_field(sqrt(2), sqrt(3))
+    >>> K
+    QQ<sqrt(2) + sqrt(3)>
+    >>> K.ext
+    sqrt(2) + sqrt(3)
+    >>> K.orig_ext
+    (sqrt(2), sqrt(3))
+    >>> K.mod  # doctest: +SKIP
+    DMP_Python([1, 0, -10, 0, 1], QQ)
+
+    The `discriminant`_ of the field can be obtained from the
+    :py:meth:`~.discriminant` method, and an `integral basis`_ from the
+    :py:meth:`~.integral_basis` method. The latter returns a list of
+    :py:class:`~.ANP` instances by default, but can be made to return instances
+    of :py:class:`~.Expr` or :py:class:`~.AlgebraicNumber` by passing a ``fmt``
+    argument. The maximal order, or ring of integers, of the field can also be
+    obtained from the :py:meth:`~.maximal_order` method, as a
+    :py:class:`~sympy.polys.numberfields.modules.Submodule`.
+
+    >>> zeta5 = exp(2*I*pi/5)
+    >>> K = QQ.algebraic_field(zeta5)
+    >>> K
+    QQ<exp(2*I*pi/5)>
+    >>> K.discriminant()
+    125
+    >>> K = QQ.algebraic_field(sqrt(5))
+    >>> K
+    QQ<sqrt(5)>
+    >>> K.integral_basis(fmt='sympy')
+    [1, 1/2 + sqrt(5)/2]
+    >>> K.maximal_order()
+    Submodule[[2, 0], [1, 1]]/2
+
+    The factorization of a rational prime into prime ideals of the field is
+    computed by the :py:meth:`~.primes_above` method, which returns a list
+    of :py:class:`~sympy.polys.numberfields.primes.PrimeIdeal` instances.
+
+    >>> zeta7 = exp(2*I*pi/7)
+    >>> K = QQ.algebraic_field(zeta7)
+    >>> K
+    QQ<exp(2*I*pi/7)>
+    >>> K.primes_above(11)
+    [(11, _x**3 + 5*_x**2 + 4*_x - 1), (11, _x**3 - 4*_x**2 - 5*_x - 1)]
+
+    The Galois group of the Galois closure of the field can be computed (when
+    the minimal polynomial of the field is of sufficiently small degree).
+
+    >>> K.galois_group(by_name=True)[0]
+    S6TransitiveSubgroups.C6
+
+    Notes
+    =====
+
+    It is not currently possible to generate an algebraic extension over any
+    domain other than :ref:`QQ`. Ideally it would be possible to generate
+    extensions like ``QQ(x)(sqrt(x**2 - 2))``. This is equivalent to the
+    quotient ring ``QQ(x)[y]/(y**2 - x**2 + 2)`` and there are two
+    implementations of this kind of quotient ring/extension in the
+    :py:class:`~.QuotientRing` and :py:class:`~.MonogenicFiniteExtension`
+    classes.  Each of those implementations needs some work to make them fully
+    usable though.
+
+    .. _algebraic number field: https://en.wikipedia.org/wiki/Algebraic_number_field
+    .. _algebraic numbers: https://en.wikipedia.org/wiki/Algebraic_number
+    .. _discriminant: https://en.wikipedia.org/wiki/Discriminant_of_an_algebraic_number_field
+    .. _integral basis: https://en.wikipedia.org/wiki/Algebraic_number_field#Integral_basis
+    .. _minimal polynomial: https://en.wikipedia.org/wiki/Minimal_polynomial_(field_theory)
+    .. _primitive element: https://en.wikipedia.org/wiki/Primitive_element_theorem
+    """
+
+    dtype = ANP
+
+    is_AlgebraicField = is_Algebraic = True
+    is_Numerical = True
+
+    has_assoc_Ring = False
+    has_assoc_Field = True
+
+    def __init__(self, dom, *ext, alias=None):
+        r"""
+        Parameters
+        ==========
+
+        dom : :py:class:`~.Domain`
+            The base field over which this is an extension field.
+            Currently only :ref:`QQ` is accepted.
+
+        *ext : One or more :py:class:`~.Expr`
+            Generators of the extension. These should be expressions that are
+            algebraic over `\mathbb{Q}`.
+
+        alias : str, :py:class:`~.Symbol`, None, optional (default=None)
+            If provided, this will be used as the alias symbol for the
+            primitive element of the :py:class:`~.AlgebraicField`.
+            If ``None``, while ``ext`` consists of exactly one
+            :py:class:`~.AlgebraicNumber`, its alias (if any) will be used.
+        """
+        if not dom.is_QQ:
+            raise DomainError("ground domain must be a rational field")
+
+        from sympy.polys.numberfields import to_number_field
+        if len(ext) == 1 and isinstance(ext[0], tuple):
+            orig_ext = ext[0][1:]
+        else:
+            orig_ext = ext
+
+        if alias is None and len(ext) == 1:
+            alias = getattr(ext[0], 'alias', None)
+
+        self.orig_ext = orig_ext
+        """
+        Original elements given to generate the extension.
+
+        >>> from sympy import QQ, sqrt
+        >>> K = QQ.algebraic_field(sqrt(2), sqrt(3))
+        >>> K.orig_ext
+        (sqrt(2), sqrt(3))
+        """
+
+        self.ext = to_number_field(ext, alias=alias)
+        """
+        Primitive element used for the extension.
+
+        >>> from sympy import QQ, sqrt
+        >>> K = QQ.algebraic_field(sqrt(2), sqrt(3))
+        >>> K.ext
+        sqrt(2) + sqrt(3)
+        """
+
+        self.mod = self.ext.minpoly.rep
+        """
+        Minimal polynomial for the primitive element of the extension.
+
+        >>> from sympy import QQ, sqrt
+        >>> K = QQ.algebraic_field(sqrt(2))
+        >>> K.mod
+        DMP([1, 0, -2], QQ)
+        """
+
+        self.domain = self.dom = dom
+
+        self.ngens = 1
+        self.symbols = self.gens = (self.ext,)
+        self.unit = self([dom(1), dom(0)])
+
+        self.zero = self.dtype.zero(self.mod.to_list(), dom)
+        self.one = self.dtype.one(self.mod.to_list(), dom)
+
+        self._maximal_order = None
+        self._discriminant = None
+        self._nilradicals_mod_p = {}
+
+    def new(self, element):
+        return self.dtype(element, self.mod.to_list(), self.dom)
+
+    def __str__(self):
+        return str(self.dom) + '<' + str(self.ext) + '>'
+
+    def __hash__(self):
+        return hash((self.__class__.__name__, self.dtype, self.dom, self.ext))
+
+    def __eq__(self, other):
+        """Returns ``True`` if two domains are equivalent. """
+        if isinstance(other, AlgebraicField):
+            return self.dtype == other.dtype and self.ext == other.ext
+        else:
+            return NotImplemented
+
+    def algebraic_field(self, *extension, alias=None):
+        r"""Returns an algebraic field, i.e. `\mathbb{Q}(\alpha, \ldots)`. """
+        return AlgebraicField(self.dom, *((self.ext,) + extension), alias=alias)
+
+    def to_alg_num(self, a):
+        """Convert ``a`` of ``dtype`` to an :py:class:`~.AlgebraicNumber`. """
+        return self.ext.field_element(a)
+
+    def to_sympy(self, a):
+        """Convert ``a`` of ``dtype`` to a SymPy object. """
+        # Precompute a converter to be reused:
+        if not hasattr(self, '_converter'):
+            self._converter = _make_converter(self)
+
+        return self._converter(a)
+
+    def from_sympy(self, a):
+        """Convert SymPy's expression to ``dtype``. """
+        try:
+            return self([self.dom.from_sympy(a)])
+        except CoercionFailed:
+            pass
+
+        from sympy.polys.numberfields import to_number_field
+
+        try:
+            return self(to_number_field(a, self.ext).native_coeffs())
+        except (NotAlgebraic, IsomorphismFailed):
+            raise CoercionFailed(
+                "%s is not a valid algebraic number in %s" % (a, self))
+
+    def from_ZZ(K1, a, K0):
+        """Convert a Python ``int`` object to ``dtype``. """
+        return K1(K1.dom.convert(a, K0))
+
+    def from_ZZ_python(K1, a, K0):
+        """Convert a Python ``int`` object to ``dtype``. """
+        return K1(K1.dom.convert(a, K0))
+
+    def from_QQ(K1, a, K0):
+        """Convert a Python ``Fraction`` object to ``dtype``. """
+        return K1(K1.dom.convert(a, K0))
+
+    def from_QQ_python(K1, a, K0):
+        """Convert a Python ``Fraction`` object to ``dtype``. """
+        return K1(K1.dom.convert(a, K0))
+
+    def from_ZZ_gmpy(K1, a, K0):
+        """Convert a GMPY ``mpz`` object to ``dtype``. """
+        return K1(K1.dom.convert(a, K0))
+
+    def from_QQ_gmpy(K1, a, K0):
+        """Convert a GMPY ``mpq`` object to ``dtype``. """
+        return K1(K1.dom.convert(a, K0))
+
+    def from_RealField(K1, a, K0):
+        """Convert a mpmath ``mpf`` object to ``dtype``. """
+        return K1(K1.dom.convert(a, K0))
+
+    def get_ring(self):
+        """Returns a ring associated with ``self``. """
+        raise DomainError('there is no ring associated with %s' % self)
+
+    def is_positive(self, a):
+        """Returns True if ``a`` is positive. """
+        return self.dom.is_positive(a.LC())
+
+    def is_negative(self, a):
+        """Returns True if ``a`` is negative. """
+        return self.dom.is_negative(a.LC())
+
+    def is_nonpositive(self, a):
+        """Returns True if ``a`` is non-positive. """
+        return self.dom.is_nonpositive(a.LC())
+
+    def is_nonnegative(self, a):
+        """Returns True if ``a`` is non-negative. """
+        return self.dom.is_nonnegative(a.LC())
+
+    def numer(self, a):
+        """Returns numerator of ``a``. """
+        return a
+
+    def denom(self, a):
+        """Returns denominator of ``a``. """
+        return self.one
+
+    def from_AlgebraicField(K1, a, K0):
+        """Convert AlgebraicField element 'a' to another AlgebraicField """
+        return K1.from_sympy(K0.to_sympy(a))
+
+    def from_GaussianIntegerRing(K1, a, K0):
+        """Convert a GaussianInteger element 'a' to ``dtype``. """
+        return K1.from_sympy(K0.to_sympy(a))
+
+    def from_GaussianRationalField(K1, a, K0):
+        """Convert a GaussianRational element 'a' to ``dtype``. """
+        return K1.from_sympy(K0.to_sympy(a))
+
+    def _do_round_two(self):
+        from sympy.polys.numberfields.basis import round_two
+        ZK, dK = round_two(self, radicals=self._nilradicals_mod_p)
+        self._maximal_order = ZK
+        self._discriminant = dK
+
+    def maximal_order(self):
+        """
+        Compute the maximal order, or ring of integers, of the field.
+
+        Returns
+        =======
+
+        :py:class:`~sympy.polys.numberfields.modules.Submodule`.
+
+        See Also
+        ========
+
+        integral_basis
+
+        """
+        if self._maximal_order is None:
+            self._do_round_two()
+        return self._maximal_order
+
+    def integral_basis(self, fmt=None):
+        r"""
+        Get an integral basis for the field.
+
+        Parameters
+        ==========
+
+        fmt : str, None, optional (default=None)
+            If ``None``, return a list of :py:class:`~.ANP` instances.
+            If ``"sympy"``, convert each element of the list to an
+            :py:class:`~.Expr`, using ``self.to_sympy()``.
+            If ``"alg"``, convert each element of the list to an
+            :py:class:`~.AlgebraicNumber`, using ``self.to_alg_num()``.
+
+        Examples
+        ========
+
+        >>> from sympy import QQ, AlgebraicNumber, sqrt
+        >>> alpha = AlgebraicNumber(sqrt(5), alias='alpha')
+        >>> k = QQ.algebraic_field(alpha)
+        >>> B0 = k.integral_basis()
+        >>> B1 = k.integral_basis(fmt='sympy')
+        >>> B2 = k.integral_basis(fmt='alg')
+        >>> print(B0[1])  # doctest: +SKIP
+        ANP([mpq(1,2), mpq(1,2)], [mpq(1,1), mpq(0,1), mpq(-5,1)], QQ)
+        >>> print(B1[1])
+        1/2 + alpha/2
+        >>> print(B2[1])
+        alpha/2 + 1/2
+
+        In the last two cases we get legible expressions, which print somewhat
+        differently because of the different types involved:
+
+        >>> print(type(B1[1]))
+        <class 'sympy.core.add.Add'>
+        >>> print(type(B2[1]))
+        <class 'sympy.core.numbers.AlgebraicNumber'>
+
+        See Also
+        ========
+
+        to_sympy
+        to_alg_num
+        maximal_order
+        """
+        ZK = self.maximal_order()
+        M = ZK.QQ_matrix
+        n = M.shape[1]
+        B = [self.new(list(reversed(M[:, j].flat()))) for j in range(n)]
+        if fmt == 'sympy':
+            return [self.to_sympy(b) for b in B]
+        elif fmt == 'alg':
+            return [self.to_alg_num(b) for b in B]
+        return B
+
+    def discriminant(self):
+        """Get the discriminant of the field."""
+        if self._discriminant is None:
+            self._do_round_two()
+        return self._discriminant
+
+    def primes_above(self, p):
+        """Compute the prime ideals lying above a given rational prime *p*."""
+        from sympy.polys.numberfields.primes import prime_decomp
+        ZK = self.maximal_order()
+        dK = self.discriminant()
+        rad = self._nilradicals_mod_p.get(p)
+        return prime_decomp(p, ZK=ZK, dK=dK, radical=rad)
+
+    def galois_group(self, by_name=False, max_tries=30, randomize=False):
+        """
+        Compute the Galois group of the Galois closure of this field.
+
+        Examples
+        ========
+
+        If the field is Galois, the order of the group will equal the degree
+        of the field:
+
+        >>> from sympy import QQ
+        >>> from sympy.abc import x
+        >>> k = QQ.alg_field_from_poly(x**4 + 1)
+        >>> G, _ = k.galois_group()
+        >>> G.order()
+        4
+
+        If the field is not Galois, then its Galois closure is a proper
+        extension, and the order of the Galois group will be greater than the
+        degree of the field:
+
+        >>> k = QQ.alg_field_from_poly(x**4 - 2)
+        >>> G, _ = k.galois_group()
+        >>> G.order()
+        8
+
+        See Also
+        ========
+
+        sympy.polys.numberfields.galoisgroups.galois_group
+
+        """
+        return self.ext.minpoly_of_element().galois_group(
+            by_name=by_name, max_tries=max_tries, randomize=randomize)
+
+
+def _make_converter(K):
+    """Construct the converter to convert back to Expr"""
+    # Precompute the effect of converting to SymPy and expanding expressions
+    # like (sqrt(2) + sqrt(3))**2. Asking Expr to do the expansion on every
+    # conversion from K to Expr is slow. Here we compute the expansions for
+    # each power of the generator and collect together the resulting algebraic
+    # terms and the rational coefficients into a matrix.
+
+    gen = K.ext.as_expr()
+    todom = K.dom.from_sympy
+
+    # We'll let Expr compute the expansions. We won't make any presumptions
+    # about what this results in except that it is QQ-linear in some terms
+    # that we will call algebraics. The final result will be expressed in
+    # terms of those.
+    powers = [S.One, gen]
+    for n in range(2, K.mod.degree()):
+        powers.append((gen * powers[-1]).expand())
+
+    # Collect the rational coefficients and algebraic Expr that can
+    # map the ANP coefficients into an expanded SymPy expression
+    terms = [dict(t.as_coeff_Mul()[::-1] for t in Add.make_args(p)) for p in powers]
+    algebraics = set().union(*terms)
+    matrix = [[todom(t.get(a, S.Zero)) for t in terms] for a in algebraics]
+
+    # Create a function to do the conversion efficiently:
+
+    def converter(a):
+        """Convert a to Expr using converter"""
+        ai = a.to_list()[::-1]
+        tosympy = K.dom.to_sympy
+        coeffs_dom = [sum(mij*aj for mij, aj in zip(mi, ai)) for mi in matrix]
+        coeffs_sympy = [tosympy(c) for c in coeffs_dom]
+        res = Add(*(Mul(c, a) for c, a in zip(coeffs_sympy, algebraics)))
+        return res
+
+    return converter

lib/python3.12/site-packages/sympy/polys/domains/characteristiczero.py

@@ -0,0 +1,15 @@
+"""Implementation of :class:`CharacteristicZero` class. """
+
+
+from sympy.polys.domains.domain import Domain
+from sympy.utilities import public
+
+@public
+class CharacteristicZero(Domain):
+    """Domain that has infinite number of elements. """
+
+    has_CharacteristicZero = True
+
+    def characteristic(self):
+        """Return the characteristic of this domain. """
+        return 0

lib/python3.12/site-packages/sympy/polys/domains/complexfield.py

@@ -0,0 +1,185 @@
+"""Implementation of :class:`ComplexField` class. """
+
+
+from sympy.external.gmpy import SYMPY_INTS
+from sympy.core.numbers import Float, I
+from sympy.polys.domains.characteristiczero import CharacteristicZero
+from sympy.polys.domains.field import Field
+from sympy.polys.domains.gaussiandomains import QQ_I
+from sympy.polys.domains.mpelements import MPContext
+from sympy.polys.domains.simpledomain import SimpleDomain
+from sympy.polys.polyerrors import DomainError, CoercionFailed
+from sympy.utilities import public
+
+@public
+class ComplexField(Field, CharacteristicZero, SimpleDomain):
+    """Complex numbers up to the given precision. """
+
+    rep = 'CC'
+
+    is_ComplexField = is_CC = True
+
+    is_Exact = False
+    is_Numerical = True
+
+    has_assoc_Ring = False
+    has_assoc_Field = True
+
+    _default_precision = 53
+
+    @property
+    def has_default_precision(self):
+        return self.precision == self._default_precision
+
+    @property
+    def precision(self):
+        return self._context.prec
+
+    @property
+    def dps(self):
+        return self._context.dps
+
+    @property
+    def tolerance(self):
+        return self._context.tolerance
+
+    def __init__(self, prec=_default_precision, dps=None, tol=None):
+        context = MPContext(prec, dps, tol, False)
+        context._parent = self
+        self._context = context
+
+        self._dtype = context.mpc
+        self.zero = self.dtype(0)
+        self.one = self.dtype(1)
+
+    @property
+    def tp(self):
+        # XXX: Domain treats tp as an alis of dtype. Here we need to two
+        # separate things: dtype is a callable to make/convert instances.
+        # We use tp with isinstance to check if an object is an instance
+        # of the domain already.
+        return self._dtype
+
+    def dtype(self, x, y=0):
+        # XXX: This is needed because mpmath does not recognise fmpz.
+        # It might be better to add conversion routines to mpmath and if that
+        # happens then this can be removed.
+        if isinstance(x, SYMPY_INTS):
+            x = int(x)
+        if isinstance(y, SYMPY_INTS):
+            y = int(y)
+        return self._dtype(x, y)
+
+    def __eq__(self, other):
+        return (isinstance(other, ComplexField)
+           and self.precision == other.precision
+           and self.tolerance == other.tolerance)
+
+    def __hash__(self):
+        return hash((self.__class__.__name__, self._dtype, self.precision, self.tolerance))
+
+    def to_sympy(self, element):
+        """Convert ``element`` to SymPy number. """
+        return Float(element.real, self.dps) + I*Float(element.imag, self.dps)
+
+    def from_sympy(self, expr):
+        """Convert SymPy's number to ``dtype``. """
+        number = expr.evalf(n=self.dps)
+        real, imag = number.as_real_imag()
+
+        if real.is_Number and imag.is_Number:
+            return self.dtype(real, imag)
+        else:
+            raise CoercionFailed("expected complex number, got %s" % expr)
+
+    def from_ZZ(self, element, base):
+        return self.dtype(element)
+
+    def from_ZZ_gmpy(self, element, base):
+        return self.dtype(int(element))
+
+    def from_ZZ_python(self, element, base):
+        return self.dtype(element)
+
+    def from_QQ(self, element, base):
+        return self.dtype(int(element.numerator)) / int(element.denominator)
+
+    def from_QQ_python(self, element, base):
+        return self.dtype(element.numerator) / element.denominator
+
+    def from_QQ_gmpy(self, element, base):
+        return self.dtype(int(element.numerator)) / int(element.denominator)
+
+    def from_GaussianIntegerRing(self, element, base):
+        return self.dtype(int(element.x), int(element.y))
+
+    def from_GaussianRationalField(self, element, base):
+        x = element.x
+        y = element.y
+        return (self.dtype(int(x.numerator)) / int(x.denominator) +
+                self.dtype(0, int(y.numerator)) / int(y.denominator))
+
+    def from_AlgebraicField(self, element, base):
+        return self.from_sympy(base.to_sympy(element).evalf(self.dps))
+
+    def from_RealField(self, element, base):
+        return self.dtype(element)
+
+    def from_ComplexField(self, element, base):
+        if self == base:
+            return element
+        else:
+            return self.dtype(element)
+
+    def get_ring(self):
+        """Returns a ring associated with ``self``. """
+        raise DomainError("there is no ring associated with %s" % self)
+
+    def get_exact(self):
+        """Returns an exact domain associated with ``self``. """
+        return QQ_I
+
+    def is_negative(self, element):
+        """Returns ``False`` for any ``ComplexElement``. """
+        return False
+
+    def is_positive(self, element):
+        """Returns ``False`` for any ``ComplexElement``. """
+        return False
+
+    def is_nonnegative(self, element):
+        """Returns ``False`` for any ``ComplexElement``. """
+        return False
+
+    def is_nonpositive(self, element):
+        """Returns ``False`` for any ``ComplexElement``. """
+        return False
+
+    def gcd(self, a, b):
+        """Returns GCD of ``a`` and ``b``. """
+        return self.one
+
+    def lcm(self, a, b):
+        """Returns LCM of ``a`` and ``b``. """
+        return a*b
+
+    def almosteq(self, a, b, tolerance=None):
+        """Check if ``a`` and ``b`` are almost equal. """
+        return self._context.almosteq(a, b, tolerance)
+
+    def is_square(self, a):
+        """Returns ``True``. Every complex number has a complex square root."""
+        return True
+
+    def exsqrt(self, a):
+        r"""Returns the principal complex square root of ``a``.
+
+        Explanation
+        ===========
+        The argument of the principal square root is always within
+        $(-\frac{\pi}{2}, \frac{\pi}{2}]$. The square root may be
+        slightly inaccurate due to floating point rounding error.
+        """
+        return a ** 0.5
+
+CC = ComplexField()

lib/python3.12/site-packages/sympy/polys/domains/compositedomain.py

@@ -0,0 +1,52 @@
+"""Implementation of :class:`CompositeDomain` class. """
+
+
+from sympy.polys.domains.domain import Domain
+from sympy.polys.polyerrors import GeneratorsError
+
+from sympy.utilities import public
+
+@public
+class CompositeDomain(Domain):
+    """Base class for composite domains, e.g. ZZ[x], ZZ(X). """
+
+    is_Composite = True
+
+    gens, ngens, symbols, domain = [None]*4
+
+    def inject(self, *symbols):
+        """Inject generators into this domain.  """
+        if not (set(self.symbols) & set(symbols)):
+            return self.__class__(self.domain, self.symbols + symbols, self.order)
+        else:
+            raise GeneratorsError("common generators in %s and %s" % (self.symbols, symbols))
+
+    def drop(self, *symbols):
+        """Drop generators from this domain. """
+        symset = set(symbols)
+        newsyms = tuple(s for s in self.symbols if s not in symset)
+        domain = self.domain.drop(*symbols)
+        if not newsyms:
+            return domain
+        else:
+            return self.__class__(domain, newsyms, self.order)
+
+    def set_domain(self, domain):
+        """Set the ground domain of this domain. """
+        return self.__class__(domain, self.symbols, self.order)
+
+    @property
+    def is_Exact(self):
+        """Returns ``True`` if this domain is exact. """
+        return self.domain.is_Exact
+
+    def get_exact(self):
+        """Returns an exact version of this domain. """
+        return self.set_domain(self.domain.get_exact())
+
+    @property
+    def has_CharacteristicZero(self):
+        return self.domain.has_CharacteristicZero
+
+    def characteristic(self):
+        return self.domain.characteristic()

lib/python3.12/site-packages/sympy/polys/domains/domain.py

@@ -0,0 +1,1372 @@
+"""Implementation of :class:`Domain` class. """
+
+from __future__ import annotations
+from typing import Any
+
+from sympy.core.numbers import AlgebraicNumber
+from sympy.core import Basic, sympify
+from sympy.core.sorting import ordered
+from sympy.external.gmpy import GROUND_TYPES
+from sympy.polys.domains.domainelement import DomainElement
+from sympy.polys.orderings import lex
+from sympy.polys.polyerrors import UnificationFailed, CoercionFailed, DomainError
+from sympy.polys.polyutils import _unify_gens, _not_a_coeff
+from sympy.utilities import public
+from sympy.utilities.iterables import is_sequence
+
+
+@public
+class Domain:
+    """Superclass for all domains in the polys domains system.
+
+    See :ref:`polys-domainsintro` for an introductory explanation of the
+    domains system.
+
+    The :py:class:`~.Domain` class is an abstract base class for all of the
+    concrete domain types. There are many different :py:class:`~.Domain`
+    subclasses each of which has an associated ``dtype`` which is a class
+    representing the elements of the domain. The coefficients of a
+    :py:class:`~.Poly` are elements of a domain which must be a subclass of
+    :py:class:`~.Domain`.
+
+    Examples
+    ========
+
+    The most common example domains are the integers :ref:`ZZ` and the
+    rationals :ref:`QQ`.
+
+    >>> from sympy import Poly, symbols, Domain
+    >>> x, y = symbols('x, y')
+    >>> p = Poly(x**2 + y)
+    >>> p
+    Poly(x**2 + y, x, y, domain='ZZ')
+    >>> p.domain
+    ZZ
+    >>> isinstance(p.domain, Domain)
+    True
+    >>> Poly(x**2 + y/2)
+    Poly(x**2 + 1/2*y, x, y, domain='QQ')
+
+    The domains can be used directly in which case the domain object e.g.
+    (:ref:`ZZ` or :ref:`QQ`) can be used as a constructor for elements of
+    ``dtype``.
+
+    >>> from sympy import ZZ, QQ
+    >>> ZZ(2)
+    2
+    >>> ZZ.dtype  # doctest: +SKIP
+    <class 'int'>
+    >>> type(ZZ(2))  # doctest: +SKIP
+    <class 'int'>
+    >>> QQ(1, 2)
+    1/2
+    >>> type(QQ(1, 2))  # doctest: +SKIP
+    <class 'sympy.polys.domains.pythonrational.PythonRational'>
+
+    The corresponding domain elements can be used with the arithmetic
+    operations ``+,-,*,**`` and depending on the domain some combination of
+    ``/,//,%`` might be usable. For example in :ref:`ZZ` both ``//`` (floor
+    division) and ``%`` (modulo division) can be used but ``/`` (true
+    division) cannot. Since :ref:`QQ` is a :py:class:`~.Field` its elements
+    can be used with ``/`` but ``//`` and ``%`` should not be used. Some
+    domains have a :py:meth:`~.Domain.gcd` method.
+
+    >>> ZZ(2) + ZZ(3)
+    5
+    >>> ZZ(5) // ZZ(2)
+    2
+    >>> ZZ(5) % ZZ(2)
+    1
+    >>> QQ(1, 2) / QQ(2, 3)
+    3/4
+    >>> ZZ.gcd(ZZ(4), ZZ(2))
+    2
+    >>> QQ.gcd(QQ(2,7), QQ(5,3))
+    1/21
+    >>> ZZ.is_Field
+    False
+    >>> QQ.is_Field
+    True
+
+    There are also many other domains including:
+
+        1. :ref:`GF(p)` for finite fields of prime order.
+        2. :ref:`RR` for real (floating point) numbers.
+        3. :ref:`CC` for complex (floating point) numbers.
+        4. :ref:`QQ(a)` for algebraic number fields.
+        5. :ref:`K[x]` for polynomial rings.
+        6. :ref:`K(x)` for rational function fields.
+        7. :ref:`EX` for arbitrary expressions.
+
+    Each domain is represented by a domain object and also an implementation
+    class (``dtype``) for the elements of the domain. For example the
+    :ref:`K[x]` domains are represented by a domain object which is an
+    instance of :py:class:`~.PolynomialRing` and the elements are always
+    instances of :py:class:`~.PolyElement`. The implementation class
+    represents particular types of mathematical expressions in a way that is
+    more efficient than a normal SymPy expression which is of type
+    :py:class:`~.Expr`. The domain methods :py:meth:`~.Domain.from_sympy` and
+    :py:meth:`~.Domain.to_sympy` are used to convert from :py:class:`~.Expr`
+    to a domain element and vice versa.
+
+    >>> from sympy import Symbol, ZZ, Expr
+    >>> x = Symbol('x')
+    >>> K = ZZ[x]           # polynomial ring domain
+    >>> K
+    ZZ[x]
+    >>> type(K)             # class of the domain
+    <class 'sympy.polys.domains.polynomialring.PolynomialRing'>
+    >>> K.dtype             # class of the elements
+    <class 'sympy.polys.rings.PolyElement'>
+    >>> p_expr = x**2 + 1   # Expr
+    >>> p_expr
+    x**2 + 1
+    >>> type(p_expr)
+    <class 'sympy.core.add.Add'>
+    >>> isinstance(p_expr, Expr)
+    True
+    >>> p_domain = K.from_sympy(p_expr)
+    >>> p_domain            # domain element
+    x**2 + 1
+    >>> type(p_domain)
+    <class 'sympy.polys.rings.PolyElement'>
+    >>> K.to_sympy(p_domain) == p_expr
+    True
+
+    The :py:meth:`~.Domain.convert_from` method is used to convert domain
+    elements from one domain to another.
+
+    >>> from sympy import ZZ, QQ
+    >>> ez = ZZ(2)
+    >>> eq = QQ.convert_from(ez, ZZ)
+    >>> type(ez)  # doctest: +SKIP
+    <class 'int'>
+    >>> type(eq)  # doctest: +SKIP
+    <class 'sympy.polys.domains.pythonrational.PythonRational'>
+
+    Elements from different domains should not be mixed in arithmetic or other
+    operations: they should be converted to a common domain first.  The domain
+    method :py:meth:`~.Domain.unify` is used to find a domain that can
+    represent all the elements of two given domains.
+
+    >>> from sympy import ZZ, QQ, symbols
+    >>> x, y = symbols('x, y')
+    >>> ZZ.unify(QQ)
+    QQ
+    >>> ZZ[x].unify(QQ)
+    QQ[x]
+    >>> ZZ[x].unify(QQ[y])
+    QQ[x,y]
+
+    If a domain is a :py:class:`~.Ring` then is might have an associated
+    :py:class:`~.Field` and vice versa. The :py:meth:`~.Domain.get_field` and
+    :py:meth:`~.Domain.get_ring` methods will find or create the associated
+    domain.
+
+    >>> from sympy import ZZ, QQ, Symbol
+    >>> x = Symbol('x')
+    >>> ZZ.has_assoc_Field
+    True
+    >>> ZZ.get_field()
+    QQ
+    >>> QQ.has_assoc_Ring
+    True
+    >>> QQ.get_ring()
+    ZZ
+    >>> K = QQ[x]
+    >>> K
+    QQ[x]
+    >>> K.get_field()
+    QQ(x)
+
+    See also
+    ========
+
+    DomainElement: abstract base class for domain elements
+    construct_domain: construct a minimal domain for some expressions
+
+    """
+
+    dtype: type | None = None
+    """The type (class) of the elements of this :py:class:`~.Domain`:
+
+    >>> from sympy import ZZ, QQ, Symbol
+    >>> ZZ.dtype
+    <class 'int'>
+    >>> z = ZZ(2)
+    >>> z
+    2
+    >>> type(z)
+    <class 'int'>
+    >>> type(z) == ZZ.dtype
+    True
+
+    Every domain has an associated **dtype** ("datatype") which is the
+    class of the associated domain elements.
+
+    See also
+    ========
+
+    of_type
+    """
+
+    zero: Any = None
+    """The zero element of the :py:class:`~.Domain`:
+
+    >>> from sympy import QQ
+    >>> QQ.zero
+    0
+    >>> QQ.of_type(QQ.zero)
+    True
+
+    See also
+    ========
+
+    of_type
+    one
+    """
+
+    one: Any = None
+    """The one element of the :py:class:`~.Domain`:
+
+    >>> from sympy import QQ
+    >>> QQ.one
+    1
+    >>> QQ.of_type(QQ.one)
+    True
+
+    See also
+    ========
+
+    of_type
+    zero
+    """
+
+    is_Ring = False
+    """Boolean flag indicating if the domain is a :py:class:`~.Ring`.
+
+    >>> from sympy import ZZ
+    >>> ZZ.is_Ring
+    True
+
+    Basically every :py:class:`~.Domain` represents a ring so this flag is
+    not that useful.
+
+    See also
+    ========
+
+    is_PID
+    is_Field
+    get_ring
+    has_assoc_Ring
+    """
+
+    is_Field = False
+    """Boolean flag indicating if the domain is a :py:class:`~.Field`.
+
+    >>> from sympy import ZZ, QQ
+    >>> ZZ.is_Field
+    False
+    >>> QQ.is_Field
+    True
+
+    See also
+    ========
+
+    is_PID
+    is_Ring
+    get_field
+    has_assoc_Field
+    """
+
+    has_assoc_Ring = False
+    """Boolean flag indicating if the domain has an associated
+    :py:class:`~.Ring`.
+
+    >>> from sympy import QQ
+    >>> QQ.has_assoc_Ring
+    True
+    >>> QQ.get_ring()
+    ZZ
+
+    See also
+    ========
+
+    is_Field
+    get_ring
+    """
+
+    has_assoc_Field = False
+    """Boolean flag indicating if the domain has an associated
+    :py:class:`~.Field`.
+
+    >>> from sympy import ZZ
+    >>> ZZ.has_assoc_Field
+    True
+    >>> ZZ.get_field()
+    QQ
+
+    See also
+    ========
+
+    is_Field
+    get_field
+    """
+
+    is_FiniteField = is_FF = False
+    is_IntegerRing = is_ZZ = False
+    is_RationalField = is_QQ = False
+    is_GaussianRing = is_ZZ_I = False
+    is_GaussianField = is_QQ_I = False
+    is_RealField = is_RR = False
+    is_ComplexField = is_CC = False
+    is_AlgebraicField = is_Algebraic = False
+    is_PolynomialRing = is_Poly = False
+    is_FractionField = is_Frac = False
+    is_SymbolicDomain = is_EX = False
+    is_SymbolicRawDomain = is_EXRAW = False
+    is_FiniteExtension = False
+
+    is_Exact = True
+    is_Numerical = False
+
+    is_Simple = False
+    is_Composite = False
+
+    is_PID = False
+    """Boolean flag indicating if the domain is a `principal ideal domain`_.
+
+    >>> from sympy import ZZ
+    >>> ZZ.has_assoc_Field
+    True
+    >>> ZZ.get_field()
+    QQ
+
+    .. _principal ideal domain: https://en.wikipedia.org/wiki/Principal_ideal_domain
+
+    See also
+    ========
+
+    is_Field
+    get_field
+    """
+
+    has_CharacteristicZero = False
+
+    rep: str | None = None
+    alias: str | None = None
+
+    def __init__(self):
+        raise NotImplementedError
+
+    def __str__(self):
+        return self.rep
+
+    def __repr__(self):
+        return str(self)
+
+    def __hash__(self):
+        return hash((self.__class__.__name__, self.dtype))
+
+    def new(self, *args):
+        return self.dtype(*args)
+
+    @property
+    def tp(self):
+        """Alias for :py:attr:`~.Domain.dtype`"""
+        return self.dtype
+
+    def __call__(self, *args):
+        """Construct an element of ``self`` domain from ``args``. """
+        return self.new(*args)
+
+    def normal(self, *args):
+        return self.dtype(*args)
+
+    def convert_from(self, element, base):
+        """Convert ``element`` to ``self.dtype`` given the base domain. """
+        if base.alias is not None:
+            method = "from_" + base.alias
+        else:
+            method = "from_" + base.__class__.__name__
+
+        _convert = getattr(self, method)
+
+        if _convert is not None:
+            result = _convert(element, base)
+
+            if result is not None:
+                return result
+
+        raise CoercionFailed("Cannot convert %s of type %s from %s to %s" % (element, type(element), base, self))
+
+    def convert(self, element, base=None):
+        """Convert ``element`` to ``self.dtype``. """
+
+        if base is not None:
+            if _not_a_coeff(element):
+                raise CoercionFailed('%s is not in any domain' % element)
+            return self.convert_from(element, base)
+
+        if self.of_type(element):
+            return element
+
+        if _not_a_coeff(element):
+            raise CoercionFailed('%s is not in any domain' % element)
+
+        from sympy.polys.domains import ZZ, QQ, RealField, ComplexField
+
+        if ZZ.of_type(element):
+            return self.convert_from(element, ZZ)
+
+        if isinstance(element, int):
+            return self.convert_from(ZZ(element), ZZ)
+
+        if GROUND_TYPES != 'python':
+            if isinstance(element, ZZ.tp):
+                return self.convert_from(element, ZZ)
+            if isinstance(element, QQ.tp):
+                return self.convert_from(element, QQ)
+
+        if isinstance(element, float):
+            parent = RealField(tol=False)
+            return self.convert_from(parent(element), parent)
+
+        if isinstance(element, complex):
+            parent = ComplexField(tol=False)
+            return self.convert_from(parent(element), parent)
+
+        if isinstance(element, DomainElement):
+            return self.convert_from(element, element.parent())
+
+        # TODO: implement this in from_ methods
+        if self.is_Numerical and getattr(element, 'is_ground', False):
+            return self.convert(element.LC())
+
+        if isinstance(element, Basic):
+            try:
+                return self.from_sympy(element)
+            except (TypeError, ValueError):
+                pass
+        else: # TODO: remove this branch
+            if not is_sequence(element):
+                try:
+                    element = sympify(element, strict=True)
+                    if isinstance(element, Basic):
+                        return self.from_sympy(element)
+                except (TypeError, ValueError):
+                    pass
+
+        raise CoercionFailed("Cannot convert %s of type %s to %s" % (element, type(element), self))
+
+    def of_type(self, element):
+        """Check if ``a`` is of type ``dtype``. """
+        return isinstance(element, self.tp) # XXX: this isn't correct, e.g. PolyElement
+
+    def __contains__(self, a):
+        """Check if ``a`` belongs to this domain. """
+        try:
+            if _not_a_coeff(a):
+                raise CoercionFailed
+            self.convert(a)  # this might raise, too
+        except CoercionFailed:
+            return False
+
+        return True
+
+    def to_sympy(self, a):
+        """Convert domain element *a* to a SymPy expression (Expr).
+
+        Explanation
+        ===========
+
+        Convert a :py:class:`~.Domain` element *a* to :py:class:`~.Expr`. Most
+        public SymPy functions work with objects of type :py:class:`~.Expr`.
+        The elements of a :py:class:`~.Domain` have a different internal
+        representation. It is not possible to mix domain elements with
+        :py:class:`~.Expr` so each domain has :py:meth:`~.Domain.to_sympy` and
+        :py:meth:`~.Domain.from_sympy` methods to convert its domain elements
+        to and from :py:class:`~.Expr`.
+
+        Parameters
+        ==========
+
+        a: domain element
+            An element of this :py:class:`~.Domain`.
+
+        Returns
+        =======
+
+        expr: Expr
+            A normal SymPy expression of type :py:class:`~.Expr`.
+
+        Examples
+        ========
+
+        Construct an element of the :ref:`QQ` domain and then convert it to
+        :py:class:`~.Expr`.
+
+        >>> from sympy import QQ, Expr
+        >>> q_domain = QQ(2)
+        >>> q_domain
+        2
+        >>> q_expr = QQ.to_sympy(q_domain)
+        >>> q_expr
+        2
+
+        Although the printed forms look similar these objects are not of the
+        same type.
+
+        >>> isinstance(q_domain, Expr)
+        False
+        >>> isinstance(q_expr, Expr)
+        True
+
+        Construct an element of :ref:`K[x]` and convert to
+        :py:class:`~.Expr`.
+
+        >>> from sympy import Symbol
+        >>> x = Symbol('x')
+        >>> K = QQ[x]
+        >>> x_domain = K.gens[0]  # generator x as a domain element
+        >>> p_domain = x_domain**2/3 + 1
+        >>> p_domain
+        1/3*x**2 + 1
+        >>> p_expr = K.to_sympy(p_domain)
+        >>> p_expr
+        x**2/3 + 1
+
+        The :py:meth:`~.Domain.from_sympy` method is used for the opposite
+        conversion from a normal SymPy expression to a domain element.
+
+        >>> p_domain == p_expr
+        False
+        >>> K.from_sympy(p_expr) == p_domain
+        True
+        >>> K.to_sympy(p_domain) == p_expr
+        True
+        >>> K.from_sympy(K.to_sympy(p_domain)) == p_domain
+        True
+        >>> K.to_sympy(K.from_sympy(p_expr)) == p_expr
+        True
+
+        The :py:meth:`~.Domain.from_sympy` method makes it easier to construct
+        domain elements interactively.
+
+        >>> from sympy import Symbol
+        >>> x = Symbol('x')
+        >>> K = QQ[x]
+        >>> K.from_sympy(x**2/3 + 1)
+        1/3*x**2 + 1
+
+        See also
+        ========
+
+        from_sympy
+        convert_from
+        """
+        raise NotImplementedError
+
+    def from_sympy(self, a):
+        """Convert a SymPy expression to an element of this domain.
+
+        Explanation
+        ===========
+
+        See :py:meth:`~.Domain.to_sympy` for explanation and examples.
+
+        Parameters
+        ==========
+
+        expr: Expr
+            A normal SymPy expression of type :py:class:`~.Expr`.
+
+        Returns
+        =======
+
+        a: domain element
+            An element of this :py:class:`~.Domain`.
+
+        See also
+        ========
+
+        to_sympy
+        convert_from
+        """
+        raise NotImplementedError
+
+    def sum(self, args):
+        return sum(args, start=self.zero)
+
+    def from_FF(K1, a, K0):
+        """Convert ``ModularInteger(int)`` to ``dtype``. """
+        return None
+
+    def from_FF_python(K1, a, K0):
+        """Convert ``ModularInteger(int)`` to ``dtype``. """
+        return None
+
+    def from_ZZ_python(K1, a, K0):
+        """Convert a Python ``int`` object to ``dtype``. """
+        return None
+
+    def from_QQ_python(K1, a, K0):
+        """Convert a Python ``Fraction`` object to ``dtype``. """
+        return None
+
+    def from_FF_gmpy(K1, a, K0):
+        """Convert ``ModularInteger(mpz)`` to ``dtype``. """
+        return None
+
+    def from_ZZ_gmpy(K1, a, K0):
+        """Convert a GMPY ``mpz`` object to ``dtype``. """
+        return None
+
+    def from_QQ_gmpy(K1, a, K0):
+        """Convert a GMPY ``mpq`` object to ``dtype``. """
+        return None
+
+    def from_RealField(K1, a, K0):
+        """Convert a real element object to ``dtype``. """
+        return None
+
+    def from_ComplexField(K1, a, K0):
+        """Convert a complex element to ``dtype``. """
+        return None
+
+    def from_AlgebraicField(K1, a, K0):
+        """Convert an algebraic number to ``dtype``. """
+        return None
+
+    def from_PolynomialRing(K1, a, K0):
+        """Convert a polynomial to ``dtype``. """
+        if a.is_ground:
+            return K1.convert(a.LC, K0.dom)
+
+    def from_FractionField(K1, a, K0):
+        """Convert a rational function to ``dtype``. """
+        return None
+
+    def from_MonogenicFiniteExtension(K1, a, K0):
+        """Convert an ``ExtensionElement`` to ``dtype``. """
+        return K1.convert_from(a.rep, K0.ring)
+
+    def from_ExpressionDomain(K1, a, K0):
+        """Convert a ``EX`` object to ``dtype``. """
+        return K1.from_sympy(a.ex)
+
+    def from_ExpressionRawDomain(K1, a, K0):
+        """Convert a ``EX`` object to ``dtype``. """
+        return K1.from_sympy(a)
+
+    def from_GlobalPolynomialRing(K1, a, K0):
+        """Convert a polynomial to ``dtype``. """
+        if a.degree() <= 0:
+            return K1.convert(a.LC(), K0.dom)
+
+    def from_GeneralizedPolynomialRing(K1, a, K0):
+        return K1.from_FractionField(a, K0)
+
+    def unify_with_symbols(K0, K1, symbols):
+        if (K0.is_Composite and (set(K0.symbols) & set(symbols))) or (K1.is_Composite and (set(K1.symbols) & set(symbols))):
+            raise UnificationFailed("Cannot unify %s with %s, given %s generators" % (K0, K1, tuple(symbols)))
+
+        return K0.unify(K1)
+
+    def unify_composite(K0, K1):
+        """Unify two domains where at least one is composite."""
+        K0_ground = K0.dom if K0.is_Composite else K0
+        K1_ground = K1.dom if K1.is_Composite else K1
+
+        K0_symbols = K0.symbols if K0.is_Composite else ()
+        K1_symbols = K1.symbols if K1.is_Composite else ()
+
+        domain = K0_ground.unify(K1_ground)
+        symbols = _unify_gens(K0_symbols, K1_symbols)
+        order = K0.order if K0.is_Composite else K1.order
+
+        # E.g. ZZ[x].unify(QQ.frac_field(x)) -> ZZ.frac_field(x)
+        if ((K0.is_FractionField and K1.is_PolynomialRing or
+             K1.is_FractionField and K0.is_PolynomialRing) and
+             (not K0_ground.is_Field or not K1_ground.is_Field) and domain.is_Field
+             and domain.has_assoc_Ring):
+            domain = domain.get_ring()
+
+        if K0.is_Composite and (not K1.is_Composite or K0.is_FractionField or K1.is_PolynomialRing):
+            cls = K0.__class__
+        else:
+            cls = K1.__class__
+
+        # Here cls might be PolynomialRing, FractionField, GlobalPolynomialRing
+        # (dense/old Polynomialring) or dense/old FractionField.
+
+        from sympy.polys.domains.old_polynomialring import GlobalPolynomialRing
+        if cls == GlobalPolynomialRing:
+            return cls(domain, symbols)
+
+        return cls(domain, symbols, order)
+
+    def unify(K0, K1, symbols=None):
+        """
+        Construct a minimal domain that contains elements of ``K0`` and ``K1``.
+
+        Known domains (from smallest to largest):
+
+        - ``GF(p)``
+        - ``ZZ``
+        - ``QQ``
+        - ``RR(prec, tol)``
+        - ``CC(prec, tol)``
+        - ``ALG(a, b, c)``
+        - ``K[x, y, z]``
+        - ``K(x, y, z)``
+        - ``EX``
+
+        """
+        if symbols is not None:
+            return K0.unify_with_symbols(K1, symbols)
+
+        if K0 == K1:
+            return K0
+
+        if not (K0.has_CharacteristicZero and K1.has_CharacteristicZero):
+            # Reject unification of domains with different characteristics.
+            if K0.characteristic() != K1.characteristic():
+                raise UnificationFailed("Cannot unify %s with %s" % (K0, K1))
+
+            # We do not get here if K0 == K1. The two domains have the same
+            # characteristic but are unequal so at least one is composite and
+            # we are unifying something like GF(3).unify(GF(3)[x]).
+            return K0.unify_composite(K1)
+
+        # From here we know both domains have characteristic zero and it can be
+        # acceptable to fall back on EX.
+
+        if K0.is_EXRAW:
+            return K0
+        if K1.is_EXRAW:
+            return K1
+
+        if K0.is_EX:
+            return K0
+        if K1.is_EX:
+            return K1
+
+        if K0.is_FiniteExtension or K1.is_FiniteExtension:
+            if K1.is_FiniteExtension:
+                K0, K1 = K1, K0
+            if K1.is_FiniteExtension:
+                # Unifying two extensions.
+                # Try to ensure that K0.unify(K1) == K1.unify(K0)
+                if list(ordered([K0.modulus, K1.modulus]))[1] == K0.modulus:
+                    K0, K1 = K1, K0
+                return K1.set_domain(K0)
+            else:
+                # Drop the generator from other and unify with the base domain
+                K1 = K1.drop(K0.symbol)
+                K1 = K0.domain.unify(K1)
+                return K0.set_domain(K1)
+
+        if K0.is_Composite or K1.is_Composite:
+            return K0.unify_composite(K1)
+
+        def mkinexact(cls, K0, K1):
+            prec = max(K0.precision, K1.precision)
+            tol = max(K0.tolerance, K1.tolerance)
+            return cls(prec=prec, tol=tol)
+
+        if K1.is_ComplexField:
+            K0, K1 = K1, K0
+        if K0.is_ComplexField:
+            if K1.is_ComplexField or K1.is_RealField:
+                return mkinexact(K0.__class__, K0, K1)
+            else:
+                return K0
+
+        if K1.is_RealField:
+            K0, K1 = K1, K0
+        if K0.is_RealField:
+            if K1.is_RealField:
+                return mkinexact(K0.__class__, K0, K1)
+            elif K1.is_GaussianRing or K1.is_GaussianField:
+                from sympy.polys.domains.complexfield import ComplexField
+                return ComplexField(prec=K0.precision, tol=K0.tolerance)
+            else:
+                return K0
+
+        if K1.is_AlgebraicField:
+            K0, K1 = K1, K0
+        if K0.is_AlgebraicField:
+            if K1.is_GaussianRing:
+                K1 = K1.get_field()
+            if K1.is_GaussianField:
+                K1 = K1.as_AlgebraicField()
+            if K1.is_AlgebraicField:
+                return K0.__class__(K0.dom.unify(K1.dom), *_unify_gens(K0.orig_ext, K1.orig_ext))
+            else:
+                return K0
+
+        if K0.is_GaussianField:
+            return K0
+        if K1.is_GaussianField:
+            return K1
+
+        if K0.is_GaussianRing:
+            if K1.is_RationalField:
+                K0 = K0.get_field()
+            return K0
+        if K1.is_GaussianRing:
+            if K0.is_RationalField:
+                K1 = K1.get_field()
+            return K1
+
+        if K0.is_RationalField:
+            return K0
+        if K1.is_RationalField:
+            return K1
+
+        if K0.is_IntegerRing:
+            return K0
+        if K1.is_IntegerRing:
+            return K1
+
+        from sympy.polys.domains import EX
+        return EX
+
+    def __eq__(self, other):
+        """Returns ``True`` if two domains are equivalent. """
+        # XXX: Remove this.
+        return isinstance(other, Domain) and self.dtype == other.dtype
+
+    def __ne__(self, other):
+        """Returns ``False`` if two domains are equivalent. """
+        return not self == other
+
+    def map(self, seq):
+        """Rersively apply ``self`` to all elements of ``seq``. """
+        result = []
+
+        for elt in seq:
+            if isinstance(elt, list):
+                result.append(self.map(elt))
+            else:
+                result.append(self(elt))
+
+        return result
+
+    def get_ring(self):
+        """Returns a ring associated with ``self``. """
+        raise DomainError('there is no ring associated with %s' % self)
+
+    def get_field(self):
+        """Returns a field associated with ``self``. """
+        raise DomainError('there is no field associated with %s' % self)
+
+    def get_exact(self):
+        """Returns an exact domain associated with ``self``. """
+        return self
+
+    def __getitem__(self, symbols):
+        """The mathematical way to make a polynomial ring. """
+        if hasattr(symbols, '__iter__'):
+            return self.poly_ring(*symbols)
+        else:
+            return self.poly_ring(symbols)
+
+    def poly_ring(self, *symbols, order=lex):
+        """Returns a polynomial ring, i.e. `K[X]`. """
+        from sympy.polys.domains.polynomialring import PolynomialRing
+        return PolynomialRing(self, symbols, order)
+
+    def frac_field(self, *symbols, order=lex):
+        """Returns a fraction field, i.e. `K(X)`. """
+        from sympy.polys.domains.fractionfield import FractionField
+        return FractionField(self, symbols, order)
+
+    def old_poly_ring(self, *symbols, **kwargs):
+        """Returns a polynomial ring, i.e. `K[X]`. """
+        from sympy.polys.domains.old_polynomialring import PolynomialRing
+        return PolynomialRing(self, *symbols, **kwargs)
+
+    def old_frac_field(self, *symbols, **kwargs):
+        """Returns a fraction field, i.e. `K(X)`. """
+        from sympy.polys.domains.old_fractionfield import FractionField
+        return FractionField(self, *symbols, **kwargs)
+
+    def algebraic_field(self, *extension, alias=None):
+        r"""Returns an algebraic field, i.e. `K(\alpha, \ldots)`. """
+        raise DomainError("Cannot create algebraic field over %s" % self)
+
+    def alg_field_from_poly(self, poly, alias=None, root_index=-1):
+        r"""
+        Convenience method to construct an algebraic extension on a root of a
+        polynomial, chosen by root index.
+
+        Parameters
+        ==========
+
+        poly : :py:class:`~.Poly`
+            The polynomial whose root generates the extension.
+        alias : str, optional (default=None)
+            Symbol name for the generator of the extension.
+            E.g. "alpha" or "theta".
+        root_index : int, optional (default=-1)
+            Specifies which root of the polynomial is desired. The ordering is
+            as defined by the :py:class:`~.ComplexRootOf` class. The default of
+            ``-1`` selects the most natural choice in the common cases of
+            quadratic and cyclotomic fields (the square root on the positive
+            real or imaginary axis, resp. $\mathrm{e}^{2\pi i/n}$).
+
+        Examples
+        ========
+
+        >>> from sympy import QQ, Poly
+        >>> from sympy.abc import x
+        >>> f = Poly(x**2 - 2)
+        >>> K = QQ.alg_field_from_poly(f)
+        >>> K.ext.minpoly == f
+        True
+        >>> g = Poly(8*x**3 - 6*x - 1)
+        >>> L = QQ.alg_field_from_poly(g, "alpha")
+        >>> L.ext.minpoly == g
+        True
+        >>> L.to_sympy(L([1, 1, 1]))
+        alpha**2 + alpha + 1
+
+        """
+        from sympy.polys.rootoftools import CRootOf
+        root = CRootOf(poly, root_index)
+        alpha = AlgebraicNumber(root, alias=alias)
+        return self.algebraic_field(alpha, alias=alias)
+
+    def cyclotomic_field(self, n, ss=False, alias="zeta", gen=None, root_index=-1):
+        r"""
+        Convenience method to construct a cyclotomic field.
+
+        Parameters
+        ==========
+
+        n : int
+            Construct the nth cyclotomic field.
+        ss : boolean, optional (default=False)
+            If True, append *n* as a subscript on the alias string.
+        alias : str, optional (default="zeta")
+            Symbol name for the generator.
+        gen : :py:class:`~.Symbol`, optional (default=None)
+            Desired variable for the cyclotomic polynomial that defines the
+            field. If ``None``, a dummy variable will be used.
+        root_index : int, optional (default=-1)
+            Specifies which root of the polynomial is desired. The ordering is
+            as defined by the :py:class:`~.ComplexRootOf` class. The default of
+            ``-1`` selects the root $\mathrm{e}^{2\pi i/n}$.
+
+        Examples
+        ========
+
+        >>> from sympy import QQ, latex
+        >>> K = QQ.cyclotomic_field(5)
+        >>> K.to_sympy(K([-1, 1]))
+        1 - zeta
+        >>> L = QQ.cyclotomic_field(7, True)
+        >>> a = L.to_sympy(L([-1, 1]))
+        >>> print(a)
+        1 - zeta7
+        >>> print(latex(a))
+        1 - \zeta_{7}
+
+        """
+        from sympy.polys.specialpolys import cyclotomic_poly
+        if ss:
+            alias += str(n)
+        return self.alg_field_from_poly(cyclotomic_poly(n, gen), alias=alias,
+                                        root_index=root_index)
+
+    def inject(self, *symbols):
+        """Inject generators into this domain. """
+        raise NotImplementedError
+
+    def drop(self, *symbols):
+        """Drop generators from this domain. """
+        if self.is_Simple:
+            return self
+        raise NotImplementedError  # pragma: no cover
+
+    def is_zero(self, a):
+        """Returns True if ``a`` is zero. """
+        return not a
+
+    def is_one(self, a):
+        """Returns True if ``a`` is one. """
+        return a == self.one
+
+    def is_positive(self, a):
+        """Returns True if ``a`` is positive. """
+        return a > 0
+
+    def is_negative(self, a):
+        """Returns True if ``a`` is negative. """
+        return a < 0
+
+    def is_nonpositive(self, a):
+        """Returns True if ``a`` is non-positive. """
+        return a <= 0
+
+    def is_nonnegative(self, a):
+        """Returns True if ``a`` is non-negative. """
+        return a >= 0
+
+    def canonical_unit(self, a):
+        if self.is_negative(a):
+            return -self.one
+        else:
+            return self.one
+
+    def abs(self, a):
+        """Absolute value of ``a``, implies ``__abs__``. """
+        return abs(a)
+
+    def neg(self, a):
+        """Returns ``a`` negated, implies ``__neg__``. """
+        return -a
+
+    def pos(self, a):
+        """Returns ``a`` positive, implies ``__pos__``. """
+        return +a
+
+    def add(self, a, b):
+        """Sum of ``a`` and ``b``, implies ``__add__``.  """
+        return a + b
+
+    def sub(self, a, b):
+        """Difference of ``a`` and ``b``, implies ``__sub__``.  """
+        return a - b
+
+    def mul(self, a, b):
+        """Product of ``a`` and ``b``, implies ``__mul__``.  """
+        return a * b
+
+    def pow(self, a, b):
+        """Raise ``a`` to power ``b``, implies ``__pow__``.  """
+        return a ** b
+
+    def exquo(self, a, b):
+        """Exact quotient of *a* and *b*. Analogue of ``a / b``.
+
+        Explanation
+        ===========
+
+        This is essentially the same as ``a / b`` except that an error will be
+        raised if the division is inexact (if there is any remainder) and the
+        result will always be a domain element. When working in a
+        :py:class:`~.Domain` that is not a :py:class:`~.Field` (e.g. :ref:`ZZ`
+        or :ref:`K[x]`) ``exquo`` should be used instead of ``/``.
+
+        The key invariant is that if ``q = K.exquo(a, b)`` (and ``exquo`` does
+        not raise an exception) then ``a == b*q``.
+
+        Examples
+        ========
+
+        We can use ``K.exquo`` instead of ``/`` for exact division.
+
+        >>> from sympy import ZZ
+        >>> ZZ.exquo(ZZ(4), ZZ(2))
+        2
+        >>> ZZ.exquo(ZZ(5), ZZ(2))
+        Traceback (most recent call last):
+            ...
+        ExactQuotientFailed: 2 does not divide 5 in ZZ
+
+        Over a :py:class:`~.Field` such as :ref:`QQ`, division (with nonzero
+        divisor) is always exact so in that case ``/`` can be used instead of
+        :py:meth:`~.Domain.exquo`.
+
+        >>> from sympy import QQ
+        >>> QQ.exquo(QQ(5), QQ(2))
+        5/2
+        >>> QQ(5) / QQ(2)
+        5/2
+
+        Parameters
+        ==========
+
+        a: domain element
+            The dividend
+        b: domain element
+            The divisor
+
+        Returns
+        =======
+
+        q: domain element
+            The exact quotient
+
+        Raises
+        ======
+
+        ExactQuotientFailed: if exact division is not possible.
+        ZeroDivisionError: when the divisor is zero.
+
+        See also
+        ========
+
+        quo: Analogue of ``a // b``
+        rem: Analogue of ``a % b``
+        div: Analogue of ``divmod(a, b)``
+
+        Notes
+        =====
+
+        Since the default :py:attr:`~.Domain.dtype` for :ref:`ZZ` is ``int``
+        (or ``mpz``) division as ``a / b`` should not be used as it would give
+        a ``float`` which is not a domain element.
+
+        >>> ZZ(4) / ZZ(2) # doctest: +SKIP
+        2.0
+        >>> ZZ(5) / ZZ(2) # doctest: +SKIP
+        2.5
+
+        On the other hand with `SYMPY_GROUND_TYPES=flint` elements of :ref:`ZZ`
+        are ``flint.fmpz`` and division would raise an exception:
+
+        >>> ZZ(4) / ZZ(2) # doctest: +SKIP
+        Traceback (most recent call last):
+        ...
+        TypeError: unsupported operand type(s) for /: 'fmpz' and 'fmpz'
+
+        Using ``/`` with :ref:`ZZ` will lead to incorrect results so
+        :py:meth:`~.Domain.exquo` should be used instead.
+
+        """
+        raise NotImplementedError
+
+    def quo(self, a, b):
+        """Quotient of *a* and *b*. Analogue of ``a // b``.
+
+        ``K.quo(a, b)`` is equivalent to ``K.div(a, b)[0]``. See
+        :py:meth:`~.Domain.div` for more explanation.
+
+        See also
+        ========
+
+        rem: Analogue of ``a % b``
+        div: Analogue of ``divmod(a, b)``
+        exquo: Analogue of ``a / b``
+        """
+        raise NotImplementedError
+
+    def rem(self, a, b):
+        """Modulo division of *a* and *b*. Analogue of ``a % b``.
+
+        ``K.rem(a, b)`` is equivalent to ``K.div(a, b)[1]``. See
+        :py:meth:`~.Domain.div` for more explanation.
+
+        See also
+        ========
+
+        quo: Analogue of ``a // b``
+        div: Analogue of ``divmod(a, b)``
+        exquo: Analogue of ``a / b``
+        """
+        raise NotImplementedError
+
+    def div(self, a, b):
+        """Quotient and remainder for *a* and *b*. Analogue of ``divmod(a, b)``
+
+        Explanation
+        ===========
+
+        This is essentially the same as ``divmod(a, b)`` except that is more
+        consistent when working over some :py:class:`~.Field` domains such as
+        :ref:`QQ`. When working over an arbitrary :py:class:`~.Domain` the
+        :py:meth:`~.Domain.div` method should be used instead of ``divmod``.
+
+        The key invariant is that if ``q, r = K.div(a, b)`` then
+        ``a == b*q + r``.
+
+        The result of ``K.div(a, b)`` is the same as the tuple
+        ``(K.quo(a, b), K.rem(a, b))`` except that if both quotient and
+        remainder are needed then it is more efficient to use
+        :py:meth:`~.Domain.div`.
+
+        Examples
+        ========
+
+        We can use ``K.div`` instead of ``divmod`` for floor division and
+        remainder.
+
+        >>> from sympy import ZZ, QQ
+        >>> ZZ.div(ZZ(5), ZZ(2))
+        (2, 1)
+
+        If ``K`` is a :py:class:`~.Field` then the division is always exact
+        with a remainder of :py:attr:`~.Domain.zero`.
+
+        >>> QQ.div(QQ(5), QQ(2))
+        (5/2, 0)
+
+        Parameters
+        ==========
+
+        a: domain element
+            The dividend
+        b: domain element
+            The divisor
+
+        Returns
+        =======
+
+        (q, r): tuple of domain elements
+            The quotient and remainder
+
+        Raises
+        ======
+
+        ZeroDivisionError: when the divisor is zero.
+
+        See also
+        ========
+
+        quo: Analogue of ``a // b``
+        rem: Analogue of ``a % b``
+        exquo: Analogue of ``a / b``
+
+        Notes
+        =====
+
+        If ``gmpy`` is installed then the ``gmpy.mpq`` type will be used as
+        the :py:attr:`~.Domain.dtype` for :ref:`QQ`. The ``gmpy.mpq`` type
+        defines ``divmod`` in a way that is undesirable so
+        :py:meth:`~.Domain.div` should be used instead of ``divmod``.
+
+        >>> a = QQ(1)
+        >>> b = QQ(3, 2)
+        >>> a               # doctest: +SKIP
+        mpq(1,1)
+        >>> b               # doctest: +SKIP
+        mpq(3,2)
+        >>> divmod(a, b)    # doctest: +SKIP
+        (mpz(0), mpq(1,1))
+        >>> QQ.div(a, b)    # doctest: +SKIP
+        (mpq(2,3), mpq(0,1))
+
+        Using ``//`` or ``%`` with :ref:`QQ` will lead to incorrect results so
+        :py:meth:`~.Domain.div` should be used instead.
+
+        """
+        raise NotImplementedError
+
+    def invert(self, a, b):
+        """Returns inversion of ``a mod b``, implies something. """
+        raise NotImplementedError
+
+    def revert(self, a):
+        """Returns ``a**(-1)`` if possible. """
+        raise NotImplementedError
+
+    def numer(self, a):
+        """Returns numerator of ``a``. """
+        raise NotImplementedError
+
+    def denom(self, a):
+        """Returns denominator of ``a``. """
+        raise NotImplementedError
+
+    def half_gcdex(self, a, b):
+        """Half extended GCD of ``a`` and ``b``. """
+        s, t, h = self.gcdex(a, b)
+        return s, h
+
+    def gcdex(self, a, b):
+        """Extended GCD of ``a`` and ``b``. """
+        raise NotImplementedError
+
+    def cofactors(self, a, b):
+        """Returns GCD and cofactors of ``a`` and ``b``. """
+        gcd = self.gcd(a, b)
+        cfa = self.quo(a, gcd)
+        cfb = self.quo(b, gcd)
+        return gcd, cfa, cfb
+
+    def gcd(self, a, b):
+        """Returns GCD of ``a`` and ``b``. """
+        raise NotImplementedError
+
+    def lcm(self, a, b):
+        """Returns LCM of ``a`` and ``b``. """
+        raise NotImplementedError
+
+    def log(self, a, b):
+        """Returns b-base logarithm of ``a``. """
+        raise NotImplementedError
+
+    def sqrt(self, a):
+        """Returns a (possibly inexact) square root of ``a``.
+
+        Explanation
+        ===========
+        There is no universal definition of "inexact square root" for all
+        domains. It is not recommended to implement this method for domains
+        other then :ref:`ZZ`.
+
+        See also
+        ========
+        exsqrt
+        """
+        raise NotImplementedError
+
+    def is_square(self, a):
+        """Returns whether ``a`` is a square in the domain.
+
+        Explanation
+        ===========
+        Returns ``True`` if there is an element ``b`` in the domain such that
+        ``b * b == a``, otherwise returns ``False``. For inexact domains like
+        :ref:`RR` and :ref:`CC`, a tiny difference in this equality can be
+        tolerated.
+
+        See also
+        ========
+        exsqrt
+        """
+        raise NotImplementedError
+
+    def exsqrt(self, a):
+        """Principal square root of a within the domain if ``a`` is square.
+
+        Explanation
+        ===========
+        The implementation of this method should return an element ``b`` in the
+        domain such that ``b * b == a``, or ``None`` if there is no such ``b``.
+        For inexact domains like :ref:`RR` and :ref:`CC`, a tiny difference in
+        this equality can be tolerated. The choice of a "principal" square root
+        should follow a consistent rule whenever possible.
+
+        See also
+        ========
+        sqrt, is_square
+        """
+        raise NotImplementedError
+
+    def evalf(self, a, prec=None, **options):
+        """Returns numerical approximation of ``a``. """
+        return self.to_sympy(a).evalf(prec, **options)
+
+    n = evalf
+
+    def real(self, a):
+        return a
+
+    def imag(self, a):
+        return self.zero
+
+    def almosteq(self, a, b, tolerance=None):
+        """Check if ``a`` and ``b`` are almost equal. """
+        return a == b
+
+    def characteristic(self):
+        """Return the characteristic of this domain. """
+        raise NotImplementedError('characteristic()')
+
+
+__all__ = ['Domain']

lib/python3.12/site-packages/sympy/polys/domains/domainelement.py

@@ -0,0 +1,38 @@
+"""Trait for implementing domain elements. """
+
+
+from sympy.utilities import public
+
+@public
+class DomainElement:
+    """
+    Represents an element of a domain.
+
+    Mix in this trait into a class whose instances should be recognized as
+    elements of a domain. Method ``parent()`` gives that domain.
+    """
+
+    __slots__ = ()
+
+    def parent(self):
+        """Get the domain associated with ``self``
+
+        Examples
+        ========
+
+        >>> from sympy import ZZ, symbols
+        >>> x, y = symbols('x, y')
+        >>> K = ZZ[x,y]
+        >>> p = K(x)**2 + K(y)**2
+        >>> p
+        x**2 + y**2
+        >>> p.parent()
+        ZZ[x,y]
+
+        Notes
+        =====
+
+        This is used by :py:meth:`~.Domain.convert` to identify the domain
+        associated with a domain element.
+        """
+        raise NotImplementedError("abstract method")

lib/python3.12/site-packages/sympy/polys/domains/expressiondomain.py

@@ -0,0 +1,278 @@
+"""Implementation of :class:`ExpressionDomain` class. """
+
+
+from sympy.core import sympify, SympifyError
+from sympy.polys.domains.domainelement import DomainElement
+from sympy.polys.domains.characteristiczero import CharacteristicZero
+from sympy.polys.domains.field import Field
+from sympy.polys.domains.simpledomain import SimpleDomain
+from sympy.polys.polyutils import PicklableWithSlots
+from sympy.utilities import public
+
+eflags = {"deep": False, "mul": True, "power_exp": False, "power_base": False,
+              "basic": False, "multinomial": False, "log": False}
+
+@public
+class ExpressionDomain(Field, CharacteristicZero, SimpleDomain):
+    """A class for arbitrary expressions. """
+
+    is_SymbolicDomain = is_EX = True
+
+    class Expression(DomainElement, PicklableWithSlots):
+        """An arbitrary expression. """
+
+        __slots__ = ('ex',)
+
+        def __init__(self, ex):
+            if not isinstance(ex, self.__class__):
+                self.ex = sympify(ex)
+            else:
+                self.ex = ex.ex
+
+        def __repr__(f):
+            return 'EX(%s)' % repr(f.ex)
+
+        def __str__(f):
+            return 'EX(%s)' % str(f.ex)
+
+        def __hash__(self):
+            return hash((self.__class__.__name__, self.ex))
+
+        def parent(self):
+            return EX
+
+        def as_expr(f):
+            return f.ex
+
+        def numer(f):
+            return f.__class__(f.ex.as_numer_denom()[0])
+
+        def denom(f):
+            return f.__class__(f.ex.as_numer_denom()[1])
+
+        def simplify(f, ex):
+            return f.__class__(ex.cancel().expand(**eflags))
+
+        def __abs__(f):
+            return f.__class__(abs(f.ex))
+
+        def __neg__(f):
+            return f.__class__(-f.ex)
+
+        def _to_ex(f, g):
+            try:
+                return f.__class__(g)
+            except SympifyError:
+                return None
+
+        def __lt__(f, g):
+            return f.ex.sort_key() < g.ex.sort_key()
+
+        def __add__(f, g):
+            g = f._to_ex(g)
+
+            if g is None:
+                return NotImplemented
+            elif g == EX.zero:
+                return f
+            elif f == EX.zero:
+                return g
+            else:
+                return f.simplify(f.ex + g.ex)
+
+        def __radd__(f, g):
+            return f.simplify(f.__class__(g).ex + f.ex)
+
+        def __sub__(f, g):
+            g = f._to_ex(g)
+
+            if g is None:
+                return NotImplemented
+            elif g == EX.zero:
+                return f
+            elif f == EX.zero:
+                return -g
+            else:
+                return f.simplify(f.ex - g.ex)
+
+        def __rsub__(f, g):
+            return f.simplify(f.__class__(g).ex - f.ex)
+
+        def __mul__(f, g):
+            g = f._to_ex(g)
+
+            if g is None:
+                return NotImplemented
+
+            if EX.zero in (f, g):
+                return EX.zero
+            elif f.ex.is_Number and g.ex.is_Number:
+                return f.__class__(f.ex*g.ex)
+
+            return f.simplify(f.ex*g.ex)
+
+        def __rmul__(f, g):
+            return f.simplify(f.__class__(g).ex*f.ex)
+
+        def __pow__(f, n):
+            n = f._to_ex(n)
+
+            if n is not None:
+                return f.simplify(f.ex**n.ex)
+            else:
+                return NotImplemented
+
+        def __truediv__(f, g):
+            g = f._to_ex(g)
+
+            if g is not None:
+                return f.simplify(f.ex/g.ex)
+            else:
+                return NotImplemented
+
+        def __rtruediv__(f, g):
+            return f.simplify(f.__class__(g).ex/f.ex)
+
+        def __eq__(f, g):
+            return f.ex == f.__class__(g).ex
+
+        def __ne__(f, g):
+            return not f == g
+
+        def __bool__(f):
+            return not f.ex.is_zero
+
+        def gcd(f, g):
+            from sympy.polys import gcd
+            return f.__class__(gcd(f.ex, f.__class__(g).ex))
+
+        def lcm(f, g):
+            from sympy.polys import lcm
+            return f.__class__(lcm(f.ex, f.__class__(g).ex))
+
+    dtype = Expression
+
+    zero = Expression(0)
+    one = Expression(1)
+
+    rep = 'EX'
+
+    has_assoc_Ring = False
+    has_assoc_Field = True
+
+    def __init__(self):
+        pass
+
+    def __eq__(self, other):
+        if isinstance(other, ExpressionDomain):
+            return True
+        else:
+            return NotImplemented
+
+    def __hash__(self):
+        return hash("EX")
+
+    def to_sympy(self, a):
+        """Convert ``a`` to a SymPy object. """
+        return a.as_expr()
+
+    def from_sympy(self, a):
+        """Convert SymPy's expression to ``dtype``. """
+        return self.dtype(a)
+
+    def from_ZZ(K1, a, K0):
+        """Convert a Python ``int`` object to ``dtype``. """
+        return K1(K0.to_sympy(a))
+
+    def from_ZZ_python(K1, a, K0):
+        """Convert a Python ``int`` object to ``dtype``. """
+        return K1(K0.to_sympy(a))
+
+    def from_QQ(K1, a, K0):
+        """Convert a Python ``Fraction`` object to ``dtype``. """
+        return K1(K0.to_sympy(a))
+
+    def from_QQ_python(K1, a, K0):
+        """Convert a Python ``Fraction`` object to ``dtype``. """
+        return K1(K0.to_sympy(a))
+
+    def from_ZZ_gmpy(K1, a, K0):
+        """Convert a GMPY ``mpz`` object to ``dtype``. """
+        return K1(K0.to_sympy(a))
+
+    def from_QQ_gmpy(K1, a, K0):
+        """Convert a GMPY ``mpq`` object to ``dtype``. """
+        return K1(K0.to_sympy(a))
+
+    def from_GaussianIntegerRing(K1, a, K0):
+        """Convert a ``GaussianRational`` object to ``dtype``. """
+        return K1(K0.to_sympy(a))
+
+    def from_GaussianRationalField(K1, a, K0):
+        """Convert a ``GaussianRational`` object to ``dtype``. """
+        return K1(K0.to_sympy(a))
+
+    def from_AlgebraicField(K1, a, K0):
+        """Convert an ``ANP`` object to ``dtype``. """
+        return K1(K0.to_sympy(a))
+
+    def from_RealField(K1, a, K0):
+        """Convert a mpmath ``mpf`` object to ``dtype``. """
+        return K1(K0.to_sympy(a))
+
+    def from_ComplexField(K1, a, K0):
+        """Convert a mpmath ``mpc`` object to ``dtype``. """
+        return K1(K0.to_sympy(a))
+
+    def from_PolynomialRing(K1, a, K0):
+        """Convert a ``DMP`` object to ``dtype``. """
+        return K1(K0.to_sympy(a))
+
+    def from_FractionField(K1, a, K0):
+        """Convert a ``DMF`` object to ``dtype``. """
+        return K1(K0.to_sympy(a))
+
+    def from_ExpressionDomain(K1, a, K0):
+        """Convert a ``EX`` object to ``dtype``. """
+        return a
+
+    def get_ring(self):
+        """Returns a ring associated with ``self``. """
+        return self  # XXX: EX is not a ring but we don't have much choice here.
+
+    def get_field(self):
+        """Returns a field associated with ``self``. """
+        return self
+
+    def is_positive(self, a):
+        """Returns True if ``a`` is positive. """
+        return a.ex.as_coeff_mul()[0].is_positive
+
+    def is_negative(self, a):
+        """Returns True if ``a`` is negative. """
+        return a.ex.could_extract_minus_sign()
+
+    def is_nonpositive(self, a):
+        """Returns True if ``a`` is non-positive. """
+        return a.ex.as_coeff_mul()[0].is_nonpositive
+
+    def is_nonnegative(self, a):
+        """Returns True if ``a`` is non-negative. """
+        return a.ex.as_coeff_mul()[0].is_nonnegative
+
+    def numer(self, a):
+        """Returns numerator of ``a``. """
+        return a.numer()
+
+    def denom(self, a):
+        """Returns denominator of ``a``. """
+        return a.denom()
+
+    def gcd(self, a, b):
+        return self(1)
+
+    def lcm(self, a, b):
+        return a.lcm(b)
+
+
+EX = ExpressionDomain()

lib/python3.12/site-packages/sympy/polys/domains/expressionrawdomain.py

@@ -0,0 +1,57 @@
+"""Implementation of :class:`ExpressionRawDomain` class. """
+
+
+from sympy.core import Expr, S, sympify, Add
+from sympy.polys.domains.characteristiczero import CharacteristicZero
+from sympy.polys.domains.field import Field
+from sympy.polys.domains.simpledomain import SimpleDomain
+from sympy.polys.polyerrors import CoercionFailed
+from sympy.utilities import public
+
+
+@public
+class ExpressionRawDomain(Field, CharacteristicZero, SimpleDomain):
+    """A class for arbitrary expressions but without automatic simplification. """
+
+    is_SymbolicRawDomain = is_EXRAW = True
+
+    dtype = Expr
+
+    zero = S.Zero
+    one = S.One
+
+    rep = 'EXRAW'
+
+    has_assoc_Ring = False
+    has_assoc_Field = True
+
+    def __init__(self):
+        pass
+
+    @classmethod
+    def new(self, a):
+        return sympify(a)
+
+    def to_sympy(self, a):
+        """Convert ``a`` to a SymPy object. """
+        return a
+
+    def from_sympy(self, a):
+        """Convert SymPy's expression to ``dtype``. """
+        if not isinstance(a, Expr):
+            raise CoercionFailed(f"Expecting an Expr instance but found: {type(a).__name__}")
+        return a
+
+    def convert_from(self, a, K):
+        """Convert a domain element from another domain to EXRAW"""
+        return K.to_sympy(a)
+
+    def get_field(self):
+        """Returns a field associated with ``self``. """
+        return self
+
+    def sum(self, items):
+        return Add(*items)
+
+
+EXRAW = ExpressionRawDomain()

lib/python3.12/site-packages/sympy/polys/domains/field.py

@@ -0,0 +1,104 @@
+"""Implementation of :class:`Field` class. """
+
+
+from sympy.polys.domains.ring import Ring
+from sympy.polys.polyerrors import NotReversible, DomainError
+from sympy.utilities import public
+
+@public
+class Field(Ring):
+    """Represents a field domain. """
+
+    is_Field = True
+    is_PID = True
+
+    def get_ring(self):
+        """Returns a ring associated with ``self``. """
+        raise DomainError('there is no ring associated with %s' % self)
+
+    def get_field(self):
+        """Returns a field associated with ``self``. """
+        return self
+
+    def exquo(self, a, b):
+        """Exact quotient of ``a`` and ``b``, implies ``__truediv__``.  """
+        return a / b
+
+    def quo(self, a, b):
+        """Quotient of ``a`` and ``b``, implies ``__truediv__``. """
+        return a / b
+
+    def rem(self, a, b):
+        """Remainder of ``a`` and ``b``, implies nothing.  """
+        return self.zero
+
+    def div(self, a, b):
+        """Division of ``a`` and ``b``, implies ``__truediv__``. """
+        return a / b, self.zero
+
+    def gcd(self, a, b):
+        """
+        Returns GCD of ``a`` and ``b``.
+
+        This definition of GCD over fields allows to clear denominators
+        in `primitive()`.
+
+        Examples
+        ========
+
+        >>> from sympy.polys.domains import QQ
+        >>> from sympy import S, gcd, primitive
+        >>> from sympy.abc import x
+
+        >>> QQ.gcd(QQ(2, 3), QQ(4, 9))
+        2/9
+        >>> gcd(S(2)/3, S(4)/9)
+        2/9
+        >>> primitive(2*x/3 + S(4)/9)
+        (2/9, 3*x + 2)
+
+        """
+        try:
+            ring = self.get_ring()
+        except DomainError:
+            return self.one
+
+        p = ring.gcd(self.numer(a), self.numer(b))
+        q = ring.lcm(self.denom(a), self.denom(b))
+
+        return self.convert(p, ring)/q
+
+    def lcm(self, a, b):
+        """
+        Returns LCM of ``a`` and ``b``.
+
+        >>> from sympy.polys.domains import QQ
+        >>> from sympy import S, lcm
+
+        >>> QQ.lcm(QQ(2, 3), QQ(4, 9))
+        4/3
+        >>> lcm(S(2)/3, S(4)/9)
+        4/3
+
+        """
+
+        try:
+            ring = self.get_ring()
+        except DomainError:
+            return a*b
+
+        p = ring.lcm(self.numer(a), self.numer(b))
+        q = ring.gcd(self.denom(a), self.denom(b))
+
+        return self.convert(p, ring)/q
+
+    def revert(self, a):
+        """Returns ``a**(-1)`` if possible. """
+        if a:
+            return 1/a
+        else:
+            raise NotReversible('zero is not reversible')
+
+    def is_unit(self, a):
+        """Return true if ``a`` is a invertible"""
+        return bool(a)

lib/python3.12/site-packages/sympy/polys/domains/finitefield.py

@@ -0,0 +1,328 @@
+"""Implementation of :class:`FiniteField` class. """
+
+import operator
+
+from sympy.external.gmpy import GROUND_TYPES
+from sympy.utilities.decorator import doctest_depends_on
+
+from sympy.core.numbers import int_valued
+from sympy.polys.domains.field import Field
+
+from sympy.polys.domains.modularinteger import ModularIntegerFactory
+from sympy.polys.domains.simpledomain import SimpleDomain
+from sympy.polys.galoistools import gf_zassenhaus, gf_irred_p_rabin
+from sympy.polys.polyerrors import CoercionFailed
+from sympy.utilities import public
+from sympy.polys.domains.groundtypes import SymPyInteger
+
+
+if GROUND_TYPES == 'flint':
+    __doctest_skip__ = ['FiniteField']
+
+
+if GROUND_TYPES == 'flint':
+    import flint
+    # Don't use python-flint < 0.5.0 because nmod was missing some features in
+    # previous versions of python-flint and fmpz_mod was not yet added.
+    _major, _minor, *_ = flint.__version__.split('.')
+    if (int(_major), int(_minor)) < (0, 5):
+        flint = None
+else:
+    flint = None
+
+
+def _modular_int_factory(mod, dom, symmetric, self):
+
+    # Use flint if available
+    if flint is not None:
+
+        nmod = flint.nmod
+        fmpz_mod_ctx = flint.fmpz_mod_ctx
+        index = operator.index
+
+        try:
+            mod = dom.convert(mod)
+        except CoercionFailed:
+            raise ValueError('modulus must be an integer, got %s' % mod)
+
+        # mod might be e.g. Integer
+        try:
+            fmpz_mod_ctx(mod)
+        except TypeError:
+            mod = index(mod)
+
+        # flint's nmod is only for moduli up to 2^64-1 (on a 64-bit machine)
+        try:
+            nmod(0, mod)
+        except OverflowError:
+            # Use fmpz_mod
+            fctx = fmpz_mod_ctx(mod)
+
+            def ctx(x):
+                try:
+                    return fctx(x)
+                except TypeError:
+                    # x might be Integer
+                    return fctx(index(x))
+        else:
+            # Use nmod
+            def ctx(x):
+                try:
+                    return nmod(x, mod)
+                except TypeError:
+                    return nmod(index(x), mod)
+
+        return ctx
+
+    # Use the Python implementation
+    return ModularIntegerFactory(mod, dom, symmetric, self)
+
+
+@public
+@doctest_depends_on(modules=['python', 'gmpy'])
+class FiniteField(Field, SimpleDomain):
+    r"""Finite field of prime order :ref:`GF(p)`
+
+    A :ref:`GF(p)` domain represents a `finite field`_ `\mathbb{F}_p` of prime
+    order as :py:class:`~.Domain` in the domain system (see
+    :ref:`polys-domainsintro`).
+
+    A :py:class:`~.Poly` created from an expression with integer
+    coefficients will have the domain :ref:`ZZ`. However, if the ``modulus=p``
+    option is given then the domain will be a finite field instead.
+
+    >>> from sympy import Poly, Symbol
+    >>> x = Symbol('x')
+    >>> p = Poly(x**2 + 1)
+    >>> p
+    Poly(x**2 + 1, x, domain='ZZ')
+    >>> p.domain
+    ZZ
+    >>> p2 = Poly(x**2 + 1, modulus=2)
+    >>> p2
+    Poly(x**2 + 1, x, modulus=2)
+    >>> p2.domain
+    GF(2)
+
+    It is possible to factorise a polynomial over :ref:`GF(p)` using the
+    modulus argument to :py:func:`~.factor` or by specifying the domain
+    explicitly. The domain can also be given as a string.
+
+    >>> from sympy import factor, GF
+    >>> factor(x**2 + 1)
+    x**2 + 1
+    >>> factor(x**2 + 1, modulus=2)
+    (x + 1)**2
+    >>> factor(x**2 + 1, domain=GF(2))
+    (x + 1)**2
+    >>> factor(x**2 + 1, domain='GF(2)')
+    (x + 1)**2
+
+    It is also possible to use :ref:`GF(p)` with the :py:func:`~.cancel`
+    and :py:func:`~.gcd` functions.
+
+    >>> from sympy import cancel, gcd
+    >>> cancel((x**2 + 1)/(x + 1))
+    (x**2 + 1)/(x + 1)
+    >>> cancel((x**2 + 1)/(x + 1), domain=GF(2))
+    x + 1
+    >>> gcd(x**2 + 1, x + 1)
+    1
+    >>> gcd(x**2 + 1, x + 1, domain=GF(2))
+    x + 1
+
+    When using the domain directly :ref:`GF(p)` can be used as a constructor
+    to create instances which then support the operations ``+,-,*,**,/``
+
+    >>> from sympy import GF
+    >>> K = GF(5)
+    >>> K
+    GF(5)
+    >>> x = K(3)
+    >>> y = K(2)
+    >>> x
+    3 mod 5
+    >>> y
+    2 mod 5
+    >>> x * y
+    1 mod 5
+    >>> x / y
+    4 mod 5
+
+    Notes
+    =====
+
+    It is also possible to create a :ref:`GF(p)` domain of **non-prime**
+    order but the resulting ring is **not** a field: it is just the ring of
+    the integers modulo ``n``.
+
+    >>> K = GF(9)
+    >>> z = K(3)
+    >>> z
+    3 mod 9
+    >>> z**2
+    0 mod 9
+
+    It would be good to have a proper implementation of prime power fields
+    (``GF(p**n)``) but these are not yet implemented in SymPY.
+
+    .. _finite field: https://en.wikipedia.org/wiki/Finite_field
+    """
+
+    rep = 'FF'
+    alias = 'FF'
+
+    is_FiniteField = is_FF = True
+    is_Numerical = True
+
+    has_assoc_Ring = False
+    has_assoc_Field = True
+
+    dom = None
+    mod = None
+
+    def __init__(self, mod, symmetric=True):
+        from sympy.polys.domains import ZZ
+        dom = ZZ
+
+        if mod <= 0:
+            raise ValueError('modulus must be a positive integer, got %s' % mod)
+
+        self.dtype = _modular_int_factory(mod, dom, symmetric, self)
+        self.zero = self.dtype(0)
+        self.one = self.dtype(1)
+        self.dom = dom
+        self.mod = mod
+        self.sym = symmetric
+        self._tp = type(self.zero)
+
+    @property
+    def tp(self):
+        return self._tp
+
+    def __str__(self):
+        return 'GF(%s)' % self.mod
+
+    def __hash__(self):
+        return hash((self.__class__.__name__, self.dtype, self.mod, self.dom))
+
+    def __eq__(self, other):
+        """Returns ``True`` if two domains are equivalent. """
+        return isinstance(other, FiniteField) and \
+            self.mod == other.mod and self.dom == other.dom
+
+    def characteristic(self):
+        """Return the characteristic of this domain. """
+        return self.mod
+
+    def get_field(self):
+        """Returns a field associated with ``self``. """
+        return self
+
+    def to_sympy(self, a):
+        """Convert ``a`` to a SymPy object. """
+        return SymPyInteger(self.to_int(a))
+
+    def from_sympy(self, a):
+        """Convert SymPy's Integer to SymPy's ``Integer``. """
+        if a.is_Integer:
+            return self.dtype(self.dom.dtype(int(a)))
+        elif int_valued(a):
+            return self.dtype(self.dom.dtype(int(a)))
+        else:
+            raise CoercionFailed("expected an integer, got %s" % a)
+
+    def to_int(self, a):
+        """Convert ``val`` to a Python ``int`` object. """
+        aval = int(a)
+        if self.sym and aval > self.mod // 2:
+            aval -= self.mod
+        return aval
+
+    def is_positive(self, a):
+        """Returns True if ``a`` is positive. """
+        return bool(a)
+
+    def is_nonnegative(self, a):
+        """Returns True if ``a`` is non-negative. """
+        return True
+
+    def is_negative(self, a):
+        """Returns True if ``a`` is negative. """
+        return False
+
+    def is_nonpositive(self, a):
+        """Returns True if ``a`` is non-positive. """
+        return not a
+
+    def from_FF(K1, a, K0=None):
+        """Convert ``ModularInteger(int)`` to ``dtype``. """
+        return K1.dtype(K1.dom.from_ZZ(int(a), K0.dom))
+
+    def from_FF_python(K1, a, K0=None):
+        """Convert ``ModularInteger(int)`` to ``dtype``. """
+        return K1.dtype(K1.dom.from_ZZ_python(int(a), K0.dom))
+
+    def from_ZZ(K1, a, K0=None):
+        """Convert Python's ``int`` to ``dtype``. """
+        return K1.dtype(K1.dom.from_ZZ_python(a, K0))
+
+    def from_ZZ_python(K1, a, K0=None):
+        """Convert Python's ``int`` to ``dtype``. """
+        return K1.dtype(K1.dom.from_ZZ_python(a, K0))
+
+    def from_QQ(K1, a, K0=None):
+        """Convert Python's ``Fraction`` to ``dtype``. """
+        if a.denominator == 1:
+            return K1.from_ZZ_python(a.numerator)
+
+    def from_QQ_python(K1, a, K0=None):
+        """Convert Python's ``Fraction`` to ``dtype``. """
+        if a.denominator == 1:
+            return K1.from_ZZ_python(a.numerator)
+
+    def from_FF_gmpy(K1, a, K0=None):
+        """Convert ``ModularInteger(mpz)`` to ``dtype``. """
+        return K1.dtype(K1.dom.from_ZZ_gmpy(a.val, K0.dom))
+
+    def from_ZZ_gmpy(K1, a, K0=None):
+        """Convert GMPY's ``mpz`` to ``dtype``. """
+        return K1.dtype(K1.dom.from_ZZ_gmpy(a, K0))
+
+    def from_QQ_gmpy(K1, a, K0=None):
+        """Convert GMPY's ``mpq`` to ``dtype``. """
+        if a.denominator == 1:
+            return K1.from_ZZ_gmpy(a.numerator)
+
+    def from_RealField(K1, a, K0):
+        """Convert mpmath's ``mpf`` to ``dtype``. """
+        p, q = K0.to_rational(a)
+
+        if q == 1:
+            return K1.dtype(K1.dom.dtype(p))
+
+    def is_square(self, a):
+        """Returns True if ``a`` is a quadratic residue modulo p. """
+        # a is not a square <=> x**2-a is irreducible
+        poly = [int(x) for x in [self.one, self.zero, -a]]
+        return not gf_irred_p_rabin(poly, self.mod, self.dom)
+
+    def exsqrt(self, a):
+        """Square root modulo p of ``a`` if it is a quadratic residue.
+
+        Explanation
+        ===========
+        Always returns the square root that is no larger than ``p // 2``.
+        """
+        # x**2-a is not square-free if a=0 or the field is characteristic 2
+        if self.mod == 2 or a == 0:
+            return a
+        # Otherwise, use square-free factorization routine to factorize x**2-a
+        poly = [int(x) for x in [self.one, self.zero, -a]]
+        for factor in gf_zassenhaus(poly, self.mod, self.dom):
+            if len(factor) == 2 and factor[1] <= self.mod // 2:
+                return self.dtype(factor[1])
+        return None
+
+
+FF = GF = FiniteField

lib/python3.12/site-packages/sympy/polys/domains/fractionfield.py

@@ -0,0 +1,177 @@
+"""Implementation of :class:`FractionField` class. """
+
+
+from sympy.polys.domains.compositedomain import CompositeDomain
+from sympy.polys.domains.field import Field
+from sympy.polys.polyerrors import CoercionFailed, GeneratorsError
+from sympy.utilities import public
+
+@public
+class FractionField(Field, CompositeDomain):
+    """A class for representing multivariate rational function fields. """
+
+    is_FractionField = is_Frac = True
+
+    has_assoc_Ring = True
+    has_assoc_Field = True
+
+    def __init__(self, domain_or_field, symbols=None, order=None):
+        from sympy.polys.fields import FracField
+
+        if isinstance(domain_or_field, FracField) and symbols is None and order is None:
+            field = domain_or_field
+        else:
+            field = FracField(symbols, domain_or_field, order)
+
+        self.field = field
+        self.dtype = field.dtype
+
+        self.gens = field.gens
+        self.ngens = field.ngens
+        self.symbols = field.symbols
+        self.domain = field.domain
+
+        # TODO: remove this
+        self.dom = self.domain
+
+    def new(self, element):
+        return self.field.field_new(element)
+
+    @property
+    def zero(self):
+        return self.field.zero
+
+    @property
+    def one(self):
+        return self.field.one
+
+    @property
+    def order(self):
+        return self.field.order
+
+    def __str__(self):
+        return str(self.domain) + '(' + ','.join(map(str, self.symbols)) + ')'
+
+    def __hash__(self):
+        return hash((self.__class__.__name__, self.dtype.field, self.domain, self.symbols))
+
+    def __eq__(self, other):
+        """Returns ``True`` if two domains are equivalent. """
+        return isinstance(other, FractionField) and \
+            (self.dtype.field, self.domain, self.symbols) ==\
+            (other.dtype.field, other.domain, other.symbols)
+
+    def to_sympy(self, a):
+        """Convert ``a`` to a SymPy object. """
+        return a.as_expr()
+
+    def from_sympy(self, a):
+        """Convert SymPy's expression to ``dtype``. """
+        return self.field.from_expr(a)
+
+    def from_ZZ(K1, a, K0):
+        """Convert a Python ``int`` object to ``dtype``. """
+        return K1(K1.domain.convert(a, K0))
+
+    def from_ZZ_python(K1, a, K0):
+        """Convert a Python ``int`` object to ``dtype``. """
+        return K1(K1.domain.convert(a, K0))
+
+    def from_QQ(K1, a, K0):
+        """Convert a Python ``Fraction`` object to ``dtype``. """
+        dom = K1.domain
+        conv = dom.convert_from
+        if dom.is_ZZ:
+            return K1(conv(K0.numer(a), K0)) / K1(conv(K0.denom(a), K0))
+        else:
+            return K1(conv(a, K0))
+
+    def from_QQ_python(K1, a, K0):
+        """Convert a Python ``Fraction`` object to ``dtype``. """
+        return K1(K1.domain.convert(a, K0))
+
+    def from_ZZ_gmpy(K1, a, K0):
+        """Convert a GMPY ``mpz`` object to ``dtype``. """
+        return K1(K1.domain.convert(a, K0))
+
+    def from_QQ_gmpy(K1, a, K0):
+        """Convert a GMPY ``mpq`` object to ``dtype``. """
+        return K1(K1.domain.convert(a, K0))
+
+    def from_GaussianRationalField(K1, a, K0):
+        """Convert a ``GaussianRational`` object to ``dtype``. """
+        return K1(K1.domain.convert(a, K0))
+
+    def from_GaussianIntegerRing(K1, a, K0):
+        """Convert a ``GaussianInteger`` object to ``dtype``. """
+        return K1(K1.domain.convert(a, K0))
+
+    def from_RealField(K1, a, K0):
+        """Convert a mpmath ``mpf`` object to ``dtype``. """
+        return K1(K1.domain.convert(a, K0))
+
+    def from_ComplexField(K1, a, K0):
+        """Convert a mpmath ``mpf`` object to ``dtype``. """
+        return K1(K1.domain.convert(a, K0))
+
+    def from_AlgebraicField(K1, a, K0):
+        """Convert an algebraic number to ``dtype``. """
+        if K1.domain != K0:
+            a = K1.domain.convert_from(a, K0)
+        if a is not None:
+            return K1.new(a)
+
+    def from_PolynomialRing(K1, a, K0):
+        """Convert a polynomial to ``dtype``. """
+        if a.is_ground:
+            return K1.convert_from(a.coeff(1), K0.domain)
+        try:
+            return K1.new(a.set_ring(K1.field.ring))
+        except (CoercionFailed, GeneratorsError):
+            # XXX: We get here if K1=ZZ(x,y) and K0=QQ[x,y]
+            # and the poly a in K0 has non-integer coefficients.
+            # It seems that K1.new can handle this but K1.new doesn't work
+            # when K0.domain is an algebraic field...
+            try:
+                return K1.new(a)
+            except (CoercionFailed, GeneratorsError):
+                return None
+
+    def from_FractionField(K1, a, K0):
+        """Convert a rational function to ``dtype``. """
+        try:
+            return a.set_field(K1.field)
+        except (CoercionFailed, GeneratorsError):
+            return None
+
+    def get_ring(self):
+        """Returns a field associated with ``self``. """
+        return self.field.to_ring().to_domain()
+
+    def is_positive(self, a):
+        """Returns True if ``LC(a)`` is positive. """
+        return self.domain.is_positive(a.numer.LC)
+
+    def is_negative(self, a):
+        """Returns True if ``LC(a)`` is negative. """
+        return self.domain.is_negative(a.numer.LC)
+
+    def is_nonpositive(self, a):
+        """Returns True if ``LC(a)`` is non-positive. """
+        return self.domain.is_nonpositive(a.numer.LC)
+
+    def is_nonnegative(self, a):
+        """Returns True if ``LC(a)`` is non-negative. """
+        return self.domain.is_nonnegative(a.numer.LC)
+
+    def numer(self, a):
+        """Returns numerator of ``a``. """
+        return a.numer
+
+    def denom(self, a):
+        """Returns denominator of ``a``. """
+        return a.denom
+
+    def factorial(self, a):
+        """Returns factorial of ``a``. """
+        return self.dtype(self.domain.factorial(a))

lib/python3.12/site-packages/sympy/polys/domains/gaussiandomains.py

@@ -0,0 +1,686 @@
+"""Domains of Gaussian type."""
+
+from sympy.core.numbers import I
+from sympy.polys.polyerrors import CoercionFailed
+from sympy.polys.domains.integerring import ZZ
+from sympy.polys.domains.rationalfield import QQ
+from sympy.polys.domains.algebraicfield import AlgebraicField
+from sympy.polys.domains.domain import Domain
+from sympy.polys.domains.domainelement import DomainElement
+from sympy.polys.domains.field import Field
+from sympy.polys.domains.ring import Ring
+
+
+class GaussianElement(DomainElement):
+    """Base class for elements of Gaussian type domains."""
+    base: Domain
+    _parent: Domain
+
+    __slots__ = ('x', 'y')
+
+    def __new__(cls, x, y=0):
+        conv = cls.base.convert
+        return cls.new(conv(x), conv(y))
+
+    @classmethod
+    def new(cls, x, y):
+        """Create a new GaussianElement of the same domain."""
+        obj = super().__new__(cls)
+        obj.x = x
+        obj.y = y
+        return obj
+
+    def parent(self):
+        """The domain that this is an element of (ZZ_I or QQ_I)"""
+        return self._parent
+
+    def __hash__(self):
+        return hash((self.x, self.y))
+
+    def __eq__(self, other):
+        if isinstance(other, self.__class__):
+            return self.x == other.x and self.y == other.y
+        else:
+            return NotImplemented
+
+    def __lt__(self, other):
+        if not isinstance(other, GaussianElement):
+            return NotImplemented
+        return [self.y, self.x] < [other.y, other.x]
+
+    def __pos__(self):
+        return self
+
+    def __neg__(self):
+        return self.new(-self.x, -self.y)
+
+    def __repr__(self):
+        return "%s(%s, %s)" % (self._parent.rep, self.x, self.y)
+
+    def __str__(self):
+        return str(self._parent.to_sympy(self))
+
+    @classmethod
+    def _get_xy(cls, other):
+        if not isinstance(other, cls):
+            try:
+                other = cls._parent.convert(other)
+            except CoercionFailed:
+                return None, None
+        return other.x, other.y
+
+    def __add__(self, other):
+        x, y = self._get_xy(other)
+        if x is not None:
+            return self.new(self.x + x, self.y + y)
+        else:
+            return NotImplemented
+
+    __radd__ = __add__
+
+    def __sub__(self, other):
+        x, y = self._get_xy(other)
+        if x is not None:
+            return self.new(self.x - x, self.y - y)
+        else:
+            return NotImplemented
+
+    def __rsub__(self, other):
+        x, y = self._get_xy(other)
+        if x is not None:
+            return self.new(x - self.x, y - self.y)
+        else:
+            return NotImplemented
+
+    def __mul__(self, other):
+        x, y = self._get_xy(other)
+        if x is not None:
+            return self.new(self.x*x - self.y*y, self.x*y + self.y*x)
+        else:
+            return NotImplemented
+
+    __rmul__ = __mul__
+
+    def __pow__(self, exp):
+        if exp == 0:
+            return self.new(1, 0)
+        if exp < 0:
+            self, exp = 1/self, -exp
+        if exp == 1:
+            return self
+        pow2 = self
+        prod = self if exp % 2 else self._parent.one
+        exp //= 2
+        while exp:
+            pow2 *= pow2
+            if exp % 2:
+                prod *= pow2
+            exp //= 2
+        return prod
+
+    def __bool__(self):
+        return bool(self.x) or bool(self.y)
+
+    def quadrant(self):
+        """Return quadrant index 0-3.
+
+        0 is included in quadrant 0.
+        """
+        if self.y > 0:
+            return 0 if self.x > 0 else 1
+        elif self.y < 0:
+            return 2 if self.x < 0 else 3
+        else:
+            return 0 if self.x >= 0 else 2
+
+    def __rdivmod__(self, other):
+        try:
+            other = self._parent.convert(other)
+        except CoercionFailed:
+            return NotImplemented
+        else:
+            return other.__divmod__(self)
+
+    def __rtruediv__(self, other):
+        try:
+            other = QQ_I.convert(other)
+        except CoercionFailed:
+            return NotImplemented
+        else:
+            return other.__truediv__(self)
+
+    def __floordiv__(self, other):
+        qr = self.__divmod__(other)
+        return qr if qr is NotImplemented else qr[0]
+
+    def __rfloordiv__(self, other):
+        qr = self.__rdivmod__(other)
+        return qr if qr is NotImplemented else qr[0]
+
+    def __mod__(self, other):
+        qr = self.__divmod__(other)
+        return qr if qr is NotImplemented else qr[1]
+
+    def __rmod__(self, other):
+        qr = self.__rdivmod__(other)
+        return qr if qr is NotImplemented else qr[1]
+
+
+class GaussianInteger(GaussianElement):
+    """Gaussian integer: domain element for :ref:`ZZ_I`
+
+        >>> from sympy import ZZ_I
+        >>> z = ZZ_I(2, 3)
+        >>> z
+        (2 + 3*I)
+        >>> type(z)
+        <class 'sympy.polys.domains.gaussiandomains.GaussianInteger'>
+    """
+    base = ZZ
+
+    def __truediv__(self, other):
+        """Return a Gaussian rational."""
+        return QQ_I.convert(self)/other
+
+    def __divmod__(self, other):
+        if not other:
+            raise ZeroDivisionError('divmod({}, 0)'.format(self))
+        x, y = self._get_xy(other)
+        if x is None:
+            return NotImplemented
+
+        # multiply self and other by x - I*y
+        # self/other == (a + I*b)/c
+        a, b = self.x*x + self.y*y, -self.x*y + self.y*x
+        c = x*x + y*y
+
+        # find integers qx and qy such that
+        # |a - qx*c| <= c/2 and |b - qy*c| <= c/2
+        qx = (2*a + c) // (2*c)  # -c <= 2*a - qx*2*c < c
+        qy = (2*b + c) // (2*c)
+
+        q = GaussianInteger(qx, qy)
+        # |self/other - q| < 1 since
+        # |a/c - qx|**2 + |b/c - qy|**2 <= 1/4 + 1/4 < 1
+
+        return q, self - q*other  # |r| < |other|
+
+
+class GaussianRational(GaussianElement):
+    """Gaussian rational: domain element for :ref:`QQ_I`
+
+        >>> from sympy import QQ_I, QQ
+        >>> z = QQ_I(QQ(2, 3), QQ(4, 5))
+        >>> z
+        (2/3 + 4/5*I)
+        >>> type(z)
+        <class 'sympy.polys.domains.gaussiandomains.GaussianRational'>
+    """
+    base = QQ
+
+    def __truediv__(self, other):
+        """Return a Gaussian rational."""
+        if not other:
+            raise ZeroDivisionError('{} / 0'.format(self))
+        x, y = self._get_xy(other)
+        if x is None:
+            return NotImplemented
+        c = x*x + y*y
+
+        return GaussianRational((self.x*x + self.y*y)/c,
+                                (-self.x*y + self.y*x)/c)
+
+    def __divmod__(self, other):
+        try:
+            other = self._parent.convert(other)
+        except CoercionFailed:
+            return NotImplemented
+        if not other:
+            raise ZeroDivisionError('{} % 0'.format(self))
+        else:
+            return self/other, QQ_I.zero
+
+
+class GaussianDomain():
+    """Base class for Gaussian domains."""
+    dom = None  # type: Domain
+
+    is_Numerical = True
+    is_Exact = True
+
+    has_assoc_Ring = True
+    has_assoc_Field = True
+
+    def to_sympy(self, a):
+        """Convert ``a`` to a SymPy object. """
+        conv = self.dom.to_sympy
+        return conv(a.x) + I*conv(a.y)
+
+    def from_sympy(self, a):
+        """Convert a SymPy object to ``self.dtype``."""
+        r, b = a.as_coeff_Add()
+        x = self.dom.from_sympy(r)  # may raise CoercionFailed
+        if not b:
+            return self.new(x, 0)
+        r, b = b.as_coeff_Mul()
+        y = self.dom.from_sympy(r)
+        if b is I:
+            return self.new(x, y)
+        else:
+            raise CoercionFailed("{} is not Gaussian".format(a))
+
+    def inject(self, *gens):
+        """Inject generators into this domain. """
+        return self.poly_ring(*gens)
+
+    def canonical_unit(self, d):
+        unit = self.units[-d.quadrant()]  # - for inverse power
+        return unit
+
+    def is_negative(self, element):
+        """Returns ``False`` for any ``GaussianElement``. """
+        return False
+
+    def is_positive(self, element):
+        """Returns ``False`` for any ``GaussianElement``. """
+        return False
+
+    def is_nonnegative(self, element):
+        """Returns ``False`` for any ``GaussianElement``. """
+        return False
+
+    def is_nonpositive(self, element):
+        """Returns ``False`` for any ``GaussianElement``. """
+        return False
+
+    def from_ZZ_gmpy(K1, a, K0):
+        """Convert a GMPY mpz to ``self.dtype``."""
+        return K1(a)
+
+    def from_ZZ(K1, a, K0):
+        """Convert a ZZ_python element to ``self.dtype``."""
+        return K1(a)
+
+    def from_ZZ_python(K1, a, K0):
+        """Convert a ZZ_python element to ``self.dtype``."""
+        return K1(a)
+
+    def from_QQ(K1, a, K0):
+        """Convert a GMPY mpq to ``self.dtype``."""
+        return K1(a)
+
+    def from_QQ_gmpy(K1, a, K0):
+        """Convert a GMPY mpq to ``self.dtype``."""
+        return K1(a)
+
+    def from_QQ_python(K1, a, K0):
+        """Convert a QQ_python element to ``self.dtype``."""
+        return K1(a)
+
+    def from_AlgebraicField(K1, a, K0):
+        """Convert an element from ZZ<I> or QQ<I> to ``self.dtype``."""
+        if K0.ext.args[0] == I:
+            return K1.from_sympy(K0.to_sympy(a))
+
+
+class GaussianIntegerRing(GaussianDomain, Ring):
+    r"""Ring of Gaussian integers ``ZZ_I``
+
+    The :ref:`ZZ_I` domain represents the `Gaussian integers`_ `\mathbb{Z}[i]`
+    as a :py:class:`~.Domain` in the domain system (see
+    :ref:`polys-domainsintro`).
+
+    By default a :py:class:`~.Poly` created from an expression with
+    coefficients that are combinations of integers and ``I`` (`\sqrt{-1}`)
+    will have the domain :ref:`ZZ_I`.
+
+    >>> from sympy import Poly, Symbol, I
+    >>> x = Symbol('x')
+    >>> p = Poly(x**2 + I)
+    >>> p
+    Poly(x**2 + I, x, domain='ZZ_I')
+    >>> p.domain
+    ZZ_I
+
+    The :ref:`ZZ_I` domain can be used to factorise polynomials that are
+    reducible over the Gaussian integers.
+
+    >>> from sympy import factor
+    >>> factor(x**2 + 1)
+    x**2 + 1
+    >>> factor(x**2 + 1, domain='ZZ_I')
+    (x - I)*(x + I)
+
+    The corresponding `field of fractions`_ is the domain of the Gaussian
+    rationals :ref:`QQ_I`. Conversely :ref:`ZZ_I` is the `ring of integers`_
+    of :ref:`QQ_I`.
+
+    >>> from sympy import ZZ_I, QQ_I
+    >>> ZZ_I.get_field()
+    QQ_I
+    >>> QQ_I.get_ring()
+    ZZ_I
+
+    When using the domain directly :ref:`ZZ_I` can be used as a constructor.
+
+    >>> ZZ_I(3, 4)
+    (3 + 4*I)
+    >>> ZZ_I(5)
+    (5 + 0*I)
+
+    The domain elements of :ref:`ZZ_I` are instances of
+    :py:class:`~.GaussianInteger` which support the rings operations
+    ``+,-,*,**``.
+
+    >>> z1 = ZZ_I(5, 1)
+    >>> z2 = ZZ_I(2, 3)
+    >>> z1
+    (5 + 1*I)
+    >>> z2
+    (2 + 3*I)
+    >>> z1 + z2
+    (7 + 4*I)
+    >>> z1 * z2
+    (7 + 17*I)
+    >>> z1 ** 2
+    (24 + 10*I)
+
+    Both floor (``//``) and modulo (``%``) division work with
+    :py:class:`~.GaussianInteger` (see the :py:meth:`~.Domain.div` method).
+
+    >>> z3, z4 = ZZ_I(5), ZZ_I(1, 3)
+    >>> z3 // z4  # floor division
+    (1 + -1*I)
+    >>> z3 % z4   # modulo division (remainder)
+    (1 + -2*I)
+    >>> (z3//z4)*z4 + z3%z4 == z3
+    True
+
+    True division (``/``) in :ref:`ZZ_I` gives an element of :ref:`QQ_I`. The
+    :py:meth:`~.Domain.exquo` method can be used to divide in :ref:`ZZ_I` when
+    exact division is possible.
+
+    >>> z1 / z2
+    (1 + -1*I)
+    >>> ZZ_I.exquo(z1, z2)
+    (1 + -1*I)
+    >>> z3 / z4
+    (1/2 + -3/2*I)
+    >>> ZZ_I.exquo(z3, z4)
+    Traceback (most recent call last):
+        ...
+    ExactQuotientFailed: (1 + 3*I) does not divide (5 + 0*I) in ZZ_I
+
+    The :py:meth:`~.Domain.gcd` method can be used to compute the `gcd`_ of any
+    two elements.
+
+    >>> ZZ_I.gcd(ZZ_I(10), ZZ_I(2))
+    (2 + 0*I)
+    >>> ZZ_I.gcd(ZZ_I(5), ZZ_I(2, 1))
+    (2 + 1*I)
+
+    .. _Gaussian integers: https://en.wikipedia.org/wiki/Gaussian_integer
+    .. _gcd: https://en.wikipedia.org/wiki/Greatest_common_divisor
+
+    """
+    dom = ZZ
+    dtype = GaussianInteger
+    zero = dtype(ZZ(0), ZZ(0))
+    one = dtype(ZZ(1), ZZ(0))
+    imag_unit = dtype(ZZ(0), ZZ(1))
+    units = (one, imag_unit, -one, -imag_unit)  # powers of i
+
+    rep = 'ZZ_I'
+
+    is_GaussianRing = True
+    is_ZZ_I = True
+
+    def __init__(self):  # override Domain.__init__
+        """For constructing ZZ_I."""
+
+    def __eq__(self, other):
+        """Returns ``True`` if two domains are equivalent. """
+        if isinstance(other, GaussianIntegerRing):
+            return True
+        else:
+            return NotImplemented
+
+    def __hash__(self):
+        """Compute hash code of ``self``. """
+        return hash('ZZ_I')
+
+    @property
+    def has_CharacteristicZero(self):
+        return True
+
+    def characteristic(self):
+        return 0
+
+    def get_ring(self):
+        """Returns a ring associated with ``self``. """
+        return self
+
+    def get_field(self):
+        """Returns a field associated with ``self``. """
+        return QQ_I
+
+    def normalize(self, d, *args):
+        """Return first quadrant element associated with ``d``.
+
+        Also multiply the other arguments by the same power of i.
+        """
+        unit = self.canonical_unit(d)
+        d *= unit
+        args = tuple(a*unit for a in args)
+        return (d,) + args if args else d
+
+    def gcd(self, a, b):
+        """Greatest common divisor of a and b over ZZ_I."""
+        while b:
+            a, b = b, a % b
+        return self.normalize(a)
+
+    def lcm(self, a, b):
+        """Least common multiple of a and b over ZZ_I."""
+        return (a * b) // self.gcd(a, b)
+
+    def from_GaussianIntegerRing(K1, a, K0):
+        """Convert a ZZ_I element to ZZ_I."""
+        return a
+
+    def from_GaussianRationalField(K1, a, K0):
+        """Convert a QQ_I element to ZZ_I."""
+        return K1.new(ZZ.convert(a.x), ZZ.convert(a.y))
+
+ZZ_I = GaussianInteger._parent = GaussianIntegerRing()
+
+
+class GaussianRationalField(GaussianDomain, Field):
+    r"""Field of Gaussian rationals ``QQ_I``
+
+    The :ref:`QQ_I` domain represents the `Gaussian rationals`_ `\mathbb{Q}(i)`
+    as a :py:class:`~.Domain` in the domain system (see
+    :ref:`polys-domainsintro`).
+
+    By default a :py:class:`~.Poly` created from an expression with
+    coefficients that are combinations of rationals and ``I`` (`\sqrt{-1}`)
+    will have the domain :ref:`QQ_I`.
+
+    >>> from sympy import Poly, Symbol, I
+    >>> x = Symbol('x')
+    >>> p = Poly(x**2 + I/2)
+    >>> p
+    Poly(x**2 + I/2, x, domain='QQ_I')
+    >>> p.domain
+    QQ_I
+
+    The polys option ``gaussian=True`` can be used to specify that the domain
+    should be :ref:`QQ_I` even if the coefficients do not contain ``I`` or are
+    all integers.
+
+    >>> Poly(x**2)
+    Poly(x**2, x, domain='ZZ')
+    >>> Poly(x**2 + I)
+    Poly(x**2 + I, x, domain='ZZ_I')
+    >>> Poly(x**2/2)
+    Poly(1/2*x**2, x, domain='QQ')
+    >>> Poly(x**2, gaussian=True)
+    Poly(x**2, x, domain='QQ_I')
+    >>> Poly(x**2 + I, gaussian=True)
+    Poly(x**2 + I, x, domain='QQ_I')
+    >>> Poly(x**2/2, gaussian=True)
+    Poly(1/2*x**2, x, domain='QQ_I')
+
+    The :ref:`QQ_I` domain can be used to factorise polynomials that are
+    reducible over the Gaussian rationals.
+
+    >>> from sympy import factor, QQ_I
+    >>> factor(x**2/4 + 1)
+    (x**2 + 4)/4
+    >>> factor(x**2/4 + 1, domain='QQ_I')
+    (x - 2*I)*(x + 2*I)/4
+    >>> factor(x**2/4 + 1, domain=QQ_I)
+    (x - 2*I)*(x + 2*I)/4
+
+    It is also possible to specify the :ref:`QQ_I` domain explicitly with
+    polys functions like :py:func:`~.apart`.
+
+    >>> from sympy import apart
+    >>> apart(1/(1 + x**2))
+    1/(x**2 + 1)
+    >>> apart(1/(1 + x**2), domain=QQ_I)
+    I/(2*(x + I)) - I/(2*(x - I))
+
+    The corresponding `ring of integers`_ is the domain of the Gaussian
+    integers :ref:`ZZ_I`. Conversely :ref:`QQ_I` is the `field of fractions`_
+    of :ref:`ZZ_I`.
+
+    >>> from sympy import ZZ_I, QQ_I, QQ
+    >>> ZZ_I.get_field()
+    QQ_I
+    >>> QQ_I.get_ring()
+    ZZ_I
+
+    When using the domain directly :ref:`QQ_I` can be used as a constructor.
+
+    >>> QQ_I(3, 4)
+    (3 + 4*I)
+    >>> QQ_I(5)
+    (5 + 0*I)
+    >>> QQ_I(QQ(2, 3), QQ(4, 5))
+    (2/3 + 4/5*I)
+
+    The domain elements of :ref:`QQ_I` are instances of
+    :py:class:`~.GaussianRational` which support the field operations
+    ``+,-,*,**,/``.
+
+    >>> z1 = QQ_I(5, 1)
+    >>> z2 = QQ_I(2, QQ(1, 2))
+    >>> z1
+    (5 + 1*I)
+    >>> z2
+    (2 + 1/2*I)
+    >>> z1 + z2
+    (7 + 3/2*I)
+    >>> z1 * z2
+    (19/2 + 9/2*I)
+    >>> z2 ** 2
+    (15/4 + 2*I)
+
+    True division (``/``) in :ref:`QQ_I` gives an element of :ref:`QQ_I` and
+    is always exact.
+
+    >>> z1 / z2
+    (42/17 + -2/17*I)
+    >>> QQ_I.exquo(z1, z2)
+    (42/17 + -2/17*I)
+    >>> z1 == (z1/z2)*z2
+    True
+
+    Both floor (``//``) and modulo (``%``) division can be used with
+    :py:class:`~.GaussianRational` (see :py:meth:`~.Domain.div`)
+    but division is always exact so there is no remainder.
+
+    >>> z1 // z2
+    (42/17 + -2/17*I)
+    >>> z1 % z2
+    (0 + 0*I)
+    >>> QQ_I.div(z1, z2)
+    ((42/17 + -2/17*I), (0 + 0*I))
+    >>> (z1//z2)*z2 + z1%z2 == z1
+    True
+
+    .. _Gaussian rationals: https://en.wikipedia.org/wiki/Gaussian_rational
+    """
+    dom = QQ
+    dtype = GaussianRational
+    zero = dtype(QQ(0), QQ(0))
+    one = dtype(QQ(1), QQ(0))
+    imag_unit = dtype(QQ(0), QQ(1))
+    units = (one, imag_unit, -one, -imag_unit)  # powers of i
+
+    rep = 'QQ_I'
+
+    is_GaussianField = True
+    is_QQ_I = True
+
+    def __init__(self):  # override Domain.__init__
+        """For constructing QQ_I."""
+
+    def __eq__(self, other):
+        """Returns ``True`` if two domains are equivalent. """
+        if isinstance(other, GaussianRationalField):
+            return True
+        else:
+            return NotImplemented
+
+    def __hash__(self):
+        """Compute hash code of ``self``. """
+        return hash('QQ_I')
+
+    @property
+    def has_CharacteristicZero(self):
+        return True
+
+    def characteristic(self):
+        return 0
+
+    def get_ring(self):
+        """Returns a ring associated with ``self``. """
+        return ZZ_I
+
+    def get_field(self):
+        """Returns a field associated with ``self``. """
+        return self
+
+    def as_AlgebraicField(self):
+        """Get equivalent domain as an ``AlgebraicField``. """
+        return AlgebraicField(self.dom, I)
+
+    def numer(self, a):
+        """Get the numerator of ``a``."""
+        ZZ_I = self.get_ring()
+        return ZZ_I.convert(a * self.denom(a))
+
+    def denom(self, a):
+        """Get the denominator of ``a``."""
+        ZZ = self.dom.get_ring()
+        QQ = self.dom
+        ZZ_I = self.get_ring()
+        denom_ZZ = ZZ.lcm(QQ.denom(a.x), QQ.denom(a.y))
+        return ZZ_I(denom_ZZ, ZZ.zero)
+
+    def from_GaussianIntegerRing(K1, a, K0):
+        """Convert a ZZ_I element to QQ_I."""
+        return K1.new(a.x, a.y)
+
+    def from_GaussianRationalField(K1, a, K0):
+        """Convert a QQ_I element to QQ_I."""
+        return a
+
+    def from_ComplexField(K1, a, K0):
+        """Convert a ComplexField element to QQ_I."""
+        return K1.new(QQ.convert(a.real), QQ.convert(a.imag))
+
+
+QQ_I = GaussianRational._parent = GaussianRationalField()

lib/python3.12/site-packages/sympy/polys/domains/gmpyfinitefield.py

@@ -0,0 +1,16 @@
+"""Implementation of :class:`GMPYFiniteField` class. """
+
+
+from sympy.polys.domains.finitefield import FiniteField
+from sympy.polys.domains.gmpyintegerring import GMPYIntegerRing
+
+from sympy.utilities import public
+
+@public
+class GMPYFiniteField(FiniteField):
+    """Finite field based on GMPY integers. """
+
+    alias = 'FF_gmpy'
+
+    def __init__(self, mod, symmetric=True):
+        super().__init__(mod, GMPYIntegerRing(), symmetric)

lib/python3.12/site-packages/sympy/polys/domains/gmpyintegerring.py

@@ -0,0 +1,105 @@
+"""Implementation of :class:`GMPYIntegerRing` class. """
+
+
+from sympy.polys.domains.groundtypes import (
+    GMPYInteger, SymPyInteger,
+    factorial as gmpy_factorial,
+    gmpy_gcdex, gmpy_gcd, gmpy_lcm, sqrt as gmpy_sqrt,
+)
+from sympy.core.numbers import int_valued
+from sympy.polys.domains.integerring import IntegerRing
+from sympy.polys.polyerrors import CoercionFailed
+from sympy.utilities import public
+
+@public
+class GMPYIntegerRing(IntegerRing):
+    """Integer ring based on GMPY's ``mpz`` type.
+
+    This will be the implementation of :ref:`ZZ` if ``gmpy`` or ``gmpy2`` is
+    installed. Elements will be of type ``gmpy.mpz``.
+    """
+
+    dtype = GMPYInteger
+    zero = dtype(0)
+    one = dtype(1)
+    tp = type(one)
+    alias = 'ZZ_gmpy'
+
+    def __init__(self):
+        """Allow instantiation of this domain. """
+
+    def to_sympy(self, a):
+        """Convert ``a`` to a SymPy object. """
+        return SymPyInteger(int(a))
+
+    def from_sympy(self, a):
+        """Convert SymPy's Integer to ``dtype``. """
+        if a.is_Integer:
+            return GMPYInteger(a.p)
+        elif int_valued(a):
+            return GMPYInteger(int(a))
+        else:
+            raise CoercionFailed("expected an integer, got %s" % a)
+
+    def from_FF_python(K1, a, K0):
+        """Convert ``ModularInteger(int)`` to GMPY's ``mpz``. """
+        return K0.to_int(a)
+
+    def from_ZZ_python(K1, a, K0):
+        """Convert Python's ``int`` to GMPY's ``mpz``. """
+        return GMPYInteger(a)
+
+    def from_QQ(K1, a, K0):
+        """Convert Python's ``Fraction`` to GMPY's ``mpz``. """
+        if a.denominator == 1:
+            return GMPYInteger(a.numerator)
+
+    def from_QQ_python(K1, a, K0):
+        """Convert Python's ``Fraction`` to GMPY's ``mpz``. """
+        if a.denominator == 1:
+            return GMPYInteger(a.numerator)
+
+    def from_FF_gmpy(K1, a, K0):
+        """Convert ``ModularInteger(mpz)`` to GMPY's ``mpz``. """
+        return K0.to_int(a)
+
+    def from_ZZ_gmpy(K1, a, K0):
+        """Convert GMPY's ``mpz`` to GMPY's ``mpz``. """
+        return a
+
+    def from_QQ_gmpy(K1, a, K0):
+        """Convert GMPY ``mpq`` to GMPY's ``mpz``. """
+        if a.denominator == 1:
+            return a.numerator
+
+    def from_RealField(K1, a, K0):
+        """Convert mpmath's ``mpf`` to GMPY's ``mpz``. """
+        p, q = K0.to_rational(a)
+
+        if q == 1:
+            return GMPYInteger(p)
+
+    def from_GaussianIntegerRing(K1, a, K0):
+        if a.y == 0:
+            return a.x
+
+    def gcdex(self, a, b):
+        """Compute extended GCD of ``a`` and ``b``. """
+        h, s, t = gmpy_gcdex(a, b)
+        return s, t, h
+
+    def gcd(self, a, b):
+        """Compute GCD of ``a`` and ``b``. """
+        return gmpy_gcd(a, b)
+
+    def lcm(self, a, b):
+        """Compute LCM of ``a`` and ``b``. """
+        return gmpy_lcm(a, b)
+
+    def sqrt(self, a):
+        """Compute square root of ``a``. """
+        return gmpy_sqrt(a)
+
+    def factorial(self, a):
+        """Compute factorial of ``a``. """
+        return gmpy_factorial(a)

lib/python3.12/site-packages/sympy/polys/domains/gmpyrationalfield.py

@@ -0,0 +1,100 @@
+"""Implementation of :class:`GMPYRationalField` class. """
+
+
+from sympy.polys.domains.groundtypes import (
+    GMPYRational, SymPyRational,
+    gmpy_numer, gmpy_denom, factorial as gmpy_factorial,
+)
+from sympy.polys.domains.rationalfield import RationalField
+from sympy.polys.polyerrors import CoercionFailed
+from sympy.utilities import public
+
+@public
+class GMPYRationalField(RationalField):
+    """Rational field based on GMPY's ``mpq`` type.
+
+    This will be the implementation of :ref:`QQ` if ``gmpy`` or ``gmpy2`` is
+    installed. Elements will be of type ``gmpy.mpq``.
+    """
+
+    dtype = GMPYRational
+    zero = dtype(0)
+    one = dtype(1)
+    tp = type(one)
+    alias = 'QQ_gmpy'
+
+    def __init__(self):
+        pass
+
+    def get_ring(self):
+        """Returns ring associated with ``self``. """
+        from sympy.polys.domains import GMPYIntegerRing
+        return GMPYIntegerRing()
+
+    def to_sympy(self, a):
+        """Convert ``a`` to a SymPy object. """
+        return SymPyRational(int(gmpy_numer(a)),
+                             int(gmpy_denom(a)))
+
+    def from_sympy(self, a):
+        """Convert SymPy's Integer to ``dtype``. """
+        if a.is_Rational:
+            return GMPYRational(a.p, a.q)
+        elif a.is_Float:
+            from sympy.polys.domains import RR
+            return GMPYRational(*map(int, RR.to_rational(a)))
+        else:
+            raise CoercionFailed("expected ``Rational`` object, got %s" % a)
+
+    def from_ZZ_python(K1, a, K0):
+        """Convert a Python ``int`` object to ``dtype``. """
+        return GMPYRational(a)
+
+    def from_QQ_python(K1, a, K0):
+        """Convert a Python ``Fraction`` object to ``dtype``. """
+        return GMPYRational(a.numerator, a.denominator)
+
+    def from_ZZ_gmpy(K1, a, K0):
+        """Convert a GMPY ``mpz`` object to ``dtype``. """
+        return GMPYRational(a)
+
+    def from_QQ_gmpy(K1, a, K0):
+        """Convert a GMPY ``mpq`` object to ``dtype``. """
+        return a
+
+    def from_GaussianRationalField(K1, a, K0):
+        """Convert a ``GaussianElement`` object to ``dtype``. """
+        if a.y == 0:
+            return GMPYRational(a.x)
+
+    def from_RealField(K1, a, K0):
+        """Convert a mpmath ``mpf`` object to ``dtype``. """
+        return GMPYRational(*map(int, K0.to_rational(a)))
+
+    def exquo(self, a, b):
+        """Exact quotient of ``a`` and ``b``, implies ``__truediv__``.  """
+        return GMPYRational(a) / GMPYRational(b)
+
+    def quo(self, a, b):
+        """Quotient of ``a`` and ``b``, implies ``__truediv__``. """
+        return GMPYRational(a) / GMPYRational(b)
+
+    def rem(self, a, b):
+        """Remainder of ``a`` and ``b``, implies nothing.  """
+        return self.zero
+
+    def div(self, a, b):
+        """Division of ``a`` and ``b``, implies ``__truediv__``. """
+        return GMPYRational(a) / GMPYRational(b), self.zero
+
+    def numer(self, a):
+        """Returns numerator of ``a``. """
+        return a.numerator
+
+    def denom(self, a):
+        """Returns denominator of ``a``. """
+        return a.denominator
+
+    def factorial(self, a):
+        """Returns factorial of ``a``. """
+        return GMPYRational(gmpy_factorial(int(a)))

lib/python3.12/site-packages/sympy/polys/domains/groundtypes.py

@@ -0,0 +1,99 @@
+"""Ground types for various mathematical domains in SymPy. """
+
+import builtins
+from sympy.external.gmpy import GROUND_TYPES, factorial, sqrt, is_square, sqrtrem
+
+PythonInteger = builtins.int
+PythonReal = builtins.float
+PythonComplex = builtins.complex
+
+from .pythonrational import PythonRational
+
+from sympy.core.intfunc import (
+    igcdex as python_gcdex,
+    igcd2 as python_gcd,
+    ilcm as python_lcm,
+)
+
+from sympy.core.numbers import (Float as SymPyReal, Integer as SymPyInteger, Rational as SymPyRational)
+
+
+class _GMPYInteger:
+    def __init__(self, obj):
+        pass
+
+class _GMPYRational:
+    def __init__(self, obj):
+        pass
+
+
+if GROUND_TYPES == 'gmpy':
+
+    from gmpy2 import (
+        mpz as GMPYInteger,
+        mpq as GMPYRational,
+        numer as gmpy_numer,
+        denom as gmpy_denom,
+        gcdext as gmpy_gcdex,
+        gcd as gmpy_gcd,
+        lcm as gmpy_lcm,
+        qdiv as gmpy_qdiv,
+    )
+    gcdex = gmpy_gcdex
+    gcd = gmpy_gcd
+    lcm = gmpy_lcm
+
+elif GROUND_TYPES == 'flint':
+
+    from flint import fmpz as _fmpz
+
+    GMPYInteger = _GMPYInteger
+    GMPYRational = _GMPYRational
+    gmpy_numer = None
+    gmpy_denom = None
+    gmpy_gcdex = None
+    gmpy_gcd = None
+    gmpy_lcm = None
+    gmpy_qdiv = None
+
+    def gcd(a, b):
+        return a.gcd(b)
+
+    def gcdex(a, b):
+        x, y, g = python_gcdex(a, b)
+        return _fmpz(x), _fmpz(y), _fmpz(g)
+
+    def lcm(a, b):
+        return a.lcm(b)
+
+else:
+    GMPYInteger = _GMPYInteger
+    GMPYRational = _GMPYRational
+    gmpy_numer = None
+    gmpy_denom = None
+    gmpy_gcdex = None
+    gmpy_gcd = None
+    gmpy_lcm = None
+    gmpy_qdiv = None
+    gcdex = python_gcdex
+    gcd = python_gcd
+    lcm = python_lcm
+
+
+__all__ = [
+    'PythonInteger', 'PythonReal', 'PythonComplex',
+
+    'PythonRational',
+
+    'python_gcdex', 'python_gcd', 'python_lcm',
+
+    'SymPyReal', 'SymPyInteger', 'SymPyRational',
+
+    'GMPYInteger', 'GMPYRational', 'gmpy_numer',
+    'gmpy_denom', 'gmpy_gcdex', 'gmpy_gcd', 'gmpy_lcm',
+    'gmpy_qdiv',
+
+    'factorial', 'sqrt', 'is_square', 'sqrtrem',
+
+    'GMPYInteger', 'GMPYRational',
+]

lib/python3.12/site-packages/sympy/polys/domains/integerring.py

@@ -0,0 +1,276 @@
+"""Implementation of :class:`IntegerRing` class. """
+
+from sympy.external.gmpy import MPZ, GROUND_TYPES
+
+from sympy.core.numbers import int_valued
+from sympy.polys.domains.groundtypes import (
+    SymPyInteger,
+    factorial,
+    gcdex, gcd, lcm, sqrt, is_square, sqrtrem,
+)
+
+from sympy.polys.domains.characteristiczero import CharacteristicZero
+from sympy.polys.domains.ring import Ring
+from sympy.polys.domains.simpledomain import SimpleDomain
+from sympy.polys.polyerrors import CoercionFailed
+from sympy.utilities import public
+
+import math
+
+@public
+class IntegerRing(Ring, CharacteristicZero, SimpleDomain):
+    r"""The domain ``ZZ`` representing the integers `\mathbb{Z}`.
+
+    The :py:class:`IntegerRing` class represents the ring of integers as a
+    :py:class:`~.Domain` in the domain system. :py:class:`IntegerRing` is a
+    super class of :py:class:`PythonIntegerRing` and
+    :py:class:`GMPYIntegerRing` one of which will be the implementation for
+    :ref:`ZZ` depending on whether or not ``gmpy`` or ``gmpy2`` is installed.
+
+    See also
+    ========
+
+    Domain
+    """
+
+    rep = 'ZZ'
+    alias = 'ZZ'
+    dtype = MPZ
+    zero = dtype(0)
+    one = dtype(1)
+    tp = type(one)
+
+
+    is_IntegerRing = is_ZZ = True
+    is_Numerical = True
+    is_PID = True
+
+    has_assoc_Ring = True
+    has_assoc_Field = True
+
+    def __init__(self):
+        """Allow instantiation of this domain. """
+
+    def __eq__(self, other):
+        """Returns ``True`` if two domains are equivalent. """
+        if isinstance(other, IntegerRing):
+            return True
+        else:
+            return NotImplemented
+
+    def __hash__(self):
+        """Compute a hash value for this domain. """
+        return hash('ZZ')
+
+    def to_sympy(self, a):
+        """Convert ``a`` to a SymPy object. """
+        return SymPyInteger(int(a))
+
+    def from_sympy(self, a):
+        """Convert SymPy's Integer to ``dtype``. """
+        if a.is_Integer:
+            return MPZ(a.p)
+        elif int_valued(a):
+            return MPZ(int(a))
+        else:
+            raise CoercionFailed("expected an integer, got %s" % a)
+
+    def get_field(self):
+        r"""Return the associated field of fractions :ref:`QQ`
+
+        Returns
+        =======
+
+        :ref:`QQ`:
+            The associated field of fractions :ref:`QQ`, a
+            :py:class:`~.Domain` representing the rational numbers
+            `\mathbb{Q}`.
+
+        Examples
+        ========
+
+        >>> from sympy import ZZ
+        >>> ZZ.get_field()
+        QQ
+        """
+        from sympy.polys.domains import QQ
+        return QQ
+
+    def algebraic_field(self, *extension, alias=None):
+        r"""Returns an algebraic field, i.e. `\mathbb{Q}(\alpha, \ldots)`.
+
+        Parameters
+        ==========
+
+        *extension : One or more :py:class:`~.Expr`.
+            Generators of the extension. These should be expressions that are
+            algebraic over `\mathbb{Q}`.
+
+        alias : str, :py:class:`~.Symbol`, None, optional (default=None)
+            If provided, this will be used as the alias symbol for the
+            primitive element of the returned :py:class:`~.AlgebraicField`.
+
+        Returns
+        =======
+
+        :py:class:`~.AlgebraicField`
+            A :py:class:`~.Domain` representing the algebraic field extension.
+
+        Examples
+        ========
+
+        >>> from sympy import ZZ, sqrt
+        >>> ZZ.algebraic_field(sqrt(2))
+        QQ<sqrt(2)>
+        """
+        return self.get_field().algebraic_field(*extension, alias=alias)
+
+    def from_AlgebraicField(K1, a, K0):
+        """Convert a :py:class:`~.ANP` object to :ref:`ZZ`.
+
+        See :py:meth:`~.Domain.convert`.
+        """
+        if a.is_ground:
+            return K1.convert(a.LC(), K0.dom)
+
+    def log(self, a, b):
+        r"""Logarithm of *a* to the base *b*.
+
+        Parameters
+        ==========
+
+        a: number
+        b: number
+
+        Returns
+        =======
+
+        $\\lfloor\log(a, b)\\rfloor$:
+            Floor of the logarithm of *a* to the base *b*
+
+        Examples
+        ========
+
+        >>> from sympy import ZZ
+        >>> ZZ.log(ZZ(8), ZZ(2))
+        3
+        >>> ZZ.log(ZZ(9), ZZ(2))
+        3
+
+        Notes
+        =====
+
+        This function uses ``math.log`` which is based on ``float`` so it will
+        fail for large integer arguments.
+        """
+        return self.dtype(int(math.log(int(a), b)))
+
+    def from_FF(K1, a, K0):
+        """Convert ``ModularInteger(int)`` to GMPY's ``mpz``. """
+        return MPZ(K0.to_int(a))
+
+    def from_FF_python(K1, a, K0):
+        """Convert ``ModularInteger(int)`` to GMPY's ``mpz``. """
+        return MPZ(K0.to_int(a))
+
+    def from_ZZ(K1, a, K0):
+        """Convert Python's ``int`` to GMPY's ``mpz``. """
+        return MPZ(a)
+
+    def from_ZZ_python(K1, a, K0):
+        """Convert Python's ``int`` to GMPY's ``mpz``. """
+        return MPZ(a)
+
+    def from_QQ(K1, a, K0):
+        """Convert Python's ``Fraction`` to GMPY's ``mpz``. """
+        if a.denominator == 1:
+            return MPZ(a.numerator)
+
+    def from_QQ_python(K1, a, K0):
+        """Convert Python's ``Fraction`` to GMPY's ``mpz``. """
+        if a.denominator == 1:
+            return MPZ(a.numerator)
+
+    def from_FF_gmpy(K1, a, K0):
+        """Convert ``ModularInteger(mpz)`` to GMPY's ``mpz``. """
+        return MPZ(K0.to_int(a))
+
+    def from_ZZ_gmpy(K1, a, K0):
+        """Convert GMPY's ``mpz`` to GMPY's ``mpz``. """
+        return a
+
+    def from_QQ_gmpy(K1, a, K0):
+        """Convert GMPY ``mpq`` to GMPY's ``mpz``. """
+        if a.denominator == 1:
+            return a.numerator
+
+    def from_RealField(K1, a, K0):
+        """Convert mpmath's ``mpf`` to GMPY's ``mpz``. """
+        p, q = K0.to_rational(a)
+
+        if q == 1:
+            # XXX: If MPZ is flint.fmpz and p is a gmpy2.mpz, then we need
+            # to convert via int because fmpz and mpz do not know about each
+            # other.
+            return MPZ(int(p))
+
+    def from_GaussianIntegerRing(K1, a, K0):
+        if a.y == 0:
+            return a.x
+
+    def from_EX(K1, a, K0):
+        """Convert ``Expression`` to GMPY's ``mpz``. """
+        if a.is_Integer:
+            return K1.from_sympy(a)
+
+    def gcdex(self, a, b):
+        """Compute extended GCD of ``a`` and ``b``. """
+        h, s, t = gcdex(a, b)
+        # XXX: This conditional logic should be handled somewhere else.
+        if GROUND_TYPES == 'gmpy':
+            return s, t, h
+        else:
+            return h, s, t
+
+    def gcd(self, a, b):
+        """Compute GCD of ``a`` and ``b``. """
+        return gcd(a, b)
+
+    def lcm(self, a, b):
+        """Compute LCM of ``a`` and ``b``. """
+        return lcm(a, b)
+
+    def sqrt(self, a):
+        """Compute square root of ``a``. """
+        return sqrt(a)
+
+    def is_square(self, a):
+        """Return ``True`` if ``a`` is a square.
+
+        Explanation
+        ===========
+        An integer is a square if and only if there exists an integer
+        ``b`` such that ``b * b == a``.
+        """
+        return is_square(a)
+
+    def exsqrt(self, a):
+        """Non-negative square root of ``a`` if ``a`` is a square.
+
+        See also
+        ========
+        is_square
+        """
+        if a < 0:
+            return None
+        root, rem = sqrtrem(a)
+        if rem != 0:
+            return None
+        return root
+
+    def factorial(self, a):
+        """Compute factorial of ``a``. """
+        return factorial(a)
+
+
+ZZ = IntegerRing()