Solving a system of equations using modular arithmetic

Question

I am trying to implement a solver for the game lights out. You have a grid of lights, when you click on one of them the light you clicked and its four neighbours change colour, with the light starting over when it runs out of colours. The aim is to get all the lights to a particular colour.

Because the way the colours change is cyclic, I thought I could implement it as a system of equations and do all calculations mod n (n being the number of available colours).

This method worked for some puzzles but I got stuck in others.

I am representing the game as a system of equations in an augmented matrix and reducing it to reduced row echelon form using gaussian elimination. As I said in some cases this worked well. However there are some cases (example to follow) where I end up with a line which I cannot reduce completely, reason being that, since I'm using modular arithmetic, some values don't have a multiplicative inverse so I get stuck.

Here is an example: The game shown here represents a 4x4 grid with 4 colours available. Here is the matrix as it started out:

1100100000000000 3
1110010000000000 3
0111001000000000 2
0011000100000000 3
1000110010000000 3
0100111001000000 3
0010011100100000 0
0001001100010000 0
0000100011001000 0
0000010011100100 0
0000001001110010 2
0000000100110001 1
0000000010001100 3
0000000001001110 1
0000000000100111 0
0000000000010011 0

and this is as far as I've managed to reduce it:

1000000000000333 2
0100000000003323 3
0010000000003233 0
0001000000003330 0
0000100000001232 2
0000010000002003 2
0000001000003002 3
0000000100002321 3
0000000010001320 1
0000000001003332 1
0000000000102333 0
0000000000010231 2
0000000000000220 2
0000000000002222 0
0000000000002222 0
0000000000000220 2

In this case the last line, for example, is ...220 2 and I cannot figure out how I can reduce it since I cannot simply divide by 2 (2 has no multiplicative inverse in z4). Whatever I tried, I always ended up with a leading 2 in a row. I am absolutely certain a solution exists for the game (I did solve it correctly myself) but I'm not sure if I'm doing missing something here, or just that the method simply does not work for all cases.

Thanks

Edit: Fixed the matrices in the example. There were inconsistencies at the end. I've figured out that any inconsistencies there mean there is no solution. In the updated example a solution does exist for sure, and this results in no inconsistencies at the end. However I still cannot solve it

Answer 1

Here is a brief sketch of how one can solve linear systems modulo any fixed integer$~n>0$. If $n$ is prime, usual Gaussian reduction of the system works. If $n$ is a product of distinct relatively prime factors$~n_i$, then by the Chinese remainder theorem the system is equivalent to the collection of systems obtained from it by reducing modulo$~n_i$ for every $i$. Supposing each of these systems can be solved (or proven inconsistent, which I consider a form of solving) one can combine the solutions to get the solution for the original system via the Chinese remainder theorem. In detail: if the system reduced modulo any$~n_i$ is inconsistent, so it the original system; otherwise there are parametrised solutions modulo each$~n_i$, one may need to introduce dummy parameters (occurring with coefficient$~0$ everywhere) in some of them to force the same number of free parameters for each$~n_i$, and then the CRT can be applied to a particular solution and to each of the parameters to lift each of them to a value modulo$~n$, which gives a parametrised solution over $\def\Z{\Bbb Z}\Z/n\Z$.

One can thus reduce to the case where each $n_i$ is a prime power; one still needs to bridge the gap between prime powers $p^k$ and the prime$~p$. For this I suggest using a form of Hensel lifting* to the system. First reduce the system modulo$~p$ and solve that. Then if there are any solutions, lift one of them arbitrarily to a non-solution modulo$~p^2$. Then add an unknown multiple of$~p$ as correction term to it, and express the equation that the sum become a true solution modulo$~p^2$; since any defect of the lifted non-solution is a multiple of$~p$ one can divide a factor$~p$ from the equations obtained, and get a linear system modulo$~p$ that can again be solved using Gaussian elimination. Once lifted to$~p^2$, continue similarly to lift to higher powers of$~p$, until reaching$~p^k$.

There is a complication to this idea for which I do not have an answer right now. Solutions will in general not be unique but parametrised; in principle this means one must apply the lifting separately to a particular solution and to the parametrised part (the solution to the corresponding homogeneous system). However it looks like sometimes the question of whether the particular solution can be lifted at all to a higher power of$~p$ may depend on the choice of that particular solution. If this occurs, it would would mean that during the lifting process one would need to allow for adapting the particular solution to avoid inconsistencies in the set of equations obtained. This would be an additional headache, but maybe it can be shown that this never occurs. One thing that certainly can occur is that a solution exists modulo$~p$ but none exist modulo$~p^2$; this unlike the situation of Hensel's lemma (which only deals with one equation, but which is non-linear). A trivial example of that situation is the system of equations $(x\equiv1,x\equiv3)\pmod4$ which is clearly inconsistent but whose reduction modulo$~2$ is not.

*This is actually due to Schönemann, a relatively unknown 19th century mathematician who also first formulated "Eisenstein's criterion" for irreducibility of polynomials.

Answer 2

Here is a probably much more useful answer than my previous one. When solving a system of equations, one can use row operations freely (as usual modifying the right hand side correspondingly). One can however also allow column operations provided one keeps track of them. For instance, one could add twice the first column to the third column, which amounts to replacing the first unknown $x_1$ by a new unknown $x_1'=x_1-2x_3$, keeping all other unknowns $x_i$ (one has $ax_1+cx_3=ax'_1+(2a+c)x_3$ for any $a,c$, which explains the column operation). After solving the new system, one must not forget to recover $x_1=x'_1+2x_3$ from the solution in terms of the new variables. In matrix terms, if we modify the system $A\cdot x=b$ by column operations that amounts to right-multiplication by an invertible matrix $R$, then we get a new system $AR\cdot x'=b$ in terms of the variable $x'=R^{-1}\cdot x$, and at the end one then computes $x=R\cdot x'$. So each column operation applied to$~A$ is accompanied by the same column operation applied to$~R$ (which starts out as the identity), and at the end $R$ is applied (from the left) to the solution vector $x'$ found to obtain the true solution $x$ of the original system.

Given that one can mix row and column operations (recording the latter), one can now arrange that any pivot ultimately used divides all remaining entries (those to it bottom-right), and also the modulus$~n$. The way to achieve this is inspired by Euclid's algorithm, and very similar to the algorithm for reducing a matrix to its Smith normal form. Basically one continually reduces the size of the minimal remaining nonzero entry$~d$ by computing the remainder of the Euclidean division of another remaining entry by$~d$ as long as this is possible (i.e., while $d$ does not divide all others): if that other entry is in the same row or column as$~d$, a single row or column operation will do; if not, one performs a row operation to make a non-multiple of$~d$ appear in the column of$~d$, of which one the forms the remainder by a column operation. As a difference with the Smith normal form algorithm over$~\Bbb Z$, the entries are now all defined modulo$~n$, and so all arithmetic is done modulo$~n$. This in fact allows a significant extra operation to reduce$~d$: whenever $d$ does not divide$~n$, just form a multiple of its row so that$~d$ becomes smaller after reduction modulo$~n$; one easily sees that it can be reduced to $\gcd(d,n)$ in this manner. All in all, even though there are a large number of possibilities to be considered in general (the program is more complicated than for Gaussian elimination), in practice the reduction of$~d$ is very rapid, and most of the time (until the very end of the procedure) one will succeed to reduce it to$~1$.

With this normal form algorithm the system $A\cdot x=b\pmod n$ can be transformed into such a system in which the only nonzero entries of$~A'$ are on the main diagonal, and each one divides its successor. Once in this form, the solution of the system is obvious; to have any solution, each diagonal entry of$~A'$ should divide the corresponding entry of$~b$ (and in particular every zero row of$~A'$ should have the corresponding entry of$~b$ zero.