How to solve a system of 32 XOR equations via Gaussian Elimination?

Question

In the beginning, I needed to solve a system of linear XOR equations. I was inspired by the answer to this question:

how to solve system of linear equations of XOR operation?

, and started solving my problem with Gaussian Elimination. Now I have a system of XOR equations which is represented like this:

0:  0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 | 0
1:  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 | 1
2:  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 | 1
3:  1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 | 1
4:  0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 | 0
5:  0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 | 0
6:  0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 | 0
7:  1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 | 1
8:  0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 | 1
9:  0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 | 0
10: 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 | 1
11: 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 | 0
12: 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 | 1
13: 0 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 | 0
14: 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 | 0
15: 0 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 | 0
16: 0 0 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 | 0
17: 0 0 0 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 | 0
18: 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 | 1
19: 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 | 0
20: 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 | 0
21: 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 | 1
22: 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 | 0
23: 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 | 0
24: 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 | 1
25: 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 | 1
26: 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 | 1
27: 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 | 0
28: 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 | 0
29: 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0 | 1
30: 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 0 | 1
31: 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 | 0

Every column is the variable (32 in total), and every row is another equation (32 in total too). So, since the amount of equations is equal to the amount of variables, it's 100% solvable. I'm just not experienced with Gaussian Elimination, and maybe you could give me at least basic ideas on how to solve that amount of equations.

Disclaimer

This question is part of my original question I asked on Stack Overflow, in case you're wondering why do I even need this in the first place:

https://stackoverflow.com/questions/66607696/reverse-sha-256-sigma0-function-within-complexity-of-on

Answer 1

All calculations are in $GF(2).$ I use $+$ instead of $\oplus$ for "xor".

You can make use of the special structure of this system of equations. If $R$ is the $32\times 32$ matrix that describes the circular shift: $$ R = \begin{pmatrix} 0 & & & \cdots & & 0 & 1 \\ 1 & 0 & & & & & 0 \\ 0 & 1 & 0 & & & & \\ & 0 & 1 & \ddots & & & \vdots \\ \vdots & & 0 & \ddots & 0 & & \\ & & & \ddots & 1 & 0 & \\ 0 & & \cdots & & 0 & 1 & 0 \end{pmatrix} $$ then your matrix $A$ can be written as $A=R^3+R^7+R^{18}.$ Futhermore, we know that $R^{32}+I=0.$ You can apply the Extended Euclidean Algorithm to the polynomials $x^{18}+x^{7}+x^{3}$ and $x^{32}+1$ to find two polynomials $p$ and $q$, such that $$ \left(x^{18}+x^{7}+x^{3}\right)p(x) + \left(x^{32}+1\right)q(x)=1 $$ Then $$ \left(R^{18}+R^{7}+R^{3}\right)p(R) + \left(R^{32}+I\right)q(R)=I $$ or, as $R^{32}+I=0$, $$ A\,p(R) = I $$ which means $A^{-1} = p(R).$ Therefore, the solution of $Ax=b$ can be obtained by using $x=p(R)b.$

Note that the matrix $p(R)$ does not have to be explicitly calculated. $$ x=p(R)b = \left(R^{a_0}+R^{a_1}+R^{a_2}+\ldots \right)b = R^{a_0}b+R^{a_1}b+R^{a_2}b+\ldots $$ As $R$ is the circular shift, this means: Initialize your result with $0$. Then take $b$ and rotate it by $a_0$ bits and add/xor this to the result, then take $b$ and rotate it by $a_1$ bits and add/xor this to the result etc.

Edit

There is also a way to get $p$ without the Extended Euclidean Algorithm. In $GF(2)$, we have $(x+y)^{2^n}=x^{2^n}+y^{2^n},$ which means $$ \left(R^{18}+R^{7}+R^{3}\right)^{32} = R^{18\cdot 32}+R^{7\cdot 32}+R^{3\cdot 32} = I+I+I = I $$ If we use exponentiation by squaring and always reduce exponents as much as possible due to $R^{32}=I,$ we get $$ \left(R^{18}+R^{7}+R^{3}\right)^{-1} = \left(R^{18}+R^{7}+R^{3}\right)^{31} \\ =\left(R^{18}+R^{7}+R^{3}\right)^{16}\left(R^{18}+R^{7}+R^{3}\right)^{8}\left(R^{18}+R^{7}+R^{3}\right)^{4}\left(R^{18}+R^{7}+R^{3}\right)^{2}\left(R^{18}+R^{7}+R^{3}\right)^{1} \\ =(I)(R^{16})(R^{8}+R^{28}+R^{12})(R^{4}+R^{14}+R^{6})(R^{18}+R^7+R^3) $$ This product does not have to be expanded. It can easily be used directly to calculate $x.$