Linear Feedback Shift Registers for the Uninitiated, Part III: Multiplicative Inverse, and Blankinship's Algorithm

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):

from IPython.core.display import HTML

def binop_table(operator,operator_symbol,operand_set):
    def row(values,celltype='td'):
        return ('<'+celltype+'>'
                +('</%s><%s>' % (celltype,celltype)).join('%d' % v for v in values)
                +'</'+celltype+'>')
    return HTML('<table class="align-right"><tr><th>'+operator_symbol+'</th>'
                    +row(operand_set,celltype='th')
                    +'</tr><tr>'
                +'</tr>\n<tr>'.join('<th>%d</th>'%k
                                    +row(operator(j,k) for j in operand_set) 
                                    for k in operand_set)
                +'</tr></table>')    

def mulmod13(a,b):
    return (a*b)%13
binop_table(mulmod13,'&times;',xrange(1,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:

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

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

def gcdverbose(a,b):
    while b != 0:
        q,r = divmod(a, b)
        print a,b,q,r
        a,b = b,r
    return a

gcdverbose(225,147)

225 147 1 78
147 78 1 69
78 69 1 9
69 9 7 6
9 6 1 3
6 3 2 0

3

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:

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:

$$\begin{align} 225 &= 1\times 225 + 0\times 147, \cr 147 &=0\times225 + 1\times 147, \cr 78 &=1\times225 + -1\times 147, \cr 69 &=-1\times225 + 2\times 147, \cr \end{align} $$

and so on. You can subtract any multiple of any of these equalities from each other and you’re left with a new equality, and the last one with a nonzero \( a_k \) term is the final result.

This is easy to write in most computer languages:

def blankinship(a,b, verbose=False):
    arow = [a,1,0]
    brow = [b,0,1]
    if verbose:
        print arow
    while brow[0] != 0:
        q,r = divmod(arow[0], brow[0])
        if verbose:
            print brow, q
        if r == 0:
            break
        rrow =[r, arow[1]-q*brow[1], arow[2]-q*brow[2]]
        arow = brow
        brow = rrow
    return brow[0], brow[1], brow[2]
        
g,x,y = blankinship(225,147,True)
print g, x*225 + y*147

[225, 1, 0]
[147, 0, 1] 1
[78, 1, -1] 1
[69, -1, 2] 1
[9, 2, -3] 7
[6, -15, 23] 1
[3, 17, -26] 2
3 3

Why do we care? Because it tells us the multiplicative inverse:

blankinship(6,13)

(1, -2, 1)

This tells us that \( \gcd(6,13) = 1 = 6\times (-2) + 13\times 1 \). And that means the multiplicative inverse of \( 6 \), mod 13, is \( -2 \) which can be canonicalized to \( 13-2 = 11 \):

(6*11) % 13

1

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:

m = 268323359541617
a = 1234567
g,x,y=blankinship(a, m, verbose=True)
print "-----"
print "%d * %d = %d (mod %d)" % (x,a,(a*x) % m,m)

[1234567, 1, 0]
[268323359541617, 0, 1] 0
[1234567, 1, 0] 217342079
[1096824, -217342079, 1] 1
[137743, 217342080, -1] 7
[132623, -1738736639, 8] 1
[5120, 1956078719, -9] 25
[4623, -50640704614, 233] 1
[497, 52596783333, -242] 9
[150, -524011754611, 2411] 3
[47, 1624632047166, -7475] 3
[9, -5397907896109, 24836] 5
[2, 28614171527711, -131655] 4
[1, -119854594006953, 551456] 2
-----
-119854594006953 * 1234567 = 1 (mod 268323359541617)

The same idea works in finite fields, where elements have a multiplicative inverse that can also be calculated by Blankinship’s algorithm; we just need to know what kind of arithmetic to use here.

• For elements in multiplicative groups \( \mathbb{Z} _ {p}{}^ * \), where, given some element \( a \), we want to compute \( x \) such that \( ax = 1 \pmod p \), we run Blankinship’s algorithm in normal integer arithmetic (\( \mathbb{Z} \)) — we drop the \( \bmod p \) step and just solve \( ax + py = 1 \).
• For elements in finite fields \( GF(2)[x]/p(x) \), where given some element \( a(x) \), we want to compute \( u(x) \) such that \( a(x)u(x) = 1 \pmod {p(x)} \), we run Blankinship’s algorithm in polynomial ring arithmetic (\( GF(2)[x] \)) — we also drop the \( \bmod p(x) \) step and just solve \( a(x)u(x) + p(x)v(x) = 1 \).

Let’s take a really simple example with \( p(x) = x^4 + x + 1 \) and \( a(x) = x^3 + x \) and run Blankinship’s algorithm:

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) \):

$$ \begin{alignat}{5} x^4 &\phantom{+x^3} & &+x &+ 1 \cr x^4 & &+ x^2 & & \cr \hline \cr & & x^2&+x&+1 \end{alignat} $$

The remainder \( x^2 + x+1 \) is a smaller degree than the divisor \( x^3+x \) so we’re done.

Here’s another example with the following step, \( (x^3+x)/(x^2+x+1) \) where we subtract \( x(x^2+x+1) \) and then \( 1(x^2+x+1) \) so the quotient is \( x+1 \) and the remainder is also \( x+1 \):

$$ \begin{alignat}{3} x^3 & &+x & \cr x^3 &+ x^2 &+x & \cr \hline \cr &\phantom{+} x^2 & & \cr &\phantom{+} x^2 &+x &+1 \cr \hline \cr & & x &+1 \end{alignat} $$

Anyway, the table above tells us that \( 1 = (x^2+x+1)(x^4+x+1) + (x^3+x^2)(x^3+x) \) with arithmetic in \( GF(2)[x] \), so the multiplicative inverse of \( x^3 + x \pmod {p(x)} \) is \( x^3 + x^2 \).

We can verify this with libgf2, the Python module I wrote for \( GF(2) \) computation, which no one else actually uses:

qr = libgf2.gf2.GF2QuotientRing(0b10011)  # x^4 + x + 1
a = qr.wrap(0b1010)                       # x^3 + x
v = a.inv
v

GF2Element(0b1100,0x13)

a*v

GF2Element(0b0001,0x13)

What’s the use of a multiplicative inverse in discrete mathematics?

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.

Next time we’ll get into some good stuff: the discrete logarithm.

© 2017 Jason M. Sachs, all rights reserved.

Sign in