In a finite field, there are basically two main algorithms to compute an inverse:
The Extended Euclidean Algorithm, also known as as "extended GCD". There is a variant called Binary GCD. In a finite field of characteristic 2, where addition is bitwise XOR, the binary GCD is easier to implement.
Fermat's little theorem, which is easily generalized: if the field order is $n$ (in your case, $n = 2^8 = 256$), then the inverse of $x$ is $x^{n-2}$, because $x x^{n-2} = x^{n-1} = 1$. This holds because in a field of order $n$, the invertible elements are a multiplicative group of order $n-1$, thus exponentiation by $n-1$ must yield the multiplicative neutral element $1$.
If you already have a multiplication function (mul()), then the exponentiation can be done in relatively few calls to that function. For instance:
static inline unsigned
invPow(unsigned x)
{
unsigned x2 = mul(x, x);
unsigned x3 = mul(x, x2);
unsigned x6 = mul(x3, x3);
unsigned x12 = mul(x6, x6);
unsigned x15 = mul(x12, x3);
unsigned x30 = mul(x15, x15);
unsigned x60 = mul(x30, x30);
unsigned x120 = mul(x60, x60);
unsigned x126 = mul(x120, x6);
unsigned x127 = mul(x126, x);
unsigned x254 = mul(x127, x127);
return x254;
}
In this C function, variables are named out of their contents: x60 contains $x^{60}$. Also note that invPow() returns $0$ for an input of $0$: zero does not have any inverse, but the AES specification formally says that when computing that operation, we use zero as its own "inverse".
Now this will hardly be efficient. You could first precompute a table of inverses, but that would be a look-up table, that you want to avoid.
An optimization is to use a tower of fields: the finite field of size $256$ can be defined as an extension of degree $2$ over the finite field of size $16$. This was described in a succinct note by Rijmen, and allows computing inverses in $\mathbb{F}_{256}$ by doing the relatively simpler task of computing inverses in $\mathbb{F}_{16}$. If you go down that road, you will want to read this article by Boyar and Peralta, who describe the currently best known implementation of the AES S-box in terms of logic gates.
Going down to logic gates is expensive in software, since this means doing more than a hundred boolean operations for each S-box invocation. However, such techniques offer considerable speed-ups if used with bitslicing, especially on machines with large registers. The Käsper-Schwabe implementation is well-known, using 128-bit SSE2 registers. Other similar implementations, using only 32-bit or 64-bit registers, can be found in BearSSL (which also handles the case of the reverse S-box, which is needed for decryption with AES/CBC; the Käsper-Schwabe code supports only the forward S-box, since they target only AES/CTR).
