[Update]

Your calculation are individually correct. However, the polyonmial p4 you get at the end is almost the modular inverse you are looking for.

The steps of the extended Eulclid algorithm are: $$ \begin{array}{rcccc} p & = & 1 \times p & + & 0 \times a\\ a & = & 0 \times p & + & 1 \times a \\ r_0 & = & 1\times p & + &q_0 \times a \\ r_1 & = & q_1 \times p & + &(q_0q_1 + 1) \times a \\ r_2 & = & (q_1q_2 + 1)\times p & + & (q_2(q_0q_1 + 1) + q_0)\times a \end{array} $$ and the coefficients in front of $a$ are the polynomials $p_0$, $p_1$, $p_2$, $p_3$ and $p_4$ you computed. As you will see, the last line says that $$ p_4\times a \equiv r_2 \mod p, $$ so the inverse of $a$ is indeed $p_4 \times r_2^{-1}$ and here the value $r_2$ is {9a}.

You are only one modular inverse in $GF(2^8)$ away from finishing your calculation.

I will present an alternative method to find the inverse of the polynomial.

Let $p(x) = ax^3 + bx^2 + cx + d$ a polynomial of degree $3$ in the polynomial ring of the finite field $GF(2^8)$. We want to find $q(x) = \alpha x^3 + \beta x^2 + \gamma x + \delta$ such that $p(x)q(x) \equiv 1 \bmod x^4 + 1$.

We compute the product $p(x)q(x)$: $$ \begin{array}{rcl} p(x)q(x) & = & a\alpha x^6 + (a\beta + \alpha b) x^5 + (a\gamma + b\beta + c\alpha) x^4 + \\ & & (a\delta + b\gamma + c\beta + d\alpha) x^3 + (b\delta + c\gamma + d\beta) x^2 +\\ & & (c\delta + d\gamma) x + d\delta. \end{array} $$ But we want the product mod $x^4 + 1$, and we have $x^4 \equiv -1 \bmod x^4 + 1$, and even better since we are in a field of characteristic two, we have $x^4 \equiv 1 \bmod x^4 + 1$, so $x^5 \equiv x \bmod x^4 + 1$ and $x^6 \equiv x^2 \bmod x^4 + 1$.

Therefore we have $$ \begin{array}{rcl} p(x)q(x) & \equiv & (a\delta + b\gamma + c\beta + d\alpha) x^3 +\\ & & (b\delta + c\gamma + d\beta + a\alpha) x^2 + \\ & & (c\delta + d\gamma + a\beta + b\alpha) x + \\ & & (d\delta + a\gamma + b\beta + c\alpha) \end{array}\mod x^4 + 1 $$ Since we want $p(x)q(x) \equiv 1 \bmod x^4 + 1$, we have to solve a system of linear equations: $$ \left\{\begin{array}{rcl} a\delta + b\gamma + c\beta + d\alpha & = & 0 \\ b\delta + c\gamma + d\beta + a\alpha & = & 0 \\ c\delta + d\gamma + a\beta + b\alpha & = & 0 \\ d\delta + a\gamma + b\beta + c\alpha & = & 1, \end{array}\right. $$ which can be rewritten as $$ \begin{bmatrix} a & b & c & d \\ b & c & d & a \\ c & d & a & b \\ d & a & b & c \end{bmatrix}\cdot \begin{bmatrix}\delta\\ \gamma \\ \beta \\ \alpha\end{bmatrix} = \begin{bmatrix}0\\0\\0\\1\end{bmatrix} $$ To find the coefficients $\alpha$, $\beta$, $\gamma$ and $\delta$ of the polynomial, we have only to find the inverse of the matrix: $$ \begin{bmatrix}\delta\\ \gamma \\ \beta \\ \alpha\end{bmatrix} = \begin{bmatrix} a & b & c & d \\ b & c & d & a \\ c & d & a & b \\ d & a & b & c \end{bmatrix}^{-1}\cdot\begin{bmatrix}0\\0\\0\\1\end{bmatrix} $$ In fact, the coefficients will be the last column of this matrix.

You can compute the inverse with a method such as Gauss elimination, where all computations are in the field $GF(2^8)$.

In this specific case, the matrix keeping your notation) is: $$ \begin{bmatrix} 03 & 01 & 01 & 02 \\ 01 & 01 & 02 & 03 \\ 01 & 02 & 03 & 01 \\ 02 & 03 & 01 & 01 \end{bmatrix} $$

Whichever method you use, I hope you will get through all the calculation.

[Update]

Your calculation are individually correct. However, the polyonmial p4 you get at the end is almost the modular inverse you are looking for.

The steps of the extended Eulclid algorithm are: $$ \begin{array}{rcccc} p & = & 1 \times p & + & 0 \times a\\ a & = & 0 \times p & + & 1 \times a \\ r_0 & = & 1\times p & + &q_0 \times a \\ r_1 & = & q_1 \times p & + &(q_0q_1 + 1) \times a \\ r_2 & = & (q_1q_2 + 1)\times p & + & (q_2(q_0q_1 + 1) + q_0)\times a \end{array} $$ and the coefficients in front of $a$ are the polynomials $p_0$, $p_1$, $p_2$, $p_3$ and $p_4$ you computed. As you will see, the last line says that $$ p_4\times a \equiv r_2 \bmod p, $$ so the inverse of $a$ is indeed $p_4 \times r_2^{-1}$ and here the value $r_2$ is {9a}.

You are only one modular inverse in $GF(2^8)$ away from finishing your calculation.

I will present an alternative method to find the inverse of the polynomial.

Let $p(x) = ax^3 + bx^2 + cx + d$ a polynomial of degree $3$ in the polynomial ring of the finite field $GF(2^8)$. We want to find $q(x) = \alpha x^3 + \beta x^2 + \gamma x + \delta$ such that $p(x)q(x) \equiv 1 \bmod x^4 + 1$.

We compute the product $p(x)q(x)$: $$ \begin{array}{rcl} p(x)q(x) & = & a\alpha x^6 + (a\beta + \alpha b) x^5 + (a\gamma + b\beta + c\alpha) x^4 + \\ & & (a\delta + b\gamma + c\beta + d\alpha) x^3 + (b\delta + c\gamma + d\beta) x^2 +\\ & & (c\delta + d\gamma) x + d\delta. \end{array} $$ But we want the product mod $x^4 + 1$, and we have $x^4 \equiv -1 \bmod x^4 + 1$, and even better since we are in a field of characteristic two, we have $x^4 \equiv 1 \bmod x^4 + 1$, so $x^5 \equiv x \bmod x^4 + 1$ and $x^6 \equiv x^2 \bmod x^4 + 1$.

Therefore we have $$ \begin{array}{rcl} p(x)q(x) & \equiv & (a\delta + b\gamma + c\beta + d\alpha) x^3 +\\ & & (b\delta + c\gamma + d\beta + a\alpha) x^2 + \\ & & (c\delta + d\gamma + a\beta + b\alpha) x + \\ & & (d\delta + a\gamma + b\beta + c\alpha) \end{array}\mod x^4 + 1 $$ Since we want $p(x)q(x) \equiv 1 \bmod x^4 + 1$, we have to solve a system of linear equations: $$ \left\{\begin{array}{rcl} a\delta + b\gamma + c\beta + d\alpha & = & 0 \\ b\delta + c\gamma + d\beta + a\alpha & = & 0 \\ c\delta + d\gamma + a\beta + b\alpha & = & 0 \\ d\delta + a\gamma + b\beta + c\alpha & = & 1, \end{array}\right. $$ which can be rewritten as $$ \begin{bmatrix} a & b & c & d \\ b & c & d & a \\ c & d & a & b \\ d & a & b & c \end{bmatrix}\cdot \begin{bmatrix}\delta\\ \gamma \\ \beta \\ \alpha\end{bmatrix} = \begin{bmatrix}0\\0\\0\\1\end{bmatrix} $$ To find the coefficients $\alpha$, $\beta$, $\gamma$ and $\delta$ of the polynomial, we have only to find the inverse of the matrix: $$ \begin{bmatrix}\delta\\ \gamma \\ \beta \\ \alpha\end{bmatrix} = \begin{bmatrix} a & b & c & d \\ b & c & d & a \\ c & d & a & b \\ d & a & b & c \end{bmatrix}^{-1}\cdot\begin{bmatrix}0\\0\\0\\1\end{bmatrix} $$ In fact, the coefficients will be the last column of this matrix.

You can compute the inverse with a method such as Gauss elimination, where all computations are in the field $GF(2^8)$.

In this specific case, the matrix keeping your notation) is: $$ \begin{bmatrix} 03 & 01 & 01 & 02 \\ 01 & 01 & 02 & 03 \\ 01 & 02 & 03 & 01 \\ 02 & 03 & 01 & 01 \end{bmatrix} $$

Whichever method you use, I hope you will get through all the calculation.

I pointed out a mistake in your calculation in a comment, so maybe the remaining of your calculations can be made correct and give the right answer. However, I will present an alternative method to find the inverse of the polynomial.

Let $p(x) = ax^3 + bx^2 + cx + d$ a polynomial of degree $3$ in the polynomial ring of the finite field $GF(2^8)$. We want to find $q(x) = \alpha x^3 + \beta x^2 + \gamma x + \delta$ such that $p(x)q(x) \equiv 1 \bmod x^4 + 1$.

We compute the product $p(x)q(x)$: $$ \begin{array}{rcl} p(x)q(x) & = & a\alpha x^6 + (a\beta + \alpha b) x^5 + (a\gamma + b\beta + c\alpha) x^4 + \\ & & (a\delta + b\gamma + c\beta + d\alpha) x^3 + (b\delta + c\gamma + d\beta) x^2 +\\ & & (c\delta + d\gamma) x + d\delta. \end{array} $$ But we want the product mod $x^4 + 1$, and we have $x^4 \equiv -1 \bmod x^4 + 1$, and even better since we are in a field of characteristic two, we have $x^4 \equiv 1 \bmod x^4 + 1$, so $x^5 \equiv x \bmod x^4 + 1$ and $x^6 \equiv x^2 \bmod x^4 + 1$.

Therefore we have $$ \begin{array}{rcl} p(x)q(x) & \equiv & (a\delta + b\gamma + c\beta + d\alpha) x^3 +\\ & & (b\delta + c\gamma + d\beta + a\alpha) x^2 + \\ & & (c\delta + d\gamma + a\beta + b\alpha) x + \\ & & (d\delta + a\gamma + b\beta + c\alpha) \end{array}\mod x^4 + 1 $$ Since we want $p(x)q(x) \equiv 1 \bmod x^4 + 1$, we have to solve a system of linear equations: $$ \left\{\begin{array}{rcl} a\delta + b\gamma + c\beta + d\alpha & = & 0 \\ b\delta + c\gamma + d\beta + a\alpha & = & 0 \\ c\delta + d\gamma + a\beta + b\alpha & = & 0 \\ d\delta + a\gamma + b\beta + c\alpha & = & 1, \end{array}\right. $$ which can be rewritten as $$ \begin{bmatrix} a & b & c & d \\ b & c & d & a \\ c & d & a & b \\ d & a & b & c \end{bmatrix}\cdot \begin{bmatrix}\delta\\ \gamma \\ \beta \\ \alpha\end{bmatrix} = \begin{bmatrix}0\\0\\0\\1\end{bmatrix} $$ To find the coefficients $\alpha$, $\beta$, $\gamma$ and $\delta$ of the polynomial, we have only to find the inverse of the matrix: $$ \begin{bmatrix}\delta\\ \gamma \\ \beta \\ \alpha\end{bmatrix} = \begin{bmatrix} a & b & c & d \\ b & c & d & a \\ c & d & a & b \\ d & a & b & c \end{bmatrix}^{-1}\cdot\begin{bmatrix}0\\0\\0\\1\end{bmatrix} $$ In fact, the coefficients will be the last column of this matrix.

You can compute the inverse with a method such as Gauss elimination, where all computations are in the field $GF(2^8)$.

In this specific case, the matrix keeping your notation) is: $$ \begin{bmatrix} 03 & 01 & 01 & 02 \\ 01 & 01 & 02 & 03 \\ 01 & 02 & 03 & 01 \\ 02 & 03 & 01 & 01 \end{bmatrix} $$

Whichever method you use, I hope you will get through all the calculation.

[Update]

Your calculation are individually correct. However, the polyonmial p4 you get at the end is almost the modular inverse you are looking for.

The steps of the extended Eulclid algorithm are: $$ \begin{array}{rcccc} p & = & 1 \times p & + & 0 \times a\\ a & = & 0 \times p & + & 1 \times a \\ r_0 & = & 1\times p & + &q_0 \times a \\ r_1 & = & q_1 \times p & + &(q_0q_1 + 1) \times a \\ r_2 & = & (q_1q_2 + 1)\times p & + & (q_2(q_0q_1 + 1) + q_0)\times a \end{array} $$ and the coefficients in front of $a$ are the polynomials $p_0$, $p_1$, $p_2$, $p_3$ and $p_4$ you computed. As you will see, the last line says that $$ p_4\times a \equiv r_2 \mod p, $$ so the inverse of $a$ is indeed $p_4 \times r_2^{-1}$ and here the value $r_2$ is {9a}.

You are only one modular inverse in $GF(2^8)$ away from finishing your calculation.

I will present an alternative method to find the inverse of the polynomial.

Let $p(x) = ax^3 + bx^2 + cx + d$ a polynomial of degree $3$ in the polynomial ring of the finite field $GF(2^8)$. We want to find $q(x) = \alpha x^3 + \beta x^2 + \gamma x + \delta$ such that $p(x)q(x) \equiv 1 \bmod x^4 + 1$.

We compute the product $p(x)q(x)$: $$ \begin{array}{rcl} p(x)q(x) & = & a\alpha x^6 + (a\beta + \alpha b) x^5 + (a\gamma + b\beta + c\alpha) x^4 + \\ & & (a\delta + b\gamma + c\beta + d\alpha) x^3 + (b\delta + c\gamma + d\beta) x^2 +\\ & & (c\delta + d\gamma) x + d\delta. \end{array} $$ But we want the product mod $x^4 + 1$, and we have $x^4 \equiv -1 \bmod x^4 + 1$, and even better since we are in a field of characteristic two, we have $x^4 \equiv 1 \bmod x^4 + 1$, so $x^5 \equiv x \bmod x^4 + 1$ and $x^6 \equiv x^2 \bmod x^4 + 1$.

Therefore we have $$ \begin{array}{rcl} p(x)q(x) & \equiv & (a\delta + b\gamma + c\beta + d\alpha) x^3 +\\ & & (b\delta + c\gamma + d\beta + a\alpha) x^2 + \\ & & (c\delta + d\gamma + a\beta + b\alpha) x + \\ & & (d\delta + a\gamma + b\beta + c\alpha) \end{array}\mod x^4 + 1 $$ Since we want $p(x)q(x) \equiv 1 \bmod x^4 + 1$, we have to solve a system of linear equations: $$ \left\{\begin{array}{rcl} a\delta + b\gamma + c\beta + d\alpha & = & 0 \\ b\delta + c\gamma + d\beta + a\alpha & = & 0 \\ c\delta + d\gamma + a\beta + b\alpha & = & 0 \\ d\delta + a\gamma + b\beta + c\alpha & = & 1, \end{array}\right. $$ which can be rewritten as $$ \begin{bmatrix} a & b & c & d \\ b & c & d & a \\ c & d & a & b \\ d & a & b & c \end{bmatrix}\cdot \begin{bmatrix}\delta\\ \gamma \\ \beta \\ \alpha\end{bmatrix} = \begin{bmatrix}0\\0\\0\\1\end{bmatrix} $$ To find the coefficients $\alpha$, $\beta$, $\gamma$ and $\delta$ of the polynomial, we have only to find the inverse of the matrix: $$ \begin{bmatrix}\delta\\ \gamma \\ \beta \\ \alpha\end{bmatrix} = \begin{bmatrix} a & b & c & d \\ b & c & d & a \\ c & d & a & b \\ d & a & b & c \end{bmatrix}^{-1}\cdot\begin{bmatrix}0\\0\\0\\1\end{bmatrix} $$ In fact, the coefficients will be the last column of this matrix.

You can compute the inverse with a method such as Gauss elimination, where all computations are in the field $GF(2^8)$.

In this specific case, the matrix keeping your notation) is: $$ \begin{bmatrix} 03 & 01 & 01 & 02 \\ 01 & 01 & 02 & 03 \\ 01 & 02 & 03 & 01 \\ 02 & 03 & 01 & 01 \end{bmatrix} $$

Whichever method you use, I hope you will get through all the calculation.