<center> # How to build a modular arithmetic library in Python <big> **We’ll use class operator overloading and redefining NumPy’s built-in functions to solve the Lights Out Game.** </big> *Written by Alejandro Sánchez Yalí. Originally published 2022-10-05 on the [Monadical blog](https://monadical.com/blog.html).* </center> Years ago, I was building my version of the [Lights Out Game](https://www.lightsout.ir/) 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.  <center> Image [source](http://cau.ac.kr/~mhhgtx/courses/LinearAlgebra/references/MadsenLightsOut.pdf): Madsen, 2010. </center> 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](http://cau.ac.kr/~mhhgtx/courses/LinearAlgebra/references/MadsenLightsOut.pdf)) 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: ```python import numpy as np a = np.array([[1, 2], [3, 8]]) % 7 a ``` ```python >>> array([[1, 2], [ 3, 1]]) ``` and calculate the inverse using NumPy, we get: ```python np.linalg.inv(a) ``` ```python >>> array([[-0.2, 0.4], [ 0.6, -0.2]]) ``` This matrix is not the inverse in modular arithmetic with modulo $7$, mainly because the result is in decimals and not integers. Another way we could try is: ```python np.linalg.inv(a) % 7 ``` ```python >>> array([[6.8, 0.4], [0.6, 6.8]]) ``` Again the result is incorrect because the expected result in modular arithmetic modulo $7$ is: ```python >>> 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](https://docs.monadical.com/asanchezyali-modular-arithmetic?view#Tutorial-How-to-build-a-modular-arithmetic-library-in-Python). 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$.  <center> Image [source](https://www.cemc.uwaterloo.ca/events/mathcircles/2019-20/Fall/Junior78_Nov12_Soln.pdf): University Waterloo, 2019. </center> 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?  <center> Image [source](https://www.cemc.uwaterloo.ca/events/mathcircles/2019-20/Fall/Junior78_Nov12_Soln.pdf): University Waterloo, 2019. </center> 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: $$17\cong 5\:(mod\;12)$$ 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: $$29\cong 5\:(mod\;12)$$  <center> Image [source](https://www.cemc.uwaterloo.ca/events/mathcircles/2019-20/Fall/Junior78_Nov12_Soln.pdf): University Waterloo, 2019. </center> 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:  <center> Image [source](https://www.cemc.uwaterloo.ca/events/mathcircles/2019-20/Fall/Junior78_Nov12_Soln.pdf): University Waterloo, 2019. </center> Mathematically, we will name the hours of the clock as $\mathbb{Z}/12\mathbb{Z}$ and list it like this: $$\mathbb{Z}/12\mathbb{Z} =\{0\;(mod\;12), 1\;(mod\;12), 2\;(mod\;12), ..., 11\;(mod\;12)\}$$ 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: $$\mathbb{Z}/n\mathbb{Z} =\{0\;(mod\;12), 1\;(mod\;12), 2\;(mod\;12), ..., n-1\;(mod\;n)\}$$ 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: $$(13 + 15)\cong 28 \cong 4\;(mod\;12).$$ 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: $$13\;(mod\;12)\cong 1 \text{ and } 15 \cong 4\;(mod\;12).$$ Therefore, $$(13 + 15)\;(mod\;12)\cong 1 + 3 \cong 4\;(mod\;12).$$ 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: 1. $a\;(mod\;n) \pm b\;(mod\;n) \cong (a\pm b)\;(mod\;n)$ 2. $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](https://en.wikipedia.org/wiki/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: 1. Python 3.10 2. Familiarity with [operator overloading](https://docs.python.org/3/reference/datamodel.html?highlight=overloading#special-method-names) in Python, or read my [article](https://monadical.com/posts/operator-overloading-in-python.html) on the topic. 3. A couple of beers, and a lot of patience. ### Step 1: Build a Python class to manipulate the elements of $\mathbb{Z}/n\mathbb{Z}$ To start, let’s define a class named `Zmodn` (integers module $n$): ```python 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: ```python mod2 = Zmodn(2) ``` At this point, `mod2` is a function that we can pass any integer to. Let’s look at some examples: ```python 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](https://monadical.com/posts/operator-overloading-in-python.html), 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: ```python 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](https://en.wikipedia.org/wiki/Coprime_integers) with $n$. In that case, we first make use of the [Euclidean algorithm](https://en.wikipedia.org/wiki/Euclidean_algorithm) for remainders and then calculate the multiplicative inverse of $b$. So we have: ```python 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: ```python 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: ```python def __int__(self) -> int: return self.representative ``` Finally, our `Zmodn` class would be: ```python 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: ```python 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: ```python 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](https://numpy.org/doc/stable/reference/generated/numpy.vectorize.html) 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: ```python 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: ```python mod7_array = ZmodnArray(7) a = mod7_array([[1, 2], [3, 8]]) b = mod7_array([[10, 7], [3, 8]]) ``` The result per console would be: ```python a >>> [[1 2] [3 1]] (mod 7) b >>> [[3 0] [3 1]] (mod 7) ``` For elementary operations we have: ```python 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](https://en.wikipedia.org/wiki/Adjugate_matrix) using the formula for matrix inverse: $$A^{-1}=det(A)^{-1}Adj(A).$$ 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: ```python @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: ```python 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: ```python 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](https://monadical.com/posts/operator-overloading-in-python.html). ### 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: ```python 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: ```python 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: 1. the `HANDLED_FUNCTIONS` variable, which is used to store the new definitions for NumPy’s universal functions. 2. the implementation of the `__array_function__method`, which is in charge of managing the execution of the definitions. 3. 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](https://numpy.org/neps/nep-0013-ufunc-overrides.html#list-of-operators-and-numpy-ufuncs). For our case, we’re only going to redefine `__add__`, `__sub__` and `__mul__`. In effect, the redefinitions would be: ```python @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: ```python 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: ```python 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: ```python 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: ```python 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](https://docs.python.org/3/reference/datamodel.html?highlight=overloading#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: 1. 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. 2. Using the `__array_function__` method, we can tell NumPy how to execute the universal operations according to the needs of our objects. 3. According to the NumPy [documentation](https://numpy.org/devdocs/reference/arrays.classes.html#numpy.class.__array_ufunc__), the standard structure for building a class to redefine NumPy’s universal functions is: ```python 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: 1. [A Survey of the Game “Lights Out!”](https://link.springer.com/chapter/10.1007/978-3-642-40273-9_13) 2. [Lights Out: Solutions Using Linear Algebra](http://cau.ac.kr/~mhhgtx/courses/LinearAlgebra/references/MadsenLightsOut.pdf) 3. [Lights Out!: A Survey of Parity Domination in Grid Graphs](https://www.unf.edu/~wkloster/termpaper.pdf) 4. [Lights Out on graphs](https://arxiv.org/pdf/1903.06942.pdf) 5. [Lights Out on a Random Graph](https://journals.calstate.edu/pump/article/view/2593/2983) The code used in this blog can be found in this notebook: [Modular Arithmetic with Python](https://colab.research.google.com/drive/1o27dB5ywv2h0jS9KOfOJV4PsdPAPIi4-?usp=sharing#scrollTo=S7AgP0ydmPrK). If you enjoyed building this modular arithmetic library through my tutorial, please share it with your friends and colleagues! ## References 1. NumPy Documentation: [Standard array subclasses](https://numpy.org/devdocs/reference/arrays.classes.html#numpy.class.__array_ufunc__) 2. NumPy Documentation: [NEP 13 - A mechanism for overriding Ufuncs](https://numpy.org/neps/nep-0013-ufunc-overrides.html) 3. W3schools: [Create your own ufunc](https://www.w3schools.com/python/numpy/numpy_ufunc_create_function.asp) 4. Stackoverflow: [Python linear algebra in a finite field](https://stackoverflow.com/questions/71078626/python-linear-algebra-in-a-finite-field) --- <center> <img src="https://monadical.com/static/logo-black.png" style="height: 80px"/><br/> Monadical.com | Full-Stack Consultancy *We build software that outlasts us* </center>