Finding specific solution for linear system of equations in GF(256)

Question

not sure if this is the right place, but any help is greatly appreciated (noob here).

I am playing around with sagemath and I created the following problem ($A x = y$) in $GF(2^8)$ (which I hoped mirrored byte arithmetic which I see in C f.e). Specifically, let's say I would have given the matrix $A$ and $y$, which $x$ would solve the equation:

F.<a> = GF(2^8, modulus="minimal_weight")
A = Matrix(F, 3, 3)
x = [13, 83, 29]
y = [ 141, 40, 10 ]
m = [
    [23, 56, 2], 
    [92, 234, 182], 
    [12, 94, 132]
]
x = vector([ F.fetch_int(i) for i in x ])
y = vector([ F.fetch_int(i) for i in y ])
for i in range(len(m)):
    A.set_row(i, [ F.fetch_int(i) for i in m[i] ])

First of all, why is it in fact not mirroring the result I would expect from byte arithmetic (i. e. let's say I would implement the same thing in C using unsigned char and specifically ignore overflows):

#include <stdio.h>
int main() {
unsigned char A[3][3] = {
    { 23, 56, 2 },
    { 92, 234, 182 },
    { 12, 94, 132 }
};

unsigned char x[3] = { 13, 83, 29 };

for (int i = 0; i &lt; 3; ++i) {
    unsigned char sum = 0;
    for (int j = 0; j &lt; 3; ++j) {
        sum += x[j] * A[i][j];
    }
    printf(&quot;y[%d]: %d\n&quot;, i, (int)sum);
}

}

Which outputs [ 141, 40, 10 ]. But in sagemath the result is the following:

>>> A * x
[214, 224, 90]

Therefor obviously recovering the solution does not work either:

>>> A.inverse() * y
[194, 215, 241]

Well I was playing around, not understanding a lot, and tried the following:

F = Zmod(2^8)
A = Matrix(F, [
    [23, 56, 2], 
    [92, 234, 182], 
    [12, 94, 132]
])
x = vector(F, [13, 83, 29])
y = vector(F, [141, 40, 10 ])

With this one I at least get the expected result:

>>> A * x
[ 141, 40, 10 ]

.. however the matrix is not invertible anymore:

>>> A.is_invertible()
False

Obviously I have some kind of knowledge gap here and I am not sure where to go next. Thank you for your help.

Answer 1

Part of the answer is that the two questions in the question are in fact different questions.

The ring $\mathbb{Z}/2^8\mathbb{Z}$ of integers modulo $2^8$ is a very different ring from the field $F_{2^8}$ with $2^8$ elements. In particular it is not a field and not even an integral domain.

In Sage these two rings are Zmod(2^8) and GF(2^8).

Pursing that line of thought, given a linear system involving integers between 0 and 255, one can consider it as a system to be solved over $\mathbb{R}$, over $\mathbb{R}$, over $\mathbb{Z}/2^8\mathbb{Z}$, or over the finite field with $2^8$ elements, considering the integers as encodings of field elements.

In this answer we address solving modulo $2^8$ and solving in the finite field.

Solve a linear system modulo 256

Let's see how to define the system, find one solution, find all solutions.

The system

Define the ring $\mathbb{Z}/2^8\mathbb{Z}$, the matrix, the target vector:

sage: Z = Zmod(2^8)
sage: A = matrix(Z, 3, [[23, 56, 2],
....:                   [92, 234, 182],
....:                   [12, 94, 132]])
sage: y = vector(Z, 3, [141, 40, 10])

Check:

sage: A
[ 23  56   2]
[ 92 234 182]
[ 12  94 132]
sage: y
(141, 40, 10)

Find one solution

To find one solution to the system, i.e. a preimage of the target vector under the matrix, use the solve_right method:

sage: x = A.solve_right(y)
sage: x
(13, 83, 157)

Check:

sage: yy = A * x
sage: yy
(141, 40, 10)
sage: yy == y
True

Find all solutions

Finding all solutions is now equivalent to finding the right kernel of $A$ i.e. all vectors $x$ such that $A x$ is zero.

Sadly, finding that kernel is not implemented in Sage:

sage: A.right_kernel()
Traceback (most recent call last)
...
NotImplementedError: Echelon form not implemented
over 'Ring of integers modulo 256'.

The brute force solution is to try all vectors. There are $(2^8)^3$, i.e. $2^{24}$, i.e. 16,777,216 of them. Testing them all would take a while and waste resources.

Instead, we can observe that if $A x$ is zero in $\mathbb{Z}/2^8\mathbb{Z}$, then it is also zero modulo $\mathbb{Z}/2^k\mathbb{Z}$ for $k \in \{7, 6, 5, 4, 3, 2, 1, 0\}$.

For each of these $k$, call $A_k$ the matrix $A$ projected to $\mathbb{Z}/2^k\mathbb{Z}$.

The level $k = 0$ only has the zero matrix and the zero vector, so the kernel there is everything.

For the kernel of the projected matrix $A_1$, we can try all 8 vectors in $\mathbb{Z}/2^1\mathbb{Z}$, we find that 4 of them are part of the kernel.

For each $k$, each vector in the kernel of $A_{k-1}$ has $2^3$ lifts to $(\mathbb{Z}/2^k\mathbb{Z})^3$, which are the only candidates for forming the kernel of $A_k$.

For each starting vector in $\ker(A_1)$, we thus lift progressively, at each step there are eight vectors to try but only one of them ends up being in the kernel at the next level.

We end up computing 8 matrix times vector products in the first level (or computing a kernel basis by finding the row reduced echelon form over $\mathbb{Z}/2\mathbb{Z}$ which is a field, and then, for each of the 4 vectors of the level 1 kernel, 8 matrix times vector products in each of the next 7 levels.

In total, that means $8 + 4 \times 7 \times 8$, i.e. $225$ matrix times vector products. That's a lot more reasonable than the roughly 16.8 million matrix times vector products that brute force search would cost.

This function computes the kernel following that method:

def right_kernel_mod_p_to_the_k(A, p=None, k=None):
    r"""
    Return the kernel of ``A`` modulo ``p^k``
    """
    if p is None or k is None:
        Z = A.base_ring()
        f = Z.cardinality().factor()
        if len(f) == 1:
            p, k = f[0]
        else:
            raise ValueError("either specify p and k or use matrix with " 
                             "base ring some p^k with p prime and k > 0")
    n = A.nrows()
    W = (Zmod(p)^n).list()
    K = A.change_ring(Zmod(p)).right_kernel().list()
    for j in range(2, k + 1):
        R = Zmod(p^j)
        V = R^n
        B = A.change_ring(R)
        W = [p * V(v) for v in W]
        K = [V(x) for x in K]
        # candidates
        C = (u + v for u in K for v in W)
        K = [x for x in C if (B*x).is_zero()]
    return K

Compute the kernel with that function:

sage: K = right_kernel_mod_p_to_the_k(A)
sage: K
[(0, 0, 0), (0, 128, 0), (0, 0, 128), (0, 128, 128)]

Deduce the solutions to the original system:

sage: S = [x + u for u in K]
sage: S
[(13, 83, 157), (13, 211, 157), (13, 83, 29), (13, 211, 29)]

Check:

sage: [A*x for x in S]
[(141, 40, 10), (141, 40, 10), (141, 40, 10), (141, 40, 10)]

Solve a linear system over the finite field with 256 elements

This is a different problem. As given, it amounts to solving a linear system in $F_{2^8}$, with the extra steps of converting field elements from integer representation to polynomial representation and back.

A related question and answer are at

• Ask Sage question 56295: Referring to elements of a polynomial ring as integers

Given the severeal steps, we write a function:

def preimage(y, A):
    r"""
    Return a vector `x` over `F_{2^8}` such that `A x = y`.
Here `F_{2^8}` is the finite field with `2^8` elements
with minimal weight modulus.

INPUT:

- ``y`` -- a vector with entries in `F_{2^8}`
  where each entry is represented as an integer

- ``A`` -- a matrix with entries in `F_{2^8}`
  where each entry is represented as an integer

OUTPUT:

A vector ``x`` with entries in `F_{2^8}` such that `A x = y`.
Each entry is represented as an integer.
"""
F = GF(2^8, 'a', modulus='minimal_weight')
n = A.nrows()
AA = matrix(F, n, n, lambda i, j: F.fetch_int(A[i, j]))
yy = vector(F, n, [F.fetch_int(c) for c in y])
xx = AA.solve_right(yy)
return vector(ZZ, n, [c.integer_representation() for c in xx])

Example:

sage: A = matrix(ZZ, 3, [[23, 56, 2],
....:                    [92, 234, 182],
....:                    [12, 94, 132]])
sage: y = vector(ZZ, 3, [141, 40, 10])
sage: x = preimage(y, A)
sage: x
(194, 215, 241)

Here, finding the kernel would be easier, since the base ring of the relevant matrix is a field.