Repository: sympy/sympy
Pull request: feat(polys): Add Gaussian domains to Poly
<!-- Your title above should be a short description of what
was changed. Do not include the issue number in the title. -->

#### References to other Issues or PRs
<!-- If this pull request fixes an issue, write "Fixes #NNNN" in that exact
format, e.g. "Fixes #1234" (see
https://tinyurl.com/auto-closing for more information). Also, please
write a comment on that issue linking back to this pull request once it is
open. -->

Fixes #19579
Follows on from #15396


#### Brief description of what is fixed or changed

Add the Gaussian domains to Poly so that they are used when appropriate e.g.:
```julia
In [1]: Poly(I*x, x)                                                                                                                           
Out[1]: Poly(I*x, x, domain='ZZ_I')
```

#### Other comments

There are still some problems here most notably this line
https://github.com/sympy/sympy/blob/46469d714a095fde59f9c4d83aaa528c68b0be0f/sympy/polys/euclidtools.py#L996
because `is_negative` does not work in the Gaussian domains.

Also there needs to be an `lcm` function for Gaussian domains.

#### Release Notes

<!-- Write the release notes for this release below. See
https://github.com/sympy/sympy/wiki/Writing-Release-Notes for more information
on how to write release notes. The bot will check your release notes
automatically to see if they are formatted correctly. -->

<!-- BEGIN RELEASE NOTES -->
NO ENTRY
<!-- END RELEASE NOTES -->

Problem statement:
Gaussian integer division raises TypeError
Gaussian integer division raises error when `auto=True` is enabled, which is default behavior. 
```python3
>>> from sympy import *
>>> from sympy.polys.domains.gaussiandomains import ZZ_I, QQ_I
>>> x = Symbol('x')
>>> p = Poly(2+3*I, x, domain=ZZ_I)
>>> q = Poly(1-I, x, domain=ZZ_I)
>>> p.div(q, auto=False)
(Poly(0, x, domain='ZZ_I'), Poly(2 + 3*I, x, domain='ZZ_I'))
>>> p.div(q, auto=True)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "C:\Users\sylee\Documents\Github\sympy\sympy\polys\polytools.py", line 1663, in div
    F, G = F.to_field(), G.to_field()
  File "C:\Users\sylee\Documents\Github\sympy\sympy\polys\polyclasses.py", line 289, in to_field
    return f.convert(f.dom.get_field())
  File "C:\Users\sylee\Documents\Github\sympy\sympy\polys\polyclasses.py", line 300, in convert
    return DMP(dmp_convert(f.rep, f.lev, f.dom, dom), dom, f.lev)
  File "C:\Users\sylee\Documents\Github\sympy\sympy\polys\densebasic.py", line 562, in dmp_convert
    return dup_convert(f, K0, K1)
  File "C:\Users\sylee\Documents\Github\sympy\sympy\polys\densebasic.py", line 539, in dup_convert
    return dup_strip([ K1.convert(c, K0) for c in f ])
  File "C:\Users\sylee\Documents\Github\sympy\sympy\polys\densebasic.py", line 539, in <listcomp>
    return dup_strip([ K1.convert(c, K0) for c in f ])
TypeError: convert() takes 2 positional arguments but 3 were given
```

It works fine with Gaussian rationals. 
```python3
>>> p = Poly(2+3*I, x, domain=QQ_I)
>>> q = Poly(1-I, x, domain=QQ_I)
>>> p.div(q, auto=True)
(Poly(-1/2 + 5*I/2, x, domain='QQ_I'), Poly(0, x, domain='QQ_I'))
>>> p.div(q, auto=False)
(Poly(-1/2 + 5*I/2, x, domain='QQ_I'), Poly(0, x, domain='QQ_I'))
```

Hints from issue discussion:

Commit message: feat(polys): Add Gaussian domains to Poly

Use AlgebraicField instead of the Gaussian domains for factoring in
apart. The factorisation algorithms are not provided for Gaussian
domains yet.

Patch to review:
diff --git a/sympy/__init__.py b/sympy/__init__.py
index 4603b62f94eb..d7d1863af367 100644
--- a/sympy/__init__.py
+++ b/sympy/__init__.py
@@ -95,7 +95,7 @@ def __sympy_debug():
         GMPYFiniteField, PythonIntegerRing, GMPYIntegerRing, PythonRational,
         GMPYRationalField, AlgebraicField, PolynomialRing, FractionField,
         ExpressionDomain, FF_python, FF_gmpy, ZZ_python, ZZ_gmpy, QQ_python,
-        QQ_gmpy, GF, FF, ZZ, QQ, RR, CC, EX, construct_domain,
+        QQ_gmpy, GF, FF, ZZ, QQ, ZZ_I, QQ_I, RR, CC, EX, construct_domain,
         swinnerton_dyer_poly, cyclotomic_poly, symmetric_poly, random_poly,
         interpolating_poly, jacobi_poly, chebyshevt_poly, chebyshevu_poly,
         hermite_poly, legendre_poly, laguerre_poly, apart, apart_list,
@@ -314,7 +314,7 @@ def __sympy_debug():
     'GMPYIntegerRing', 'PythonRational', 'GMPYRationalField',
     'AlgebraicField', 'PolynomialRing', 'FractionField', 'ExpressionDomain',
     'FF_python', 'FF_gmpy', 'ZZ_python', 'ZZ_gmpy', 'QQ_python', 'QQ_gmpy',
-    'GF', 'FF', 'ZZ', 'QQ', 'RR', 'CC', 'EX', 'construct_domain',
+    'GF', 'FF', 'ZZ', 'QQ', 'ZZ_I', 'QQ_I', 'RR', 'CC', 'EX', 'construct_domain',
     'swinnerton_dyer_poly', 'cyclotomic_poly', 'symmetric_poly',
     'random_poly', 'interpolating_poly', 'jacobi_poly', 'chebyshevt_poly',
     'chebyshevu_poly', 'hermite_poly', 'legendre_poly', 'laguerre_poly',
diff --git a/sympy/polys/__init__.py b/sympy/polys/__init__.py
index 40ec0acb5879..0ad0d9014290 100644
--- a/sympy/polys/__init__.py
+++ b/sympy/polys/__init__.py
@@ -41,7 +41,8 @@
     'PythonIntegerRing', 'GMPYIntegerRing', 'PythonRational',
     'GMPYRationalField', 'AlgebraicField', 'PolynomialRing', 'FractionField',
     'ExpressionDomain', 'FF_python', 'FF_gmpy', 'ZZ_python', 'ZZ_gmpy',
-    'QQ_python', 'QQ_gmpy', 'GF', 'FF', 'ZZ', 'QQ', 'RR', 'CC', 'EX',
+    'QQ_python', 'QQ_gmpy', 'GF', 'FF', 'ZZ', 'QQ', 'ZZ_I', 'QQ_I', 'RR',
+    'CC', 'EX',
 
     'construct_domain',
 
@@ -101,7 +102,7 @@
         PythonIntegerRing, GMPYIntegerRing, PythonRational, GMPYRationalField,
         AlgebraicField, PolynomialRing, FractionField, ExpressionDomain,
         FF_python, FF_gmpy, ZZ_python, ZZ_gmpy, QQ_python, QQ_gmpy, GF, FF,
-        ZZ, QQ, RR, CC, EX)
+        ZZ, QQ, ZZ_I, QQ_I, RR, CC, EX)
 
 from .constructor import construct_domain
 
diff --git a/sympy/polys/constructor.py b/sympy/polys/constructor.py
index 12d9de96e9ae..a04306741408 100644
--- a/sympy/polys/constructor.py
+++ b/sympy/polys/constructor.py
@@ -3,7 +3,9 @@
 from __future__ import print_function, division
 
 from sympy.core import sympify
+from sympy.core.evalf import pure_complex
 from sympy.polys.domains import ZZ, QQ, EX
+from sympy.polys.domains.gaussiandomains import ZZ_I, QQ_I
 from sympy.polys.domains.realfield import RealField
 from sympy.polys.polyoptions import build_options
 from sympy.polys.polyutils import parallel_dict_from_basic
@@ -12,33 +14,44 @@
 
 def _construct_simple(coeffs, opt):
     """Handle simple domains, e.g.: ZZ, QQ, RR and algebraic domains. """
-    result, rationals, reals, algebraics = {}, False, False, False
+    rationals = reals = gaussians = algebraics = False
 
     if opt.extension is True:
         is_algebraic = lambda coeff: coeff.is_number and coeff.is_algebraic
     else:
         is_algebraic = lambda coeff: False
 
-    # XXX: add support for a + b*I coefficients
     for coeff in coeffs:
         if coeff.is_Rational:
             if not coeff.is_Integer:
                 rationals = True
         elif coeff.is_Float:
-            if not algebraics:
-                reals = True
-            else:
+            if algebraics or gaussians:
                 # there are both reals and algebraics -> EX
                 return False
-        elif is_algebraic(coeff):
-            if not reals:
-                algebraics = True
             else:
-                # there are both algebraics and reals -> EX
-                return False
+                reals = True
         else:
-            # this is a composite domain, e.g. ZZ[X], EX
-            return None
+            is_complex = pure_complex(coeff)
+            if is_complex:
+                # there are both reals and algebraics -> EX
+                if reals:
+                    return False
+                x, y = is_complex
+                if x.is_Rational and y.is_Rational:
+                    gaussians = True
+                    if not (x.is_Integer and y.is_Integer):
+                        rationals = True
+                    continue
+            if is_algebraic(coeff):
+                if not reals:
+                    algebraics = True
+                else:
+                    # there are both algebraics and reals -> EX
+                    return False
+            else:
+                # this is a composite domain, e.g. ZZ[X], EX
+                return None
 
     if algebraics:
         domain, result = _construct_algebraic(coeffs, opt)
@@ -50,9 +63,9 @@ def _construct_simple(coeffs, opt):
             domain = RealField(prec=max_prec)
         else:
             if opt.field or rationals:
-                domain = QQ
+                domain = QQ_I if gaussians else QQ
             else:
-                domain = ZZ
+                domain = ZZ_I if gaussians else ZZ
 
         result = []
 
diff --git a/sympy/polys/domains/__init__.py b/sympy/polys/domains/__init__.py
index 928938af6c9e..3263e68545ff 100644
--- a/sympy/polys/domains/__init__.py
+++ b/sympy/polys/domains/__init__.py
@@ -1,11 +1,15 @@
 """Implementation of mathematical domains. """
 
-__all__ = ['Domain', 'FiniteField', 'IntegerRing', 'RationalField',
-           'RealField', 'ComplexField', 'PythonFiniteField',
-           'GMPYFiniteField', 'PythonIntegerRing',
-           'GMPYIntegerRing', 'PythonRational',
-           'GMPYRationalField', 'AlgebraicField', 'PolynomialRing',
-           'FractionField', 'ExpressionDomain', 'PythonRational']
+__all__ = [
+    'Domain', 'FiniteField', 'IntegerRing', 'RationalField', 'RealField',
+    'ComplexField', 'PythonFiniteField', 'GMPYFiniteField',
+    'PythonIntegerRing', 'GMPYIntegerRing', 'PythonRational',
+    'GMPYRationalField', 'AlgebraicField', 'PolynomialRing', 'FractionField',
+    'ExpressionDomain', 'PythonRational',
+
+    'FF_python', 'FF_gmpy', 'ZZ_python', 'ZZ_gmpy', 'QQ_python', 'QQ_gmpy',
+    'GF', 'FF', 'ZZ', 'QQ', 'ZZ_I', 'QQ_I', 'RR', 'CC', 'EX',
+]
 
 from .domain import Domain
 from .finitefield import FiniteField
@@ -49,13 +53,9 @@
 except KeyError:
     raise ValueError("invalid ground types: %s" % GROUND_TYPES)
 
+# Needs to be after ZZ and QQ are defined:
+from .gaussiandomains import ZZ_I, QQ_I
+
 GF = FF
 
 EX = ExpressionDomain()
-
-__all__.extend([
-    "FF_python", "FF_gmpy",
-    "ZZ_python", "ZZ_gmpy",
-    "QQ_python", "QQ_gmpy",
-    "GF", "FF", "ZZ", "QQ", "RR", "CC", "EX",
-])
diff --git a/sympy/polys/domains/algebraicfield.py b/sympy/polys/domains/algebraicfield.py
index 028a29e7a811..ec1135dc1425 100644
--- a/sympy/polys/domains/algebraicfield.py
+++ b/sympy/polys/domains/algebraicfield.py
@@ -130,3 +130,7 @@ def denom(self, a):
     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 GaussianInteger element 'a' to an AlgebraicField """
+        return K1.from_sympy(K0.to_sympy(a))
diff --git a/sympy/polys/domains/domain.py b/sympy/polys/domains/domain.py
index 9af62b6a8f17..31b4f64c8145 100644
--- a/sympy/polys/domains/domain.py
+++ b/sympy/polys/domains/domain.py
@@ -30,6 +30,8 @@ class Domain(object):
     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
@@ -324,6 +326,20 @@ def mkinexact(cls, K0, K1):
         elif K1.is_AlgebraicField:
             return K1
 
+        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:
diff --git a/sympy/polys/domains/expressiondomain.py b/sympy/polys/domains/expressiondomain.py
index 9d3c88bddb11..a610775fac25 100644
--- a/sympy/polys/domains/expressiondomain.py
+++ b/sympy/polys/domains/expressiondomain.py
@@ -173,6 +173,14 @@ 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_RealField(K1, a, K0):
         """Convert a mpmath ``mpf`` object to ``dtype``. """
         return K1(K0.to_sympy(a))
diff --git a/sympy/polys/domains/fractionfield.py b/sympy/polys/domains/fractionfield.py
index 1221f5359078..53b9e9b083f1 100644
--- a/sympy/polys/domains/fractionfield.py
+++ b/sympy/polys/domains/fractionfield.py
@@ -93,6 +93,14 @@ 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))
diff --git a/sympy/polys/domains/gaussiandomains.py b/sympy/polys/domains/gaussiandomains.py
index ee1bc87f094a..80eb0728f599 100644
--- a/sympy/polys/domains/gaussiandomains.py
+++ b/sympy/polys/domains/gaussiandomains.py
@@ -1,9 +1,10 @@
 """Domains of Gaussian type."""
 
 from sympy.core.basic import Basic
-from sympy.core.numbers import I
+from sympy.core.numbers import Rational, I
 from sympy.polys.polyerrors import CoercionFailed
 from sympy.polys.domains import ZZ, QQ
+from sympy.polys.domains.algebraicfield import AlgebraicField
 from sympy.polys.domains.domainelement import DomainElement
 from sympy.polys.domains.field import Field
 from sympy.polys.domains.ring import Ring
@@ -16,11 +17,14 @@ class GaussianElement(DomainElement):
 
     __slots__ = ('x', 'y')
 
-    def __init__(self, x, y):
+    def __init__(self, x, y=0):
         conv = self.base.convert
         self.x = conv(x)
         self.y = conv(y)
 
+    def as_expr(self):
+        return Rational(self.x) + I * Rational(self.y)
+
     def parent(self):
         return self._parent
 
@@ -79,9 +83,21 @@ def __mul__(self, other):
         if x is not None:
             return self.__class__(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.__class__(1, 0)
+        if exp < 0:
+            self, exp = 1/self, -exp
+        prod = self
+        for n in range(exp-1):
+            prod *= self
+        return prod
+
     def __bool__(self):
         return bool(self.x) or bool(self.y)
 
@@ -239,7 +255,7 @@ def from_sympy(self, a):
         else:
             raise CoercionFailed("{} is not Gaussian".format(a))
 
-    def convert(self, element):
+    def convert(self, element, base=None):
         """Convert ``element`` to ``self.dtype``.
 
         Raises CoercionFailed on failure.
@@ -253,6 +269,32 @@ def convert(self, element):
         else:  # convertible to base type or failure
             return self.new(element, 0)
 
+    def inject(self, *gens):
+        """Inject generators into this domain. """
+        return self.poly_ring(*gens)
+
+    def as_AlgebraicField(self):
+        """Get equivalent domain as an ``AlgebraicField``. """
+        return AlgebraicField(self.base, I)
+
+    # Override the negative etc handlers because this isn't an ordered domain.
+
+    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
+
 
 class GaussianIntegerRing(GaussianDomain, Ring):
     """Ring of Gaussian integers."""
@@ -265,6 +307,9 @@ class GaussianIntegerRing(GaussianDomain, Ring):
 
     rep = 'ZZ_I'
 
+    is_GaussianRing = True
+    is_ZZ_I = True
+
     def __init__(self):  # override Domain.__init__
         """For constructing ZZ_I."""
 
@@ -291,6 +336,9 @@ def gcd(self, a, b):
             a, b = b, a % b
         return self.normalize(a)
 
+    def lcm(self, a, b):
+        return (a * b) // self.gcd(a, b)
+
 ZZ_I = GaussianInteger._parent = GaussianIntegerRing()
 
 
@@ -303,6 +351,9 @@ class GaussianRationalField(GaussianDomain, Field):
 
     rep = 'QQ_I'
 
+    is_GaussianField = True
+    is_QQ_I = True
+
     def __init__(self):  # override Domain.__init__
         """For constructing QQ_I."""
 
@@ -314,4 +365,12 @@ def get_field(self):
         """Returns a field associated with ``self``. """
         return self
 
+    def denom(self, a):
+        """Get the denominator of ``a``."""
+        ZZ = self.base.get_ring()
+        QQ = self.base
+        ZZ_I = self.get_ring()
+        denom_ZZ = ZZ.lcm(QQ.denom(a.x), QQ.denom(a.y))
+        return ZZ_I(denom_ZZ, ZZ.zero)
+
 QQ_I = GaussianRational._parent = GaussianRationalField()
diff --git a/sympy/polys/domains/gmpyintegerring.py b/sympy/polys/domains/gmpyintegerring.py
index 523d9a6d8fdc..20d24735557f 100644
--- a/sympy/polys/domains/gmpyintegerring.py
+++ b/sympy/polys/domains/gmpyintegerring.py
@@ -70,6 +70,10 @@ def from_RealField(K1, a, K0):
         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)
diff --git a/sympy/polys/domains/gmpyrationalfield.py b/sympy/polys/domains/gmpyrationalfield.py
index 086e659ff354..069e42aae53e 100644
--- a/sympy/polys/domains/gmpyrationalfield.py
+++ b/sympy/polys/domains/gmpyrationalfield.py
@@ -59,6 +59,11 @@ 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)))
diff --git a/sympy/polys/domains/polynomialring.py b/sympy/polys/domains/polynomialring.py
index d23bafea0ee2..cb499dce626b 100644
--- a/sympy/polys/domains/polynomialring.py
+++ b/sympy/polys/domains/polynomialring.py
@@ -92,6 +92,14 @@ def from_QQ_gmpy(K1, a, K0):
         """Convert a GMPY `mpq` 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_GaussianRationalField(K1, a, K0):
+        """Convert a `GaussianRational` 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))
diff --git a/sympy/polys/partfrac.py b/sympy/polys/partfrac.py
index ee5ffdc6403c..183759298f2e 100644
--- a/sympy/polys/partfrac.py
+++ b/sympy/polys/partfrac.py
@@ -116,6 +116,16 @@ def apart(f, x=None, full=False, **options):
                     pass
             return f.xreplace(dict(reps))
 
+    # Factorization over ZZ_I or QQ_I is not supported so we use the
+    # factorization routines for AlgebraicField instead.
+    dom = P.domain
+    if dom.is_ZZ_I or dom.is_QQ_I:
+        if dom.is_ZZ_I:
+            dom = dom.get_field()
+        dom_alg = dom.as_AlgebraicField()
+        P = P.per(P.rep.convert(dom_alg))
+        Q = Q.per(Q.rep.convert(dom_alg))
+
     if P.is_multivariate:
         fc = f.cancel()
         if fc != f:
diff --git a/sympy/polys/polyoptions.py b/sympy/polys/polyoptions.py
index 399588be0bce..d138ea06b84b 100644
--- a/sympy/polys/polyoptions.py
+++ b/sympy/polys/polyoptions.py
@@ -423,6 +423,12 @@ def preprocess(cls, domain):
             if domain in ['Q', 'QQ']:
                 return sympy.polys.domains.QQ
 
+            if domain == 'ZZ_I':
+                return sympy.polys.domains.ZZ_I
+
+            if domain == 'QQ_I':
+                return sympy.polys.domains.QQ_I
+
             if domain == 'EX':
                 return sympy.polys.domains.EX
 
diff --git a/sympy/polys/polyroots.py b/sympy/polys/polyroots.py
index 979b27483bee..9f985b176fa5 100644
--- a/sympy/polys/polyroots.py
+++ b/sympy/polys/polyroots.py
@@ -18,6 +18,7 @@
 from sympy.functions import exp, sqrt, im, cos, acos, Piecewise
 from sympy.functions.elementary.miscellaneous import root
 from sympy.ntheory import divisors, isprime, nextprime
+from sympy.polys.domains import EX
 from sympy.polys.polyerrors import (PolynomialError, GeneratorsNeeded,
     DomainError)
 from sympy.polys.polyquinticconst import PolyQuintic
@@ -970,6 +971,10 @@ def _try_heuristics(f):
 
         return result
 
+    # Convert the generators to symbols
+    dumgens = symbols('x:%d' % len(f.gens), cls=Dummy)
+    f = f.per(f.rep, dumgens)
+
     (k,), f = f.terms_gcd()
 
     if not k:
@@ -982,6 +987,10 @@ def _try_heuristics(f):
     if auto and f.get_domain().is_Ring:
         f = f.to_field()
 
+    # Use EX instead of ZZ_I or QQ_I
+    if f.get_domain().is_QQ_I:
+        f = f.per(f.rep.convert(EX))
+
     rescale_x = None
     translate_x = None
 
diff --git a/sympy/printing/tests/test_str.py b/sympy/printing/tests/test_str.py
index 911a0b10ad6c..42884149182c 100644
--- a/sympy/printing/tests/test_str.py
+++ b/sympy/printing/tests/test_str.py
@@ -372,8 +372,8 @@ def test_Poly():
     assert str(
         Poly(x**2 - 1 + y, x)) == "Poly(x**2 + y - 1, x, domain='ZZ[y]')"
 
-    assert str(Poly(x**2 + I*x, x)) == "Poly(x**2 + I*x, x, domain='EX')"
-    assert str(Poly(x**2 - I*x, x)) == "Poly(x**2 - I*x, x, domain='EX')"
+    assert str(Poly(x**2 + I*x, x)) == "Poly(x**2 + I*x, x, domain='ZZ_I')"
+    assert str(Poly(x**2 - I*x, x)) == "Poly(x**2 - I*x, x, domain='ZZ_I')"
 
     assert str(Poly(-x*y*z + x*y - 1, x, y, z)
                ) == "Poly(-x*y*z + x*y - 1, x, y, z, domain='ZZ')"


Code review comments:
diff_hunk: @@ -116,6 +116,16 @@ def apart(f, x=None, full=False, **options):
                     pass
             return f.xreplace(dict(reps))
 
+    # Factorization over ZZ_I or QQ_I is not supported so we use the
+    # factorization routines for AlgebraicField instead.
+    dom = P.domain
+    if dom.is_ZZ_I or dom.is_QQ_I:
+        if dom.is_ZZ_I:
+            dom = dom.get_field()
+        dom_alg = dom.as_AlgebraicField()
+        P = P.per(P.rep.convert(dom_alg))
+        Q = Q.per(Q.rep.convert(dom_alg))
path: sympy/polys/partfrac.py
text: @author:
Factorisation over `ZZ_I` and `QQ_I` does not work so I've prevented them from being used in `apart` here. Instead `QQ<I>` will be used.
```julia
In [16]: p1 = Poly(x**2 + 1, x)                                                                                                                

In [17]: p2 = Poly(x**2 + 1, x, domain=QQ_I.as_AlgebraicField())                                                                               

In [18]: p1.factor_list()                                                                                                                      
Out[18]: (1, [(Poly(x**2 + 1, x, domain='ZZ'), 1)])

In [19]: p2.factor_list()                                                                                                                      
Out[19]: (1, [(Poly(x - I, x, domain='QQ<I>'), 1), (Poly(x + I, x, domain='QQ<I>'), 1)])
```
Would it be better to make factorisation over `QQ_I` and `ZZ_I` work?

I wonder if it's possible to make `dup/dmp` treat `QQ_I` as an algebraic field when factorising. I guess that factorising over `ZZ_I` is harder...


@user1:
> Would it be better to make factorisation over `QQ_I` and `ZZ_I` work?

I think that `QQ_I` should be replaced by `QQ<I>` somewhere, maybe in `dup/dmp_ext_factor`.

@author:
Would it be possible to add attributes or methods to `QQ_I` so that `dup/dmp` can treat it transparently as an algebraic field without needing any conversion?

Would it make sense to have `QQ_I.is_AlgebraicField` or some other attribute that is common between `QQ_I` and the `AlgebraicField` domains?

It looks like the main attributes needed for an algebraic field are `ext` and `mod`...

@user1:
It looks like at least new versions of `dup/dmp_sqf_norm` would be needed. `QQ_I` should have `conjugate` and the norm of `f` will be `f*cf` where `cf` is `f` with all coefficients conjugated.

@author:
Would it make sense for `ZZ_I` to factorise like `QQ_I` or perhaps like `ZZ[I]`?

On master factor can factorise in some cases by treating `I` as a generator:
```julia
In [14]: e = expand((x - I)*(x + 1))                                                                                                           

In [15]: e                                                                                                                                     
Out[15]: 
 2              
x  + x - ⅈ⋅x - ⅈ

In [16]: factor(e)                                                                                                                             
Out[16]: (x + 1)⋅(x - ⅈ)
```
The PR currently gives:
```julia
In [14]: factor(e)                                                                                                                             
Out[14]: 
 2                
x  + x⋅(1 - ⅈ) - ⅈ
```

@user1:
Factorizing polynomials over `ZZ_I` should be done in two steps, factorizing the primitive part over `QQ_I` and the content in `ZZ_I`.

@author:
What would `factor_list` return in that case?

At the moment `factor_list` returns a tuple of the content and then a list of factors of the primitive:
```julia
In [21]: Poly(x**2/2 - 2, x).factor_list()                                                                                                     
Out[21]: (1/2, [(Poly(x - 2, x, domain='QQ'), 1), (Poly(x + 2, x, domain='QQ'), 1)])
```
If the content itself factorises where would that go in this data structure?

Is there an implementation for factorisation of gaussian integers anywhere?

Perhaps if not then for now `factor_list` can return the content unfactorised...

@user1:
> Is there an implementation for factorisation of gaussian integers anywhere?

No, but that would not be hard to implement if desired. First factorize the norm in `ZZ` and then factor those prime factors that split in `ZZ_I`.