Complex numbers support for modular arithmetic

Ideas
1

Seeing how Python has a support for a lot of different mathematical operations for complex numbers, I’m surprised that there is no support for modular arithmetic with complex numbers, that is z % w for z,w\in\mathbb C.

Sure, modular arithmetic for complex numbers is a bit obscure of a topic, however I don’t really see why it shouldn’t be implemented given how simple it is.

Actually I have my own implementation of it using a class method:

# Minimum reproducible example. Here I use `i` as the complex unit
# and also only use 2 decimal places of precision, which seems to
# be enough to get correct values.

class Complex(object):
    def __init__(self, real, imag):
        self.r = real
        self.i = imag

    def __sub__(self, num):
        real = self.r - num.r
        imag = self.i - num.i
        return Complex(real, imag)

    def __mul__(self, num):
        real = self.r * num.r - self.i * num.i
        imag = self.r * num.i + self.i * num.r
        return Complex(real, imag)

    def __truediv__(self,num):
        try:
            x, y, z = self.r * num.r + self.i * num.i, \
                      self.i * num.r - self.r * num.i, \
                      num.r ** 2 + num.i ** 2
            
            real = x / z
            imag = y / z
            return Complex(real, imag)
       
        except ZeroDivisionError as ZDE:
            return ZDE

    def __round__(self):
        real = (self.r + 0.5) - ((self.r + 0.5) % 1.0)
        imag = (self.i + 0.5) - ((self.i + 0.5) % 1.0)
        return Complex(real, imag)
    
    def __mod__(self, modu):
        # We can define for two complex numbers `z = a + b*i` and `w = c + d*i`
        # z % w = z - w*round(z / w). Note that the correct formula actually
        # uses `floor` and not `round` but it doesn't seem to work correctly
        # if I change the method to `floor` here in Python.
        return self - ((self / modu).__round__() * modu)

    def __str__(self):
        if self.i >= 0:
            res = "%.2f+%.2fi" % (self.r, self.i)
        else:
            res = "%.2f-%.2fi" % (self.r, abs(self.i))
        return res

# Examples: (verified with Wolfram Alpha)
# print(Complex(0.06, 9.03)   % Complex(1.02, 13.01)) # -0.96-3.98i
# print(Complex(1, 2)         % Complex(1, 3))        #  0.00-1.00i
# print(Complex(5107, 607)    % Complex(23, 27))      # -4.00+22.00i
# print(Complex(503.03, 6.56) % Complex(23.1,27.09))  # -2.59+16.85i
# print(Complex(4, 0)         % Complex(0, 1))        #  0.00+0.00i

So it seems like it could be implemented easily enough, but if this hasn’t already been implemented in Python then there must be some other reason why it currently isn’t, and that I wonder why.

2

What do you use this for?

3

Can you give any references that remainder is a well defined operation between complex numbers? I can’t really make sense of it and I am pretty sure your definition is invalid.

Also note that % should form a pair with //, i.e. floordiv, making up the function divmod, so you also need to define that one.

And e.g. the floor function doesn’t work on complex numbers for similar reasons.

4

I’m aware that remainder isn’t a well defined operation on \mathbb C like it is for integers.

The approach that I’m pretty sure Wolfram Alpha uses is to make sure that we return the remainder with the least absolute value, which is the standard approach with Gaussian integers.

Relevant Mathematics Stack Exchange link

5

Ok, so your proposal is actually to add semantic support for something like Gaussian integers into python and define at least the following functions for this:

  • floor
  • ceil
  • round
  • //
  • %
  • divmod

I am not opposed, but I do think this would need a PEP that justifies the usefulness of this.

Just filling in operations that are currently undefined is not a good enough reason. It is sometimes useful to get an error message if you are passing in incompatible types - getting potentially confusing results is not necessarily good.

Also note that floor, ceil, round are normally defined by comparison operations - which are not defined for complex numbers either (and can’t really be defined)

2 Likes
6

They’re defined that way? The way I defined round (as you can see in my class method) was by taking the choice round(a + b*i) = round(a) + round(b)*i and then doing the same things for ceil and floor.

Aside from having to check if the real part of the complex number is equivalent mod 1.0 to 0.0 (and then the same thing for the imaginary part) for ceil, I never used any comparison functions to define them.

7

Not in code, but in semantics. floor is defined as “the largest natural number n such that n <= x”, which involves two different assumptions about the order-ability of real numbers. This doesn’t really extend to complex numbers, you have to redefine what you mean with floor.

1 Like
9

Neat. But obscure (as you say yourself).

If none of numpy, scipy, sympy, and sagemath etc. implements this already, put it in your own third party library on PyPi and see if many people use it.

[Edit] Python 2 supported this.

Python 2.7.18 (v2.7.18:8d21aa21f2, Apr 20 2020, 13:25:05) [MSC v.1500 64 bit (AMD64)] on win32
Type "help", "copyright", "credits" or "license" for more information.
>>> 3+6j % 4+5j
(3+11j)
>>> (3+6j) % (4+5j)
(-1+1j)

Complex Numbers Modulus in Python - Stack Overflow I’m guessing it was intentionally deprecated.

static PyObject *
complex_divmod(PyObject *v, PyObject *w)
{
    Py_complex div, mod;
    PyObject *d, *m, *z;
    Py_complex a, b;
    TO_COMPLEX(v, a);
    TO_COMPLEX(w, b);
    if (PyErr_Warn(PyExc_DeprecationWarning,
                   "complex divmod(), // and % are deprecated") < 0)
        return NULL;

    errno = 0;
    div = c_quot(a, b); /* The raw divisor value. */
    if (errno == EDOM) {
        PyErr_SetString(PyExc_ZeroDivisionError, "complex divmod()");
        return NULL;
    }
    div.real = floor(div.real); /* Use the floor of the real part. */
    div.imag = 0.0;
    mod = c_diff(a, c_prod(b, div));
    d = PyComplex_FromCComplex(div);
    m = PyComplex_FromCComplex(mod);
    z = PyTuple_Pack(2, d, m);
    Py_XDECREF(d);
    Py_XDECREF(m);
    return z;
}
2 Likes
10

SymPy has Gaussian integers and Gaussian rationals:

>>> from sympy import ZZ_I, QQ_I
>>> divmod(ZZ_I(1,2), ZZ_I(1,1))
(ZZ_I(2, 1), ZZ_I(0, -1))
>>> divmod(QQ_I(1,2), QQ_I(1,1))
(QQ_I(3/2, 1/2), QQ_I(0, 0))
>>> ZZ_I.gcd(ZZ_I(-1, 5), ZZ_I(2))
ZZ_I(1, 1)

For the Gaussian integers there is a conventional notion of Euclidean division that can be defined to make the quotient and remainder unique. This definition can be used e.g. to compute the greatest common divisor of two Gaussian integers

The Gaussian rationals are a field so an exact quotient always exists and as implemented here in SymPy the remainder is always zero. It could have been possible to define that differently though e.g. to bring the two Gaussian rationals over a common integer denominator and then use Gaussian integer Euclidean division.

Typically field domains in SymPy just don’t define divmod etc so when using domain elements directly you need to know if you are working in a field or in a Euclidean domain or an integral domain etc. Generally the operation depends on which ring you are working in rather than being well defined for arbitrary complex numbers.

1 Like
11

Thanks - that’s a useful insight into why SymPy is structured the way it is.

1 Like