domains: Added Gaussian integers and rationals

[jksuom]

References to other Issues or PRs

#15358 (comment)

Brief description of what is fixed or changed

An implementation Gaussian domains including gcd is added.

Other comments

Release Notes

• polys
  • added new domains for Gaussian integers and rationals

[sympy-bot]

✅

Hi, I am the SymPy bot (v160). I'm here to help you write a release notes entry. Please read the guide on how to write release notes.

Your release notes are in good order.

Here is what the release notes will look like:

• polys
  • added new domains for Gaussian integers and rationals (#15396 by @jksuom and @smichr)

This will be added to https://github.com/sympy/sympy/wiki/Release-Notes-for-1.7.

Note: This comment will be updated with the latest check if you edit the pull request. You need to reload the page to see it.

Click here to see the pull request description that was parsed.

<!-- 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://github.com/blog/1506-closing-issues-via-pull-requests .-->
https://github.com/sympy/sympy/issues/15358#issuecomment-428282577

#### Brief description of what is fixed or changed
An implementation Gaussian domains including gcd is added.

#### Other comments

#### 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 -->
* polys
    * added new domains for Gaussian integers and rationals
<!-- END RELEASE NOTES -->

Update

The release notes on the wiki have been updated.

[asmeurer]

Everything looks OK. It would be good to have lots more tests, though, especially for the code that can produce wrong results from typos.

[czgdp1807]

@jksuom Is this PR still having a chance to be merged? One, option can be is to raise an issue to keep track of this PR.

[jksuom]

More tests are needed. After those pass, I think that this could be merged.

[czgdp1807]

@jksuom Please add more tests to this PR, so that we can merge this. Or, someone familiar with the module can extend the PR.

[czgdp1807]

ping @jksuom Any updates.

[czgdp1807]

Conflicts need to be resolved for this PR.

[sympy-bot]

🟠

Hi, I am the SymPy bot (v158). I've noticed that some of your commits add or delete files. Since this is sometimes done unintentionally, I wanted to alert you about it.

This is an experimental feature of SymPy Bot. If you have any feedback on it, please comment at sympy/sympy-bot#75.

The following commits add new files:

• 2abbd19:
  • sympy/polys/domains/gaussiandomains.py

If these files were added/deleted on purpose, you can ignore this message.

[codecov]

Codecov Report

Merging #15396 into master will increase coverage by 0.026%.
The diff coverage is 90.500%.

@@              Coverage Diff              @@
##            master    #15396       +/-   ##
=============================================
+ Coverage   75.660%   75.686%   +0.026%     
=============================================
  Files          647       648        +1     
  Lines       168517    168720      +203     
  Branches     39707     39729       +22     
=============================================
+ Hits        127500    127698      +198     
- Misses       35455     35456        +1     
- Partials      5562      5566        +4     

[smichr]

Only the NotImplemented are not covered.

[smichr]

@jksuom , could you please check the XXX in the file -- I had a question there.

[jksuom]

This should be the field, i.e. QQ_I.

[smichr]

The field of ZZ_I is QQ_I and the field of QQ_I is ZZ_I (as you had it before)?

[jksuom]

ZZ_I is the ring of QQ_I. It should be returned by QQ_I.get_ring().

[smichr]

In summary, here are the minor edits that I made to the orginal PR

diff --git a/sympy/polys/domains/gaussiandomains.py b/sympy/polys/domains/gaussiandomains.py
index 95d5df2..ee1bc87 100644
--- a/sympy/polys/domains/gaussiandomains.py
+++ b/sympy/polys/domains/gaussiandomains.py
@@ -1,7 +1,5 @@
 """Domains of Gaussian type."""
 
-from __future__ import print_function, division
-
 from sympy.core.basic import Basic
 from sympy.core.numbers import I
 from sympy.polys.polyerrors import CoercionFailed
@@ -118,7 +116,7 @@ def __rtruediv__(self, other):
 
     def __divmod__(self, other):
         if not other:
-            raise ZeroDivisionError('divmod(%s, 0)' % self)
+            raise ZeroDivisionError('divmod({}, 0)'.format(self))
         x, y = self._get_xy(other)
         if x is None:
             return NotImplemented
@@ -153,7 +151,7 @@ def __floordiv__(self, other):
         return qr if qr is NotImplemented else qr[0]
 
     def __rfloordiv__(self, other):
-        qr = self._rdivmod__(self, other)
+        qr = self.__rdivmod__(other)
         return qr if qr is NotImplemented else qr[0]
 
     def __mod__(self, other):
@@ -161,7 +159,7 @@ def __mod__(self, other):
         return qr if qr is NotImplemented else qr[1]
 
     def __rmod__(self, other):
-        qr = self._rdivmod__(self, other)
+        qr = self.__rdivmod__(other)
         return qr if qr is NotImplemented else qr[1]
 
 
@@ -170,7 +168,7 @@ class GaussianRational(GaussianElement):
 
     def __truediv__(self, other):
         if not other:
-            raise ZeroDivisionError('%s / 0' % self)
+            raise ZeroDivisionError('{} / 0'.format(self))
         x, y = self._get_xy(other)
         if x is None:
             return NotImplemented
@@ -184,7 +182,7 @@ def __truediv__(self, other):
     def __rtruediv__(self, other):
         try:
             other = self._parent.convert(other)
-        except CoercienFailed:
+        except CoercionFailed:
             return NotImplemented
         else:
             return other.__truediv__(self)
@@ -197,26 +195,26 @@ def __mod__(self, other):
         except CoercionFailed:
             return NotImplemented
         if not other:
-            raise ZeroDivisionError('%s % 0' % self)
+            raise ZeroDivisionError('{} % 0'.format(self))
         else:
-            return self._parent.zero
+            return self._parent.zero  # XXX always 0?
 
-    def _rmod__(self, other):
+    def __rmod__(self, other):
         try:
             other = self._parent.convert(other)
-        except CoercienFailed:
+        except CoercionFailed:
             return NotImplemented
         else:
-            other.__mod__(self)
+            return other.__mod__(self)
 
     def __divmod__(self, other):
         return self.__truediv__(other), self.__mod__(other)
 
-    def _rdivmod__(self, other):
-        return self._rtruediv__(other), self.__rmod__(other)
+    def __rdivmod__(self, other):
+        return self.__rtruediv__(other), self.__rmod__(other)
 
 
-class GaussianDomain(object):
+class GaussianDomain():
     """Base class for Gaussian domains."""
     base = None  # base domain, ZZ or QQ
 
@@ -239,7 +237,7 @@ def from_sympy(self, a):
         if b is I:
             return self.new(x, y)
         else:
-            raise CoercionFailed("%s is not Gaussian" % a)
+            raise CoercionFailed("{} is not Gaussian".format(a))
 
     def convert(self, element):
         """Convert ``element`` to ``self.dtype``.

I then added tests to bring coverage to 92% and only the NotImplemented lines are not covered explicitly.

[smichr]

@oscarbenjamin , since you work with number types a bit, do you think you could sign off on this?

[oscarbenjamin]

I'm still trying to get my head round how all of this stuff works.

Should it not be possible to do something like:

In [3]: Poly(I+x, x, domain=QQ_I)                                                                                                 
---------------------------------------------------------------------------
NameError 

In [4]: Poly(I+x, x, domain='QQ_I')                                                                                               
---------------------------------------------------------------------------
OptionError

It seems I can if I import the domain explicitly but then

In [12]: from sympy.polys.domains.gaussiandomains import QQ_I, ZZ_I                                                               

In [13]: p = Poly(I+x, x, domain=QQ_I)                                                                                            

In [14]: p                                                                                                                        
Out[14]: Poly(x + I, x, domain='QQ_I')

In [15]: p * p                                                                                                                    
---------------------------------------------------------------------------
TypeError 
...
TypeError: unsupported operand type(s) for ** or pow(): 'GaussianRational' and 'int'

[jksuom]

how all of this stuff works.

To make this work automatically with I, it is necessary to add some code to construct_domain().

[smichr]

add some code to construct_domain().

Like "if all the coefficients are integers/rationals, perhaps multiplied by I then use the ZZ_I/QQ_I domain"?

[oscarbenjamin]

Like "if all the coefficients are integers/rationals, perhaps multiplied by I then use the ZZ_I/QQ_I domain"?

Yes. Right now we get

In [1]: Poly(x+I)                                                                                                                 
Out[1]: Poly(x + I, x, I, domain='ZZ')

In [2]: Poly(x+I, x)                                                                                                              
Out[2]: Poly(x + I, x, domain='EX')

In the first example I is treated as a generator. In the second the domain is 'EX' because construct_domain hasn't found a domain that involved I.

It is possible to construct an extension:

In [3]: Poly(x+I, x, extension=True)                                                                                              
Out[3]: Poly(x + I, x, domain='QQ<I>')

That gives a an extension field but it should be possible to have the ring ZZ_I as domain.

Other domains are imported at top-level in sympy as ZZ etc but can also be passed as strings:

In [5]: Poly(x+1, x, domain='QQ')                                                                                                 
Out[5]: Poly(x + 1, x, domain='QQ')

So I guess somewhere there is code to recognise the strings.

Having the domain automatically detected in construct_domain is what will mean that it actually gets used. I think that will then lead to test failures etc and show which parts of the domain need fixing.

[jksuom]

if all the coefficients are integers/rationals, perhaps multiplied by I then use the ZZ_I/QQ_I domain

I think this how it should be. This should come before algebraics are tested so also ZZ_I could be found if the coefficients are integers. Something like this:

for coeff in coeffs:
    is_complex = pure_complex(coeff)
    if is_complex:
        # set rationals if either part is rational
        # nothing if both are integers, and
        gaussians = True
    elif coeff.is_Rational:
        # continue possibly resetting gaussians to False

At the end of _construct_simple something like this:

    else:
        if gaussians:
            if opt.field or rationals:
                domain = QQ_I
            else:
                domain = ZZ_I
        elif opt.field or rationals:
            domain = QQ
        else:
            domain = ZZ

Or possibly:

    else:
        if opt.field or rationals:
            domain = QQ_I if gaussians else QQ
        else:
            domain = ZZ_I if gaussians else ZZ

[czgdp1807]

@smichr @jksuom @oscarbenjamin Is this ready to go?

[jksuom]

There is no support for polynomials yet but it is possible to compute with the domain elements. It think that this could be merged. Polynomial support can be added later.

[oscarbenjamin]

I might try to incorporate the domain as poly and see if I can push that here. I suspect that that would give a more robust test of the implementation here.

[oscarbenjamin]

Actually we can just merge this for now and I'll do that after.

[jksuom]

Thanks!