How to find the inverse of a 3x3 MDS matrix

Question

I implemented a block cipher similar to AES. But the reason I can't decrypt is that I can't get the inverse MDS matrix. The MDS matrix I used is a 3x3 MDS matrix on $GF(2^8) \implies GF(2^8)$ like AES

\begin{bmatrix} 1 & 2 &2 \\ 2 & 2 & 1\\ 2 & 1 & 2\\ \end{bmatrix}

The encryption process is the same as that of AES mixcolumns, and the MDS matrix used is as above.

I need the MDS inverse matrix required for decryption. What is the inverse of this MDS matrix on $GF(2^8)$?

I searched for Euclid-Wallis Algorithm or extended euclidean algorithm but I did not understand it well.

How can I find the inverse of a 3x3 MDS matrix?

kelalaka · Answer 1 · 2020-09-17T18:54:31.670

Your matrix contains 2 which is not an element of $\operatorname{GF}(2^8)$ unless it means $x$. In this case, we can use SageMath to find the inverse as

R.<x> = PolynomialRing(GF(2), 'x')
S.<y> = QuotientRing(R, R.ideal(x^8+x^4+x^3+x+1))
S.is_field()
S.cardinality()
y^8 + y^4 + y^3 + y + 1
A = matrix(S,[[1,y,y],[y,y,1],[y,1,y],])
A.inverse()

That produces this output

True
256
0
[                                  1 y^7 + y^6 + y^5 + y^4 + y^2 + y + 1 y^7 + y^6 + y^5 + y^4 + y^2 + y + 1]
[y^7 + y^6 + y^5 + y^4 + y^2 + y + 1 y^7 + y^6 + y^5 + y^4 + y^2 + y + 1                                   1]
[y^7 + y^6 + y^5 + y^4 + y^2 + y + 1                                   1 y^7 + y^6 + y^5 + y^4 + y^2 + y + 1]

The matrix written in standard binary representation is:

$$\begin{bmatrix} 00000001 & 11110111 & 11110111 \\ 11110111 & 11110111 & 00000001 \\ 11110111 & 00000001 & 11110111 \end{bmatrix}$$

Stack Danny · Answer 2 · 2022-10-27T13:46:25.580

This is how I was able to compute this result.

I followed the steps on the answer on how to invert a matrix in modular arithmetic.
All operations are replaced with their $GF(2^8)$ counterpart as described by this answer on "calculating the inverse of a polynomial".

I want to show the C++ code I was able to write with these two answers. With the code we can invert any N sized matrix in the $GF(2^8)$ field. The full code can be found here.

These are the galois operations as shown by Thomas Pornin:

#include <iostream>
#include <array>
namespace galois {
    constexpr int add(int x, int y) {
        return x ^ y;
    }
constexpr int mul(int x, int y) {
    int z = 0;
    for (int i = 0; i &lt; 8; ++i) {
        z ^= x &amp; -(y &amp; 1);
        y &gt;&gt;= 1;
        x &lt;&lt;= 1;
        x ^= (0x11B &amp; -(x &gt;&gt; 8));
    }
    return z;
}

constexpr int inverse(int x) {
    int z = x;
    for (int i = 0; i &lt; 6; ++i) {
        z = galois::mul(z, z);
        z = galois::mul(z, x);
    }
    return galois::mul(z, z);
}

constexpr int divide(int x, int y) {
    return mul(x, inverse(y));
}

}

And this is the matrix_inverse as described by Jyrki Lahtonen, but with using the galois operations:

template<std::size_t N>
using matrix = std::array<std::array<int, N>, N>;
template<std::size_t N>
matrix<N> matrix_inverse(matrix<N> M) {
matrix&lt;N&gt; R{}; //fills R with 0s
for (unsigned i = 0; i &lt; N; ++i) {
    R[i][i] = 1; //makes the matrix the identity matrix
}
for (unsigned i = 0; i &lt; N; ++i) {

    //sets M[i][i] to 1, divide by row by inverse
    auto diagonal = M[i][i];
    auto inv = galois::inverse(diagonal);
    for (unsigned row = 0; row &lt; N; ++row) {
        M[i][row] = galois::mul(inv, M[i][row]);
        R[i][row] = galois::mul(inv, R[i][row]);
    }

    //pivots the column
    for (unsigned col = 0; col &lt; N; ++col) {
        if (col == i) continue;

        auto n = M[col][i];
        for (unsigned row = 0; row &lt; N; ++row) {
            M[col][row] = galois::add(galois::mul(n, M[i][row]), M[col][row]);
            R[col][row] = galois::add(galois::mul(n, R[i][row]), R[col][row]);
        }
    }
}
return R;

}

And this is how we can call matrix_inverse to test the result:

template<std::size_t N>
void print_matrix(matrix<N> m) {
    for (auto& i : m) {
        for (auto& i : i) {
            std::cout << i << ' ';
        }
        std::cout << '\n';
    }
}
int main() {
    matrix<3> m = 
    { {
        {1, 2, 2},
        {2, 2, 1},
        {2, 1, 2},
    } };
std::cout &lt;&lt; &quot;Input -\n&quot;;
print_matrix(m);

matrix&lt;3&gt; inverse = matrix_inverse(m);

std::cout &lt;&lt; &quot;\nInverse -\n&quot;;
print_matrix(inverse);

}

As can be seen in the output, we calculated the correct inverse of your matrix:

Input -
1 2 2
2 2 1
2 1 2
Inverse -
1 247 247
247 247 1
247 1 247

As a second test, we can replicate the correct AES decryption MixColumns matrix:

matrix<4> m =
{ {
    {2, 3, 1, 1},
    {1, 2, 3, 1},
    {1, 1, 2, 3},
    {3, 1, 1, 2},
} };
matrix<4> inverse = matrix_inverse(m);
std::cout << "Inverse -\n";
print_matrix(inverse);

Output:

How to find the inverse of a 3x3 MDS matrix

2 Answers2

Linked