How to compute the inverse of a polynomial under $GF(2^8)$? It's well known that AES cryptography algorithm uses Galois Field $GF(2^8)$ multiplication to process the step MixColumn, and each column of the 4x4 matrix on encrypting should multiply the polynomial $3X^3 + X^2 + X + 2$ which is usually notated as an array $\{\tt 0x03, 0x01, 0x01, 0x02\}$, while on decrypting, the matrix should multiply the inverse of the polynominal mentioned above, namely $\texttt{0x0B}X^3 + \texttt{0x0D}X^2 + \texttt{0x09}X + \texttt{0x0E}$ which can also notated as an array $\{\tt 0x0B, 0x0D, 0x09, 0x0E\}$.

Now I wonder If only one polynomial is known, e.g. $\{\tt 0x03, 0x01, 0x01, 0x02\}$, how can I compute its inverse mod $(X^4 + 1)$, namely $\{\tt 0x0B, 0x0D, 0x09, 0x0E\}$, and vice versa.
The classic generic algorithm for computing modular inverses is the Extended Euclidean Algorithm. The algorithm is primarily defined for integers, but in fact it works for all rings where you can define a notion of Euclidean division (i.e. "Euclidean domains"). In particular it works with polynomials whose coefficients are in any field.

Let $K$ be a field. In the case of AES, this is $\mathbb{F}_2^8$, the finite field with 256 elements, but this does not actually matter here. $K[X]$ is then the ring of polynomials with coefficients in $K$. Given two such polynomials $A$ and $B$, you can find polynomials $Q$ and $R$ such that $A = BQ + R$, and the degree of $R$ (the largest power of $X$ for which the corresponding coefficient of $R$ is non-zero) is strictly less than the degree of $B$. This is easily computed in about the same way you would compute a division of integers; in fact it is easier because there is no carry to propagate between coefficients.

The magic part of it is that the $Q$ and $R$ polynomials are unique; i.e. for two polynomials $A$ and $B$, there is only one pair of polynomials $(Q,R)$ that matches the conditions above. Thus, the Euclidean division of polynomials is well defined, and the Extended Euclidean algorithm works.

In practice, the Extended Euclidean Algorithm is about finding the Greatest Common Divisor. In the context of polynomials, "greatest" means "of largest degree". Note that this GCD is defined only up to a multiplicative constant; it is often taken as convention that the GCD should have its largest non-zero coefficient equal to 1 (you can always arrange for that, because if polynomial $B$ divides polynomial $A$, then $tB$ also divides $A$ for all non-zero field elements $t$).



So if you have a polynomial $A$ and want to inverse it modulo $B = X^4+1$, then you run the Extended Euclidean Algorithm to obtain polynomials $U$, $V$ and $G$ such that $G$ is the GCD of $A$ and $B$, and $AU+BV = G$.

If $G$ has degree 0 (i.e. it is the constant $1$), then $A$ is invertible modulo $B$, and its inverse is $U$.

If $G$ has a non-zero degree, then this means that $A$ and $B$ can both be divided by that non-trivial polynomial, and there is no solution: $A$ is then not invertible modulo $B$.

Note that $X^4+1$ is not irreducible over $\mathbb{F}_2^8$; it is equal to $(X+1)^4$. It follows that a polynomial $P$ of degree at most 3 in $\mathbb{F}_2^8$ is invertible modulo $X^4+1$ if and only $P$ is not divisible by $X+1$, i.e. $P(1) \neq 0$. For $P = 3·X^3+X^2+X+2$ (I am calling "$3$" what the AES specification denotes 0x03), you have $P(1) = 3 + 1 + 1 + 2 = 1$ (remember, we are working in $\mathbb{F}_2^8$, so addition is XOR), so that polynomial $P$ is indeed invertible.



Another, completely different way of computing inverses is with modular exponentiation. The polynomials which are invertible modulo some given polynomial $B$ form a multiplicative group. If that group has order $e$, then the inverse of any invertible polynomial $A$ modulo $B$ is necessarily $A^{e-1} \pmod B$. In the case of $B = X^4+1$, the invertible polynomial modulo $B$ are all the polynomials $P$ such that $P(1) \neq 0$, so $e = 255·256^3 = 4278190080$ (this is easily seen in the following way: you can choose three coefficient arbitrarily, but for the fourth coefficient, you must avoid the one value that would make $P(1) = 0$). Therefore, to compute the inverse of a given polynomial $A$ modulo $X^4+1$, you "just" have to compute $A^{4278190079} \pmod {X^4+1}$, which can be done with a square-and-multiple algorithm.

That method is less efficient than the Extended Euclidean Algorithm, but possibly simpler to implement.