Last time we talked about basic arithmetic operations in the finite field \( GF(2)[x]/p(x) \) — addition, multiplication, raising to a power, shift-left and shift-right — as well as how to determine whether a polynomial \( p(x) \) is primitive. If a polynomial \( p(x) \) is primitive, it can be used to define an LFSR with coefficients that correspond to the 1 terms in \( p(x) \), that has maximal length of \( 2^N-1 \), covering all bit patterns except the all-zero pattern.

This article and the next three go into detail about multiplicative inverse, the discrete logarithm (which gets two articles!), and an algorithm for finding the characteristic polynomial from LFSR output — and unfortunately I really have mixed feelings about them, in part because we’re drifting away from the nice simple world of embedded systems and off into the faraway land of abstract algebra. (There’s a practical reason for it, though — trust me!) Also because I need to mention some algorithms, and they belong together — they would also fit nicely into my Ten Little Algorithms series, except for the fact that I have a problem I didn’t expect, which is that now I have more than ten little algorithms to present, and that just doesn’t sit well with me. So the Four Little Algorithms are going to be presented separately from the Ten Little Algorithms, and my OCD side will just have to be okay with that.

Anyway, let’s start on that multiplicative inverse.

Multiplicative Inverse in Modular Arithmetic

In Part I I showed this Cayley table for the multiplicative group \( \mathbb{Z}_{13}{}^* \) (multiplication modulo 13):

It’s very easy to multiply modulo 13, you just subtract a multiple of 13 needed to make the result between 0 and 12, so if we have 6*8 we get 48 and divide by 13, yielding 3.something, round down to 3, then subtract 39, and we get 9.

Or just use the % operator in Python, Java, C, and other languages:

(6*8) % 13

9

What if we want to go the other way? What if I want to solve \( 6x = 9 \pmod {13} \)? I could just substitute in numbers for \( x \) until I stumble upon \( x=8 \) and find out that it works. That might be okay for small groups like \( \mathbb{Z} _ {13}{}^ * \). But what about larger multiplicative groups with a prime modulus, like \( \mathbb{Z} _ {268323359541617}{}^ * \), which has over 268 trillion elements? Trial and error isn’t going to get me anywhere.

The key here is that we can first solve for \( 6u = 1 \pmod {13} \) and then compute \( x = 9u \). Furthermore, if we do find that value of \( u \) such that \( 6u = 1 \pmod {13} \), it will mean that \( 6u + 13k = 1 \) for some value of \( k \), so all we need to find are the values \( u \) and \( k \) that work.

There’s this thing called the Extended Euclidean Algorithm which will do just that. The plain Euclidean Algorithm is used to compute greatest common divisor (gcd) and is just a matter of replacing the pair of numbers \( (a,b) \) with \( (b,r) \) where \( r = a \bmod b \), until you hit \( r=0 \), at which point \( b \) is the greatest common divisor of the numbers.

For example:

\(a\)

\(b\)

\(q = \lfloor{a/b}\rfloor\)

\(r = a\bmod b\)

225

147

1

78

147

78

1

69

78

69

1

9

69

9

7

6

9

6

1

3

6

3

2

0

which shows that \( \gcd(225,147) = 3 \).

We can write a fairly simple Python function to do the same thing for us:

That’s fine and dandy, and you don’t even need the quotient \( q \) for the plain gcd itself. The Extended Euclidean algorithm is a way of going backwards and using the \( a, b, q, r \) values to calculate Bézout’s identity, that is, values \( x \) and \( y \) such that \( ax+by = \gcd(a,b) \), but then you have to maintain a history of \( a, b, q, r \) values in an array or a stack or a list, and I can never actually remember how to do it. Bleah.

In 1963, W. A. Blankinship described a variant of this algorithm, that doesn’t require you to keep around state except for a fixed set of variables, and when you get to the GCD you also get the \( x \) and \( y \) values. It’s the iterative alternative to the original Extended Euclidean algorithm, which is recursive.

Blankinship’s Algorithm is really simple and it works like this. You have a 2x3 matrix containing an \( a \) row and a \( b \) row, and it looks like this:

\( \begin{bmatrix}a & 1 & 0 \cr b & 0 & 1\end{bmatrix} \)

Then you do the same steps just as before with \( a \) and \( b \), but instead of just computing \( r = a \bmod b \), you write it as \( r = a - bq \). Once you have \( r \) and \( q \), you take the entire \( a \) row and subtract \( q \) times the \( b \) row, then throw away the \( a \) row, promote the \( b \) row to the new \( a \) row, and then the \( r \) row becomes the \( b \) row.

Here’s the algorithm in action:

step #\( k\)

\(a_k\)

\(x_k\)

\(y_k\)

\(q_k = \lfloor{a_{k-1}/a_k}\rfloor\)

0

225

1

0

1

147

0

1

1

2

78

1

-1

1

3

69

-1

2

1

4

9

2

-3

7

5

6

-15

23

1

6

3

17

-26

2

7

0

And this tells you that \( 17\times225 - 26\times 147 = 3 \). It works because of a very simple invariant; at each step, \( a_k = ax_k + by_k \). For example:

Big deal. We didn’t need a special algorithm; we could have figured that out on our own,
by hand, but that’s only because the numbers are small.
With Blankinship’s algorihm, you can do math in \( \mathbb{Z}_{268323359541617}{}^* \) quite easily:

Just in case you didn’t follow that, the quotient calculation \( q_k = \lfloor{a_{k-1}/a_k}\rfloor \) is a little handwave-y; at each step we take \( a_{k-1}(x) \) and \( a_k(x) \) and do long polynomial division, using the quotient necessary to reduce the degree as much as possible. For example, in the first case, \( (x^4 + x+1)/(x^3+x) \) we want to get rid of the \( x^4 \) in the dividend, so we try subtracting off \( x(x^3+x) \):

Honestly, for the inverse in the finite field \( GF(2^n) \) (= quotient ring), I’m not really sure; both it and division are not really something I’ve needed to use. But it’s there in libgf2 in case you need it.

On the other hand, the “regular” version of Blankinship’s algorithm, the one that operates on plain integers rather than polynomials, is something used very commonly in discrete mathematics where you need to “undo” raising to a power (= calculating a \( k \)th root) in a group or a finite field, because even if you have the most bizarre representation of a group or field in the world, with matrices or elliptic curves or Rubik’s Cubes or whatever, raising to a power still involves plain old integer exponents. For example, suppose there is a generator element \( g \) such that \( g^m = 1 \) where \( m \) is the period of the generator. If I have an element \( x = g^k \) and I know \( k \) and \( x \) and \( m \), but I don’t know \( g \), I can calculate \( x^j = (g^k)^j \) and if \( jk = 1 \pmod m \) then \( x^j = (g^k)^j = g^{jk} = g^{bm+1} = g\cdot(g^m)^b = g \), so the required \( j \) is the multiplicative inverse of \( k \pmod m \), and all I have to do is call a,j,v = blankinship(k,m), verify that \( a=1 \), and that means \( jk = 1 \pmod m \). Then I just calculate \( g = x^j \). This is part of what makes the RSA cryptosystem possible: the encryption exponent \( e \) and the decryption exponent \( d \) are multiplicative inverses modulo \( m \), but it is essentially impossible to know \( m \) — unless you know how to factor the RSA modulus \( pq \) into its two prime factors \( p \) and \( q \), and then it’s easy to calculate \( m = \phi(pq) = (p-1)(q-1) \).

So whether you realize it or not, chances are you use multiplicative inverses every day, whenever you point a web browser at an https:// URL, because the security certificates use RSA for digital signatures to assert the web site’s identity.

Wrapup

What’s to say here? This was a relatively short article, covering the Extended Euclidean Algorithm and Blankinship’s iterative version of it that uses constant space. We looked at it both for integers and elements of a finite field, and showed how to use it to calculate the multiplicative inverse.