How to build a modular arithmetic library in Python

Written by Alejandro Sánchez Yalí. Originally published 2022-10-05 on the Monadical blog.

Years ago, I was building my version of the Lights Out Game using Python. For those who are unfamiliar with the original game, it basically consists of a grid of $5$ -by- $5$ lights. When the game starts, a random number or a stored pattern of lights is switched on. Pressing any one of the lights will toggle it and the four adjacent lights. The goal is to turn off all the lights with the fewest number of keystrokes possible.

Figura 1. Lights Out Game. Madsen, 2010. Source.

In this game, a change of state is represented by a change in colour. In the original game, the states only change between two states (two colours). When creating my version of the game, I wanted to see if the states could change between more than two states demonstrating a cyclic sequence of states.

Mathematically, the game can be modeled as a $5$ -by- $5$ matrix (Lights Out: Solutions Using Linear Algebra) where each entry represents the state of a light bulb. Modular arithmetic is an excellent way to find the winning strategy because the state of each light bulb occurs cyclically. Plus, in this case, the usual matrix operations are inefficient due to their unbounded nature. What that means is that it’s possible to get results that are outside the range of the game states, such as decimal numbers or very large numbers that don’t match the ‘on’, ‘off’, or other states in the game.

This is article one of two. In this blog, I'll demonstrate how I have successfully implemented modular arithmetic in Python. In the second post, I'll demonstrate how I used modular arithmetic to model the Lights Out Game's state change as a cyclic sequence.

In doing so, we'll learn how to:

Use matrix algebra in modular arithmetic
Use concepts from graph theory to determine when a Lights Out Game is solvable.

On the programming side, we'll learn how to:

Articulate these mathematical concepts to extend Python's capabilities
Use operator overloading and built-in functions to redefine NumPy’s universal operators.

By the end of this tutorial, we’ll have a Python library with the ability to perform computation using modular arithmetic. This tool can effectively be used in other projects that are related to cryptography, graphs, modeling of cyclic sequence processes, etc.

Problem with NumPy

Before we start our tutorial, I’d like to highlight one of the problems I encountered with NumPy while building my game. NumPy doesn't have support for performing matrix calculations with modular arithmetic, and this was exactly the type of calculation I needed to make. More specifically, I needed to solve the computation of inverse matrices in modular arithmetic in any modulo. In this case, using only the native Python mod operator (%) was not going to be enough to perform this task. Let's see why in the simple example below.

If we have, for example:

import numpy as np
a = np.array([[1, 2], [3, 8]]) % 7
a

>>> array([[1,  2],
       [ 3, 1]])

and calculate the inverse using NumPy, we get:

np.linalg.inv(a)

>>> array([[-0.2,  0.4],
       [ 0.6, -0.2]])

This matrix is not the inverse in modular arithmetic with modulo , mainly because the result is in decimals and not integers. Another way we could try is:

np.linalg.inv(a) % 7

>>> array([[6.8, 0.4],
       [0.6, 6.8]])

Again the result is incorrect because the expected result in modular arithmetic modulo 7 is:

>>> array([[4 6]
       [2 4]])

We encountered this problem because NumPy does not make use of modular arithmetic to perform its internal calculations. We’ll later solve this problem in step 3 of this tutorial by overloading the operator and redefine NumPy’s universal operator.

If you don’t understand modular arithmetic very well, then this problem could be confusing. So let’s take a moment to learn the basics of modular arithmetic in the next section before we dive into our tutorial.

Introduction to modular arithmetic

Before building our library, it’s important to understand what modular arithmetic (also called clock arithmetic) is. If you’re already familiar with modular arithmetic, feel free to skip to the next section.

In mathematics, modular arithmetic is a system of arithmetic for integers that deal primarily with operations and applications regarding remainders. In this system, numbers “cycle” or repeat when reaching a certain value, called the modulo. Let's see what this means by referring to a tool we use every day, a clock!

Let’s look at the $12$ -hour analog clock below. Suppose the clock reads $5$ o’clock. After $3$ hours the clock would read $8$ o’clock. This is determined by simply adding $5 + 3 = 8$ .

Figura 2. Clock arithmetic. University Waterloo, 2019. source.

Now, what time would it be after $12$ hours? After $12$ hours, it would be $5$ o’clock again since the clock cycles back to its original position every $12$ hours. However, if we simply add $5$ and $12$ , we get $5 + 12 = 17$ . But, $17$ isn’t on the clock! What happened?

Figura 3. Clock arithmetic. University Waterloo, 2019. Source.

On an analog clock, the numbers go from $1$ to $12$ , but when we arrive at $13$ o’clock, it appears again as $1$ on the clock. So, $13$ becomes $1$ , $14$ becomes $2$ , $15$ becomes $3$ , and so on. Every time we go past $12$ on the clock, we start counting the hours from $1$ again. Thus, we may view $17$ o’clock as the same as $5$ o’clock. We write this mathematically as:

\begin{equation} 17\cong 5\:(mod\;12) \end{equation}

We use the modular operator mod to indicate that they mean the same thing on a clock. This means that $17$ o’clock is the same as $5$ o’clock in a $12$ hour system. The $(mod\; 12)$ indicates that the clock cycles every $12$ hours.

Similarly, we can add $12$ hours again to $17$ ( $12+17= 29$ ) to get $29$ o’clock. We still understand that it is the same as $5$ o’clock. We write this as:

\begin{equation} 29\cong 5\:(mod\;12) \end{equation}

Figura 4. Clock arithmetic. University Waterloo, 2019. Source.

In general, all those numbers that leave the same remainder as $17$ divided by $12$ , represent the same thing on the clock. In fact, the numbers $... -7, 5, 17, 29, ...$ divided by $12$ leave the remainder $5$ . Therefore they all represent $5$ o’clock. This set $\{... -7, 5, 17, 29, ...\}$ is called the congruence class of $5$ modulo $12$ , which is usually represented as $5\; (mod\; 12)$ . Please note that the $12$ clock hours define $12$ different congruence classes as shown in the following figure:

Figura 5. Clock arithmetic. University Waterloo, 2019. Source.

Mathematically, we will name the hours of the clock as $\mathbb{Z}/12\mathbb{Z}$ and list it like this:

\begin{equation} \mathbb{Z}/12\mathbb{Z} =\{0\;(mod\;12), 1\;(mod\;12), 2\;(mod\;12), ..., 11\;(mod\;12)\} \end{equation}

The elements of this set are all the possible remainders left by the integers when divided by $12$ . Note that instead of listing $12\; (mod\; 12), 0\; (mod\; 12)$ has been listed because $0\cong 12\; (mod\; 12)$ .

There are many other things in our lives that repeat or cycle after a certain amount of time, like days of the week, months in a year, degrees in a circle, or seasons in a year. Modular arithmetic can be used to model any events like this that repeat. If the length of the cycle is $n$ we refer to it as modulo $n$ .

For example, since a clock cycles every $12$ hours, we refer to it as modulo $12$ . In a circle where one full revolution is $360$ degrees, it’s modulo $360$ . Since a week has $7$ days, we refer to it as modulo $7$ .

In general, for a modulo $n$ , all the remainders can be listed as the elements of the set:

\begin{equation} \mathbb{Z}/n\mathbb{Z} =\{0\;(mod\;12), 1\;(mod\;12), 2\;(mod\;12), ..., n-1\;(mod\;n)\} \end{equation}

For the Lights Out Game, each bulb will be modeled as a light that switches between a finite list of colors in a cyclic sequence. Therefore, we can use modular arithmetic $\mathbb{Z}/n\mathbb{Z}$ to tell the computer what state each of the bulbs is in.

As we saw in the clock example (in general in $\mathbb{Z}/n\mathbb{Z}$ ), it’s possible to add hours together if we take into account the cycle of the hours. For example, if we want to calculate $(13 + 15)\; (mod\; 12)$ , we’ll add $13$ and $15$ together to get $28$ . Then, we find that the number $12$ divides into $28$ a total of $2$ times with $4$ left over. We can write this as:

\begin{equation} (13 + 15)\cong 28 \cong 4\;(mod\;12). \end{equation}

Alternatively, we can find the remainder of $13$ and $15$ separately. This makes calculations easier because we are dealing with smaller numbers. Of course, this is even more useful when dealing with bigger numbers. For example, first we note that:

\begin{equation} 13\;(mod\;12)\cong 1 \text{ and } 15 \cong 4\;(mod\;12). \end{equation}

Therefore,

\begin{equation} (13 + 15)\;(mod\;12)\cong 1 + 3 \cong 4\;(mod\;12). \end{equation}

Similarly, we can also do subtraction and multiplication. In general, addition, subtraction and multiplication are defined over $\mathbb{Z}/n\mathbb{Z}$ by the following formulae:

$a\;(mod\;n) \pm b\;(mod\;n) \cong (a\pm b)\;(mod\;n)$
$a\;(mod\;n) \times b\;(mod\;n) \cong (a\times b)\;(mod\;n)$

Moreover, if $n$ is prime, in the set $\mathbb{Z}/n\mathbb{Z}$ , we can calculate the divisions using the Euclidean algorithm.

With this introduction to modular arithmetic, I hope that you have a better understanding of how it works.

Now, let's look at how I used modular arithmetic to build a library that allows me to model the lights in Lights Out as a cyclic sequence of states. Let’s get started!

Tutorial: How to build a modular arithmetic library in Python

For this tutorial we’ll need:

Python 3.10
Familiarity with operator overloading in Python, or read my article on the topic.
A couple of beers, and a lot of patience.

Step 1: Build a Python class to manipulate the elements in modular arithmetic

To start, let’s define a class named Zmodn (integers module $n$ ):

class Zmodn:
    def __init__(self, module: int = 2):
        self.representative = None
        self.module = module

    def __call__(self, integer: int):
        congruence_class = self.__class__(self.module)
        congruence_class.representative = integer % self.module
        return congruence_class

    def __repr__(self):
        return f'{self.representative} (mod {self.module})'

Note that Zmodn receives only the variable module, which allows defining the set of equivalence classes $\mathbb{Z}/n\mathbb{Z}$ according to the value of the module. To generate the elements of $\mathbb{Z}/n\mathbb{Z}$ , we redefine the __call__ method, which allows us to call our class as a function to build different instances of Zmodn by congruence_class = self.__class__(self.module), and calculate the representative (self.representative) of the class as the remainder when dividing integer by self.module. Finally, we redefine the method __repr__ to be represented by the console Zmodn class as a (mod n).

To see how this works, let’s generate the elements of $\mathbb{Z}/n\mathbb{Z} =\{0\; (mod\; 2), 1\; (mod\; 2)\}$ . To do this, we first write:

mod2 = Zmodn(2)

At this point, mod2 is a function that we can pass any integer to. Let’s look at some examples:

mod2(-3)
>>> 1 (mod 2)
mod2(4)
>>> 0 (mod 2)
mod2(5)
>>> 1 mod(2)

As we saw, the only results returned by the function are 0 (mod 2) and 1 (mod 2).

This is because any integer divided by two has a remainder of $0$ or $1$ . We can try other modular classes such as mod3 = Zmodn(3), mod4 = Zmodn(4), etc.

As we saw in my post about operator overloading, it’s possible to redefine all the basic Python operations $+,-,\times,/$ .

The next step is to redefine the methods __add__( $+$ ), __sub__( $-$ ), __mul__( $\times$ ), and __truediv__ ( $/$ ). For this, we’ll make use of what we learned in the operator overloading post. So we have:

def __add__(self, other: int):
    if not self.module == other.module and not isinstance(other, self.__class__):
        raise ValueError(
            f'Cannot add {self} and {other}, they are not in the same module'
        )
    return self.__call__(self.representative + other.representative)

def __sub__(self, other: int):
    if not self.module == other.module and not isinstance(other, self.__class__):
        raise ValueError(
            f'Cannot subtract {self} and {other}, they are not in the same module'
        )
    return self.__call__(self.representative - other.representative)

def __mul__(self, other: int):
    if not self.module == other.module and not isinstance(other, self.__class__):
        raise ValueError(
            f'Cannot multiply {self} and {other}, they are not in the same module'
        )
    return self.__call__(self.representative * other.representative)

Note that before calculating any operation, we must ensure that the elements that we’re adding are calculated with respect to the same module. For example, it doesn’t make sense to add an element of $\mathbb{Z}/2\mathbb{Z}$ with an element of $\mathbb{Z}/3\mathbb{Z}$ . We also verify that other is an instance of Zmodn by using isinstance(other, self.__class__). Then we compute the operations ( $+,-,\times, /)$ with respect to the class representatives to finally return the resulting class.

As for the division between two elements $a$ and $b$ of $\mathbb{Z}/n\mathbb{Z}$ , it only makes sense if $b$ is coprime with $n$ . In that case, we first make use of the Euclidean algorithm for remainders and then calculate the multiplicative inverse of $b$ . So we have:

def multiplicative_inverse(self):
    if self.representative == 0:
        raise ZeroDivisionError('Cannot compute the multiplicative inverse of 0')
    aux1 = 0
    aux2 = 1
    y = self.representative
    x = self.module
    while y != 0:
        q, r = divmod(x, y)
        x, y = y, r
        aux1, aux2 = aux2, aux1 - q * aux2
    if x == 1:
        return self.__call__(aux1 % self.module)
    else:
        raise ValueError(
            f'{self.representative} is not coprime to {self.module}'
            )

And after this, we can define the division as follows:

def __truediv__(self, other: int):
    if not self.module == other.module and not isinstance(other, self.__class__):
        raise ValueError(
            f'Cannot divide {self} and {other}, they are not in the same module'
        )
    return self.__call__(self.representative) * other.multiplicative_inverse()

In the next step, we’ll need the integer representation. So we define the following method:

def __int__(self) -> int:
    return self.representative

Finally, our Zmodn class would be:

class Zmodn:
    def __init__(self, module: int = 2):
        self.representative = None
        self.module = module

    def __call__(self, integer: int):
        congruence_class = self.__class__(self.module)
        congruence_class.representative = integer % self.module
        return congruence_class

    def __int__(self):
        return self.representative

    def __repr__(self):
        return f'{self.representative} (mod {self.module})'

    def __add__(self, other: int):
        if not self.module == other.module and not isinstance(other, self.__class__):
            raise ValueError(
                f'Cannot add {self} and {other}, they are not in the same module'
            )
        return self.__call__(self.representative + other.representative)

    def __sub__(self, other: int):
        if not self.module == other.module and not isinstance(other, self.__class__):
            raise ValueError(
                f'Cannot subtract {self} and {other}, they are not in the same module'
            )
        return self.__call__(self.representative - other.representative)

    def __mul__(self, other: int):
        if not self.module == other.module and not isinstance(other, self.__class__):
            raise ValueError(
                f'Cannot multiply {self} and {other}, they are not in the same module'
            )
        return self.__call__(self.representative * other.representative)

    def multiplicative_inverse(self):
        if self.representative == 0:
            raise ZeroDivisionError('Cannot compute the multiplicative inverse of 0')
        aux1 = 0
        aux2 = 1
        y = self.representative
        x = self.module
        while y != 0:
            q, r = divmod(x, y)
            x, y = y, r
            aux1, aux2 = aux2, aux1 - q * aux2
        if x == 1:
            return self.__call__(aux1 % self.module)
        else:
            raise ValueError(
                f'{self.representative} is not coprime to {self.module}'
                )

    def __truediv__(self, other: int):
        if not self.module == other.module and not isinstance(other, self.__class__):
            raise ValueError(
                f'Cannot divide {self} and {other}, they are not in the same module'
            )
        return self.__call__(self.representative) * other.multiplicative_inverse()

And therefore, we can use it to do our calculations on any set of type $\mathbb{Z}/n\mathbb{Z}$ . Let’s look at some examples:

mod5 = Zmodn(5)
a, b = mod5(7), mod5(9)
a, b
>>> (2 (mod 5), 4 (mod 5))
a + b
>>> 1 (mod 5)
a - b
>>> 3 (mod 5)
a * b
>>> 3 (mod 5)
a / b
>>> 3 (mod 5)
c = a.multiplicative_inverse()
c
>>> 3 (mod 5)
a * c
>>> 1 (mod 5)

Other methods we could implement are:__eq__, __neg__, __isub__, and __iadd__.

Step 2: Create arrays using `Zmodn` and NumPy functionalities

In this step, we’re going to look at how to create arrays using Zmodn with some of NumPy’s functionalities. For this, we’ll build the following class:

import numpy as np
class ZmodnArray:
    def __init__(self, module):
        self.module = module
        self.representatives = None
        self.congruence_class = Zmodn(module)

    def __call__(self, integers: list):
        congruence_class_array = self.__class__(self.module)
        congruence_class = np.vectorize(self.congruence_class)
        congruence_class_array.representatives = np.array(congruence_class(integers))
        return congruence_class_array

    def __repr__(self):
        return f'{self.representatives.astype(int)} (mod {self.module})'

What’s new here is the line congruence_class = np.vectorize(self.congruence_class). In essence, we’re vectorizing the self.congruence_class function using the np.vectorize method. This allows us to pass a complete array to the congruence_class function as done in the line congruence_class_array.representatives = congruence_class(integers).

To complete this class, we add the __add__, __sub__ and __mul__ methods. This has a similar structure to the methods already implemented in Zmodn, the only difference here is that the operations ( $+,-,\times$ ) are executed through NumPy, and we make use of the astype(int) method to obtain the integer representatives of the Zmodn classes. These methods would be:

def __add__(self, other: list):
    if not self.module == other.module and not isinstance(other, self.__class__):
        raise ValueError(
            f'Cannot add {self} and {other}, they are not in the same module'
        )
    return self.__call__((self.representatives + other.representatives).astype(int))

def __sub__(self, other: list):
    if not self.module == other.module and not isinstance(other, self.__class__):
        raise ValueError(
            f'Cannot subtract {self} and {other}, they are not in the same module'
        )
    return self.__call__((self.representatives - other.representatives).astype(int))

def __mul__(self, other: list):
    if not self.module == other.module and not isinstance(other, self.__class__):
        raise ValueError(
            f'Cannot multiply {self} and {other}, they are not in the same module'
        )
    return self.__call__(
        (np.dot(self.representatives, other.representatives)).astype(int)
    )

So far, we can execute some basic operations between arrays of Zmodn elements. Let’s look at some examples:

mod7_array = ZmodnArray(7)
a = mod7_array([[1, 2], [3, 8]])
b = mod7_array([[10, 7], [3, 8]])

The result per console would be:

a
>>> [[1 2]
 [3 1]] (mod 7)
b
>>> [[3 0]
 [3 1]] (mod 7)

For elementary operations we have:

a + b
>>> [[4 2]
 [6 2]] (mod 7)
a - b
>>> [[5 2]
 [0 0]] (mod 7)
a * b
>>> [[2 2]
 [5 1]] (mod 7)

Step 3: Calculate the inverse matrix with modular arithmetic using `ZmodnArray`

Now let’s see how to implement more advanced operations, such as calculating the inverse matrix. A quick way to implement this is to calculate the inverse matrix by the adjugate matrix using the formula for matrix inverse:

\begin{equation} A^{-1}=det(A)^{-1}Adj(A). \end{equation}

For the case of modular arithmetic with matrices in $\mathbb{Z}/n\mathbb{Z}$ , the determinant of a matrix must be coprime with $n$ . We don’t need to get into the technical details of the following functions. For now, what’s important to know is that they’re used to calculate the inverse of a square matrix with non-zero determinant and coprime of $n$ .

Here’s the method for calculating the adjugate matrix:

@staticmethod
def adjoint_of_matrix(matrix):
    adjoint = np.zeros(matrix.shape, dtype=np.int16)
    amount_of_rows, amount_of_columns = matrix.shape
    for i in range(amount_of_rows):
        for j in range(amount_of_columns):
            cofactor_i_j = np.delete(np.delete(matrix, i, axis=0), j, axis=1)
            determinant = int(np.linalg.det(cofactor_i_j))
            adjoint[i][j] = determinant * (-1) ** (i + j)
    return np.transpose(adjoint)

And the method for calculating the inverse matrix:

def inv(self):
    matrix = self.representatives.astype(int)
    if matrix.shape[0] != matrix.shape[1]:
        raise ValueError('Matrix is not square')
    determinant = int(np.linalg.det(matrix))
    if determinant == 0:
        raise ValueError('Matrix is not invertible')
    adjoint = ZmodnArray.adjoint_of_matrix(matrix)
    return self.__call__(
        int(self.congruence_class(1) / self.congruence_class(determinant))
        * adjoint.astype(int)
    )

Although calculating the inverse matrix by the method of the adjugate matrix is not the most efficient, it does allow us to illustrate the possibilities we have with the Zmodn classes. If everything goes well, the result we should obtain by using inv method would be:

a = mod7_array([[1, 2], [3,8]])
a * a.inv()
>>> [[1 0]
 [0 1]] (mod 7)

As a further exercise, I recommend overloading the __getitem__ method because it’ll allow us to read each entry of a Zmodn array as it is done in NumPy arrays. If this still isn’t clear, please take a second look at my post about operator overloading.

Step 4: Redefine NumPy’s universal functions

So far, we’ve seen how to build a Zmodn class that allows us to work with the $\mathbb{Z}/n\mathbb{Z}$ modular classes, and we learned how to use these elements with NumPy.

Next, we’ll solve the problem in a different way by redefining NumPy's universal functions. For this, we’ll use the special method __array_function__, which allows us to override NumPy’s native methods. The structure of this method is:

def __array_ufunc__(self, ufunc, method, *inputs, **kwargs)
	...
	return result

Here:

ufunc parameter is NumPy’s universal function to be called.
method is a string indicating how the ufunc parameter is to be called, either __call__ to indicate that it’s called directly, or one of its methods such as: reduce, accumulate, reduceat, outer or at.
*inputs is a tuple for the ufunc arguments.
**kwargs contains any optional or keyword arguments passed to the function. This includes any output arguments, which are always contained in a tuple.

Therefore, the arguments are normalized; only the necessary input arguments are passed as positional arguments. All others are passed as a dict of keyword arguments (**kwargs). In particular, if there are output arguments (positional or not) that are not None, they are passed as a tuple in the out keyword argument. This goes even for the reduce, accumulate and reduceat methods where all actual cases make sense as a single output.

With this in mind, we can build a class to work with modular arithmetic as follows:

import numpy as np
HANDLED_FUNCTIONS = dict()


class ZmodnArrays:
    def __init__(self, intergers, module):
        self.representatives = np.array(intergers) % module
        self.module = module

    def __repr__(self):
        return f'{self.representatives} (mod {self.module})'

    def __array_function__(self, func, types, args, kwargs):
        if func not in HANDLED_FUNCTIONS:
            return NotImplemented
        if not all(issubclass(t, ZmodnArrays) for t in types):
            return NotImplemented
        return HANDLED_FUNCTIONS[func](*args, **kwargs)

    def implements(numpy_function):
        def decorator(method):
            HANDLED_FUNCTIONS[numpy_function] = method
            return method
        return decorator

Note three important elements:

The HANDLED_FUNCTIONS variable, which is used to store the new definitions for NumPy’s universal functions.
The implementation of the __array_function__method, which is in charge of managing the execution of the definitions.
The implements decorator, which allows us to update the HANDLED_FUNCTIONS variable with the implementation of the new methods.

NumPy has a list of universal functions that can be overridden with this strategy. The list can be found in the documentation. For our case, we’re only going to redefine __add__, __sub__ and __mul__. In effect, the redefinitions would be:

@implements(np.add)
def __add__(self, other):
    if not self.module == other.module and not isinstance(other, self.__class__):
        raise ValueError(
            f'Cannot add {self} and {other}, they are not in the same module'
        )
    repr_sum = np.array(self.representatives) + np.array(other.representatives)
    return self.__class__(repr_sum % self.module, self.module)

@implements(np.subtract)
def __sub__(self, other):
    if not self.module == other.module and not isinstance(other, self.__class__):
        raise ValueError(
            f'Cannot subtract {self} and {other}, they are not in the same module'
        )
    repr_sub = np.array(self.representatives) - np.array(other.representatives)
    return self.__class__(repr_sub % self.module, self.module)

@implements(np.multiply)
def __mul__(self, other):
    if not self.module == other.module and not isinstance(other, self.__class__):
        raise ValueError(
            f'Cannot multiply {self} and {other}, they are not in the same module'
        )
    repr_dot = np.dot(np.array(self.representatives), np.array(other.representatives))
    return self.__class__(repr_dot % self.module, self.module)

As we can see, the structure is similar to what we already implemented in the previous sections. What’s new here is the implementation of the decorator @implements(numpy_function), which is named after the universal function we want to rewrite. Furthermore, in the redefinition, we implement addition, subtraction, and multiplication using NumPy’s native operations, but we’re careful to extract the remainder. For example, in the case of an addition, we return the following:

repr_sum = np.array(self.representatives) + np.array(other.representatives)
    return self.__class__(repr_sum % self.module, self.module)

Note that the sum is calculated through the native sum for objects of type np.array, and the result is that the remainders are computed using %. Then, an instance of the class is made with the method self.__class__. The same is true for the other two basic operations.

Our final class would be:

import numpy as np
HANDLED_FUNCTIONS = dict()


class ZmodnArrays:
    def __init__(self, intergers, module):
        self.representatives = np.array(intergers) % module
        self.module = module

    def __repr__(self):
        return f'{self.representatives} (mod {self.module})'

    def __array_function__(self, func, types, args, kwargs):
        if func not in HANDLED_FUNCTIONS:
            return NotImplemented
        if not all(issubclass(t, ZmodnArrays) for t in types):
            return NotImplemented
        return HANDLED_FUNCTIONS[func](*args, **kwargs)

    def implements(numpy_function):
        def decorator(method):
            HANDLED_FUNCTIONS[numpy_function] = method
            return method
        return decorator

    @implements(np.add)
    def __add__(self, other):
        if not self.module == other.module and not isinstance(other, self.__class__):
            raise ValueError(
                f'Cannot add {self} and {other}, they are not in the same module'
            )
        repr_sum = np.array(self.representatives) + np.array(other.representatives)
        return self.__class__(repr_sum % self.module, self.module)

    @implements(np.subtract)
    def __sub__(self, other):
        if not self.module == other.module and not isinstance(other, self.__class__):
            raise ValueError(
                f'Cannot subtract {self} and {other}, they are not in the same module'
            )
        repr_sub = np.array(self.representatives) - np.array(other.representatives)
        return self.__class__(repr_sub % self.module, self.module)

    @implements(np.multiply)
    def __mul__(self, other):
        if not self.module == other.module and not isinstance(other, self.__class__):
            raise ValueError(
                f'Cannot multiply {self} and {other}, they are not in the same module'
            )
        repr_dot = np.dot(np.array(self.representatives), np.array(other.representatives))
        return self.__class__(repr_dot % self.module, self.module)

To test this class, we can do:

mod7_array = ZmodnArrays(7)
a = mod7_array([[1, 2], [3, 8]])
b = mod7_array([[10, 7], [3, 8]])

So for the elementary operations we would have:

a + b
>>> [[4 2]
 [6 2]] (mod 7)
a - b
>>> [[5 2]
 [0 0]] (mod 7)
a * b
>>> [[2 2]
 [5 1]] (mod 7)

Finally, in the documentation of the data model we can explore what other functionalities we can add to make this class more versatile.

Conclusions

Let’s wrap up what we’ve covered today in our four steps:

The special methods __add__,__sub__, __mul__, and __call__ allow us to redefine the basic Python operators, which allows us to build other types of arithmetic we can work with. This is particularly interesting for problems related to cryptography, graphs, etc.
Using the __array_function__ method, we can tell NumPy how to execute the universal operations according to the needs of our objects.
According to the NumPy documentation, the standard structure for building a class to redefine NumPy’s universal functions is:

HANDLED_FUNCTIONS = {}

class MyArray:
    def __array_function__(self, func, types, args, kwargs):
        if func not in HANDLED_FUNCTIONS:
            return NotImplemented
        # Note: this allows subclasses that don't override
        # __array_function__ to handle MyArray objects
        if not all(issubclass(t, MyArray) for t in types):
            return NotImplemented
        return HANDLED_FUNCTIONS[func](*args, **kwargs)

def implements(numpy_function):
    """Register an __array_function__ implementation for MyArray objects."""
    def decorator(func):
        HANDLED_FUNCTIONS[numpy_function] = func
        return func
    return decorator

@implements(np.concatenate)
def concatenate(arrays, axis=0, out=None):
    ...  # implementation of concatenate for MyArray objects

@implements(np.broadcast_to)
def broadcast_to(array, shape):
    ...  # implementation of broadcast_to for MyArray objects

With this approach of overloading operators and redefining NumPy’s universal functions, we’ve built a modular arithmetic library that can be successfully used in my version of the Lights Out Game. We'll have to meet again in my next post to learn how the library functions in the game.

If you’re like me, you may be interested in building this library to learn more about topics such as graph theory, cryptography, and number theory, among others. To learn more about these topics in relation to the Lights Out Game, please check out these articles:

The code used in this blog can be found in this notebook: Modular Arithmetic with Python.

Finally, if there are any errors, omissions, or inaccuracies in this article, please feel free to contact us through the following Discord channel: Math & Code.