Linear Feedback Shift Registers for the Uninitiated, Part IX: Decimation, Trace Parity, and Cyclotomic Cosets

Jason SachsDecember 3, 2017

Last time we looked at matrix methods and how they can be used to analyze two important aspects of LFSRs:

  • time shifts
  • state recovery from LFSR output

In both cases we were able to use a finite field or bitwise approach to arrive at the same result as a matrix-based approach. The matrix approach is more expensive in terms of execution time and memory storage, but in some cases is conceptually simpler.

This article will be covering some concepts that are useful for studying the decimation of maximal-length LFSR output.

Let’s look again at our old friend, the LFSR that can be represented by the quotient ring \( H_{25} = GF(2)[x]/(x^5 + x^2 + 1) \):

Its output is 0000100101100111110001101110101, repeated over and over again. But you don’t have to take my word for it; we can compute this in Python. While we do this, let’s try seeing what happens if we take every \( j \)th output bit. In other words, if the raw output sequence of \( H_{25} \) is \( y[k] = \{0,0,0,0,1,0,0,\ldots\} \), so that we know how to generate the sequence \( y[0], y[1], y[2], y[3], \ldots \), then how can we determine the sequence \( y[0], y[2], y[4], y[6], \ldots \) or the sequence \( y[0], y[3], y[6], y[9], \ldots \) or in general the sequence \( y[0], y[j], y[2j], y[3j], \ldots \)?

from libgf2.gf2 import GF2QuotientRing

H25 = GF2QuotientRing(0x25)
e = H25.wrap(1)
def highbit(u):
    return (u.coeffs >> (u.field.degree-1)) & 1

for j in xrange(1,32):
    print 'j=%2d' % j, ''.join('01'[highbit(e << (j*k))] for k in xrange(36)) + '...'
j= 1 000010010110011111000110111010100001...
j= 2 001001011001111100011011101010000100...
j= 3 000101011010000110010011111011100010...
j= 4 010010110011111000110111010100001001...
j= 5 001101111101000100101011000011100110...
j= 6 000011001001111101110001010110100001...
j= 7 011111001001100001011010100011101111...
j= 8 001101110101000010010110011111000110...
j= 9 010001001010110000111001101111101000...
j=10 010110000111001101111101000100101011...
j=11 001011111011001110000110101001000101...
j=12 001010110100001100100111110111000101...
j=13 011010100100010111110110011100001101...
j=14 011010100011101111100100110000101101...
j=15 011101100011111001101001000010101110...
j=16 010100001001011001111100011011101010...
j=17 010000110010011111011100010101101000...
j=18 000011100110111110100010010101100001...
j=19 001110111110010011000010110101000111...
j=20 001001010110000111001101111101000100...
j=21 010010001011111011001110000110101001...
j=22 011111011001110000110101001000101111...
j=23 001111100110100100001010111011000111...
j=24 011100010101101000011001001111101110...
j=25 010110101000111011111001001100001011...
j=26 011100001101010010001011111011001110...
j=27 000010101110110001111100110100100001...
j=28 011101111100100110000101101010001110...
j=29 000101011101100011111001101001000010...
j=30 010101110110001111100110100100001010...
j=31 000000000000000000000000000000000000...

Something interesting just happened. Did you notice? Let’s look at the cases \( j=1, 2, 4 \):

j= 1 000010010110011111000110111010100001...
j= 2 001001011001111100011011101010000100...
j= 4 010010110011111000110111010100001001...

It’s the same repeating sequence, just shifted a little bit! We’ll be able to explain why this is the case a little bit later, but for now let’s take the original sequence with \( j=1 \) and show it twice, so this curious relationship is easier to notice:

  • once with every other element highlighted (the highlighted elements are the sequence for \( j=2 \))
  • then with every element shown but spread outwards
from IPython.core.display import HTML

def highlight_every_other_bit(jlist, prefixes, bitlist):
    html = """
<div><style type="text/css" scoped="scoped">
    code span.hl { font-weight: bold; }
</style><code>"""
    for j,prefix,bits in zip(jlist,prefixes,bitlist):
        html += prefix
        for k, bit in enumerate(bits):
            c = '01'[bit]
            if j == 2:
                html += c + ' '
            elif k & 1:
                html += c
            else:
                html += "<span class='hl'>" + c + "</span>"
        html += '</br>'
    html += '</code></div>'
    return HTML(html)

highlight_every_other_bit([1,2],
                          ['&nbsp;'*4,''],
                          [(highbit(e << k) for k in xrange(62)),
                           (highbit(e << k) for k in xrange(33))])
    00001001011001111100011011101010000100101100111110001101110101
0 0 0 0 1 0 0 1 0 1 1 0 0 1 1 1 1 1 0 0 0 1 1 0 1 1 1 0 1 0 1 0 0

Huh.

Since it’s rather difficult to tell at a glance that two repeating sequences are the same, let’s use the Berlekamp-Massey algorithm to identify the minimal polynomial that generates the sequence in question, for several different values of decimation ratio \( j \):

from libgf2.util import berlekamp_massey
from libgf2.gf2 import checkPeriod

polymap = {}
for j in xrange(1,32):
    poly, deg = berlekamp_massey([highbit(e << (j*k)) for k in xrange(10)])
    if poly not in polymap:
        polymap[poly] = []
    polymap[poly].append(j)
    
for k,v in polymap.iteritems():
    print 'poly=%02x, j=%s' % (k,v)
poly=01, j=[31]
poly=25, j=[1, 2, 4, 8, 16]
poly=29, j=[15, 23, 27, 29, 30]
poly=2f, j=[7, 14, 19, 25, 28]
poly=37, j=[5, 9, 10, 18, 20]
poly=3b, j=[11, 13, 21, 22, 26]
poly=3d, j=[3, 6, 12, 17, 24]

Notice that decimation by any of the powers of two yields a sequence that has the same minimal polynomial, specifically that of the original quotient ring \( H_{25} \).

It turns out that decimation ratios are partitioned into something called cyclotomic cosets, which are the sets of integers \( {j\cdot q^k \bmod q^N-1} \), where \( N \) is the polynomial degree and \( q \) is the field characteristic — that’s just \( q=2 \) for \( GF(2^N) \). For example, with \( N=5 \) we are dealing with \( \bmod 31 \) and the cyclotomic cosets are

  • \( \{0 \equiv 31\} \)
  • \( \{1,2,4,8,16\} \)
  • \( \{3, 6, 12, 24, 48\} \equiv \{3,6,12,24,17\} \) (since \( 48 \bmod 31 = 17 \))
  • \( \{5, 10, 20, 40, 80\} \equiv \{5,10,20,9,18\} \)
  • \( \{7, 14, 28, 56, 112 \} \equiv \{7,14,28,25,19\} \)
  • \( \{11, 22, 44, 88, 176\} \equiv \{11, 22, 13, 26, 21 \} \)
  • \( \{15, 30, 60, 120, 240\} \equiv \{15, 30, 29, 27, 23 \} \)

The smallest members of each coset — \( \{0,1,3,5,7,11,15\} \) — are called the coset representatives.

There’s something special as well about the polynomials we determined from feeding decimation output into Berlekamp-Massey; if the decimation ratio \( j \) is relatively prime to the period \( 2^N-1 \), then each polynomial determined by Berlekamp-Massey is the minimal polynomial of the elements \( x^j \) in each cyclotomic coset, and the decimated sequence has the same period as the original sequence. By minimal polynomial we mean the polynomial of smallest degree that yields zero when evaluated with the specified element; a primitive element is one that generates all the elements in the field — in our example of \( GF(2^5) \) an element \( x^j \) is primitive if all the \( x^{jk} \) are unique up to \( k=31 \).

For example, with \( j=1 \), we have \( p(x) = x^5 + x^2 + 1 = 0 \) because this is the characteristic polynomial of the quotient ring \( H_{25} \). But \( p(x^2) = x^{10} + x^4 + 1 \) is also equal to zero, because it is equal to \( p(x)^2 = 0^2 = 0 \). (When computing any \( (a+b)^2 = a^2 + 2ab + b^2 \) in fields of characteristic 2, the \( 2ab \) term is zero modulo 2 and drops out, leaving only the individual squared terms: \( (a+b)^2 \equiv a^2 + b^2 \). This is called the “freshman’s dream” identity.)

If the decimation ratio \( j \) is not relatively prime to the period \( 2^N-1 \), then we still get the same minimal polynomial from Berlekamp-Massey as the minimum polynomial of \( x^j \), but the decimated sequence will have a shorter period than the original sequence, and the minimal polynomial may have a degree that is less than \( N \).

e25 = H25.wrap(1)
print "p(x) = x^5 + x^2 + 1"
for j in [1,2,4,8,16]:
    xj = e25 << j
    print "j=%2d -> x^j = %s, p(x^j) = %s" % (j, xj, xj**5 + xj**2 + 1)
p(x) = x^5 + x^2 + 1
j= 1 -> x^j = GF2Element(0b00010,0x25), p(x^j) = GF2Element(0b00000,0x25)
j= 2 -> x^j = GF2Element(0b00100,0x25), p(x^j) = GF2Element(0b00000,0x25)
j= 4 -> x^j = GF2Element(0b10000,0x25), p(x^j) = GF2Element(0b00000,0x25)
j= 8 -> x^j = GF2Element(0b01101,0x25), p(x^j) = GF2Element(0b00000,0x25)
j=16 -> x^j = GF2Element(0b11011,0x25), p(x^j) = GF2Element(0b00000,0x25)

Because the period of \( H_{25} \) is a prime (31), all the decimation ratios yield primitive polynomials, and we get the same kind of calculations for all of the other cosets:

print "p(x) = x^5 + x^4 + x^3 + x^2 + 1"
for j in [3,6,12,24,17]:
    xj = e25 << j
    print "j=%2d -> x^j = %s, p(x^j) = %s" % (j, xj, xj**5 + xj**4 + xj**3 + xj**2 + 1)

print "\np(x) = x^5 + x^4 + x^2 + x + 1"
for j in [5,10,20,9,18]:
    xj = e25 << j
    print "j=%2d -> x^j = %s, p(x^j) = %s" % (j, xj, xj**5 + xj**4 + xj**2 + xj + 1)    
    
print "\np(x) = x^5 + x^3 + x^2 + x + 1"
for j in [7,14,28,25,19]:
    xj = e25 << j
    print "j=%2d -> x^j = %s, p(x^j) = %s" % (j, xj, xj**5 + xj**3 + xj**2 + xj + 1)    
    
print "\np(x) = x^5 + x^4 + x^3 + x + 1"
for j in [11, 22, 13, 26, 21]:
    xj = e25 << j
    print "j=%2d -> x^j = %s, p(x^j) = %s" % (j, xj, xj**5 + xj**4 + xj**3 + xj + 1)
        
print "\np(x) = x^5 + x^3 + 1"
for j in [15, 30, 29, 27, 23]:
    xj = e25 << j
    print "j=%2d -> x^j = %s, p(x^j) = %s" % (j, xj, xj**5 + xj**3 + 1)    
    
p(x) = x^5 + x^4 + x^3 + x^2 + 1
j= 3 -> x^j = GF2Element(0b01000,0x25), p(x^j) = GF2Element(0b00000,0x25)
j= 6 -> x^j = GF2Element(0b01010,0x25), p(x^j) = GF2Element(0b00000,0x25)
j=12 -> x^j = GF2Element(0b01110,0x25), p(x^j) = GF2Element(0b00000,0x25)
j=24 -> x^j = GF2Element(0b11110,0x25), p(x^j) = GF2Element(0b00000,0x25)
j=17 -> x^j = GF2Element(0b10011,0x25), p(x^j) = GF2Element(0b00000,0x25)

p(x) = x^5 + x^4 + x^2 + x + 1
j= 5 -> x^j = GF2Element(0b00101,0x25), p(x^j) = GF2Element(0b00000,0x25)
j=10 -> x^j = GF2Element(0b10001,0x25), p(x^j) = GF2Element(0b00000,0x25)
j=20 -> x^j = GF2Element(0b01100,0x25), p(x^j) = GF2Element(0b00000,0x25)
j= 9 -> x^j = GF2Element(0b11010,0x25), p(x^j) = GF2Element(0b00000,0x25)
j=18 -> x^j = GF2Element(0b00011,0x25), p(x^j) = GF2Element(0b00000,0x25)

p(x) = x^5 + x^3 + x^2 + x + 1
j= 7 -> x^j = GF2Element(0b10100,0x25), p(x^j) = GF2Element(0b00000,0x25)
j=14 -> x^j = GF2Element(0b11101,0x25), p(x^j) = GF2Element(0b00000,0x25)
j=28 -> x^j = GF2Element(0b10110,0x25), p(x^j) = GF2Element(0b00000,0x25)
j=25 -> x^j = GF2Element(0b11001,0x25), p(x^j) = GF2Element(0b00000,0x25)
j=19 -> x^j = GF2Element(0b00110,0x25), p(x^j) = GF2Element(0b00000,0x25)

p(x) = x^5 + x^4 + x^3 + x + 1
j=11 -> x^j = GF2Element(0b00111,0x25), p(x^j) = GF2Element(0b00000,0x25)
j=22 -> x^j = GF2Element(0b10101,0x25), p(x^j) = GF2Element(0b00000,0x25)
j=13 -> x^j = GF2Element(0b11100,0x25), p(x^j) = GF2Element(0b00000,0x25)
j=26 -> x^j = GF2Element(0b10111,0x25), p(x^j) = GF2Element(0b00000,0x25)
j=21 -> x^j = GF2Element(0b11000,0x25), p(x^j) = GF2Element(0b00000,0x25)

p(x) = x^5 + x^3 + 1
j=15 -> x^j = GF2Element(0b11111,0x25), p(x^j) = GF2Element(0b00000,0x25)
j=30 -> x^j = GF2Element(0b10010,0x25), p(x^j) = GF2Element(0b00000,0x25)
j=29 -> x^j = GF2Element(0b01001,0x25), p(x^j) = GF2Element(0b00000,0x25)
j=27 -> x^j = GF2Element(0b01011,0x25), p(x^j) = GF2Element(0b00000,0x25)
j=23 -> x^j = GF2Element(0b01111,0x25), p(x^j) = GF2Element(0b00000,0x25)

All the elements are primitive for \( N=5 \), so these are nice and symmetrical, since \( 2^N-1=31 \) is a prime; for composite order, things get a little messier:

H43 = GF2QuotientRing(0x43)
e43 = H43.wrap(1)

polymap = {}
for j in xrange(1,64):
    poly, deg = berlekamp_massey([highbit(e43 << (j*k)) for k in xrange(12)])
    if poly not in polymap:
        polymap[poly] = []
    polymap[poly].append(j)
    
for k,v in polymap.iteritems():
    print 'poly=%02x, j=%s' % (k,v)
    
print "\np(x) = x^2 + x + 1"
for j in polymap[7]:
    xj = e43 << j
    print "j=%2d -> x^j = %s, p(x^j) = %s" % (j, xj, xj**2 + xj + 1)  
    
print "\np(x) = x^6 + x^4 + x^2 + x + 1"
for j in polymap[0x57]:
    xj = e43 << j
    print "j=%2d -> x^j = %s, p(x^j) = %s" % (j, xj, xj**6 + xj**4 + xj**2 + xj + 1)  
    
print "\np(x) = x^6 + x^3 + 1"
for j in polymap[0x49]:
    xj = e43 << j
    print "j=%2d -> x^j = %s, p(x^j) = %s" % (j, xj, xj**6 + xj**3 + 1)  

print "\np(x) = x^6 + x^5 + 1"
for j in polymap[0x61]:
    xj = e43 << j
    print "j=%2d -> x^j = %s, p(x^j) = %s" % (j, xj, xj**6 + xj**5 + 1)  

print "\np(x) = x^2 + x + 1"
for j in polymap[0x7]:
    xj = e43 << j
    print "j=%2d -> x^j = %s, p(x^j) = %s" % (j, xj, xj**2 + xj + 1)  
poly=01, j=[9, 18, 27, 36, 45, 54, 63]
poly=43, j=[1, 2, 4, 8, 16, 32]
poly=61, j=[31, 47, 55, 59, 61, 62]
poly=67, j=[5, 10, 17, 20, 34, 40]
poly=49, j=[7, 14, 28, 35, 49, 56]
poly=07, j=[21, 42]
poly=6d, j=[11, 22, 25, 37, 44, 50]
poly=73, j=[23, 29, 43, 46, 53, 58]
poly=75, j=[15, 30, 39, 51, 57, 60]
poly=57, j=[3, 6, 12, 24, 33, 48]
poly=5b, j=[13, 19, 26, 38, 41, 52]

p(x) = x^2 + x + 1
j=21 -> x^j = GF2Element(0b111011,0x43), p(x^j) = GF2Element(0b000000,0x43)
j=42 -> x^j = GF2Element(0b111010,0x43), p(x^j) = GF2Element(0b000000,0x43)

p(x) = x^6 + x^4 + x^2 + x + 1
j= 3 -> x^j = GF2Element(0b001000,0x43), p(x^j) = GF2Element(0b000000,0x43)
j= 6 -> x^j = GF2Element(0b000011,0x43), p(x^j) = GF2Element(0b000000,0x43)
j=12 -> x^j = GF2Element(0b000101,0x43), p(x^j) = GF2Element(0b000000,0x43)
j=24 -> x^j = GF2Element(0b010001,0x43), p(x^j) = GF2Element(0b000000,0x43)
j=33 -> x^j = GF2Element(0b010010,0x43), p(x^j) = GF2Element(0b000000,0x43)
j=48 -> x^j = GF2Element(0b001101,0x43), p(x^j) = GF2Element(0b000000,0x43)

p(x) = x^6 + x^3 + 1
j= 7 -> x^j = GF2Element(0b000110,0x43), p(x^j) = GF2Element(0b000000,0x43)
j=14 -> x^j = GF2Element(0b010100,0x43), p(x^j) = GF2Element(0b000000,0x43)
j=28 -> x^j = GF2Element(0b011100,0x43), p(x^j) = GF2Element(0b000000,0x43)
j=35 -> x^j = GF2Element(0b001011,0x43), p(x^j) = GF2Element(0b000000,0x43)
j=49 -> x^j = GF2Element(0b011010,0x43), p(x^j) = GF2Element(0b000000,0x43)
j=56 -> x^j = GF2Element(0b011111,0x43), p(x^j) = GF2Element(0b000000,0x43)

p(x) = x^6 + x^5 + 1
j=31 -> x^j = GF2Element(0b100101,0x43), p(x^j) = GF2Element(0b000000,0x43)
j=47 -> x^j = GF2Element(0b100111,0x43), p(x^j) = GF2Element(0b000000,0x43)
j=55 -> x^j = GF2Element(0b101110,0x43), p(x^j) = GF2Element(0b000000,0x43)
j=59 -> x^j = GF2Element(0b111101,0x43), p(x^j) = GF2Element(0b000000,0x43)
j=61 -> x^j = GF2Element(0b110001,0x43), p(x^j) = GF2Element(0b000000,0x43)
j=62 -> x^j = GF2Element(0b100001,0x43), p(x^j) = GF2Element(0b000000,0x43)

p(x) = x^2 + x + 1
j=21 -> x^j = GF2Element(0b111011,0x43), p(x^j) = GF2Element(0b000000,0x43)
j=42 -> x^j = GF2Element(0b111010,0x43), p(x^j) = GF2Element(0b000000,0x43)

OK, blah blah blah minimal polynomial primitive element, blah blah blah. Big deal. Why does this matter?

Well, suppose we want to take every 47th output bit of the LFSR defined by \( H_{43} \) with polynomial 0x43 = 0b1000011 \( \Longleftrightarrow p_{43}(x) = x^6 + x + 1 \). The math tells us that the minimal polynomial of \( x^{47} \) is \( p_{61}(x) = x^6 + x^5 + 1 \Longleftrightarrow \) 0b1100001 = 0x61, so that means we can use an LFSR \( H_{61} \) with this polynomial and some particular initial state to yield the same output sequence as every 47th sample from the output of the LFSR \( H_{43} \):

from libgf2.util import state_from_reversed_output
from libgf2.gf2 import GF2Element

def decimate_LFSR(element, j):
    n = element.field.degree
    decimated_sequence = [highbit(element << (j*k)) for k in xrange(2*n)]
    # First figure out the minimal polynomial of the decimated sequence 
    poly, _ = berlekamp_massey(decimated_sequence)
    # Next figure out the LFSR state that corresponds to the decimated sequence
    state = state_from_reversed_output(poly, decimated_sequence[::-1])
    return GF2Element(state, poly)

j=47
e2 = decimate_LFSR(e, j)

print "       output from %s decimated with j=%d\nyields output from %s" %(e,j,e2)

for k in xrange(20):
    print highbit(e<<(j*k)), highbit(e2<<k), e<<(j*k), e2<<(k)
       output from GF2Element(0b00001,0x25) decimated with j=47
yields output from GF2Element(0b01011,0x25)
0 0 GF2Element(0b00001,0x25) GF2Element(0b01011,0x25)
1 1 GF2Element(0b11011,0x25) GF2Element(0b10110,0x25)
0 0 GF2Element(0b00010,0x25) GF2Element(0b01001,0x25)
1 1 GF2Element(0b10011,0x25) GF2Element(0b10010,0x25)
0 0 GF2Element(0b00100,0x25) GF2Element(0b00001,0x25)
0 0 GF2Element(0b00011,0x25) GF2Element(0b00010,0x25)
0 0 GF2Element(0b01000,0x25) GF2Element(0b00100,0x25)
0 0 GF2Element(0b00110,0x25) GF2Element(0b01000,0x25)
1 1 GF2Element(0b10000,0x25) GF2Element(0b10000,0x25)
0 0 GF2Element(0b01100,0x25) GF2Element(0b00101,0x25)
0 0 GF2Element(0b00101,0x25) GF2Element(0b01010,0x25)
1 1 GF2Element(0b11000,0x25) GF2Element(0b10100,0x25)
0 0 GF2Element(0b01010,0x25) GF2Element(0b01101,0x25)
1 1 GF2Element(0b10101,0x25) GF2Element(0b11010,0x25)
1 1 GF2Element(0b10100,0x25) GF2Element(0b10001,0x25)
0 0 GF2Element(0b01111,0x25) GF2Element(0b00111,0x25)
0 0 GF2Element(0b01101,0x25) GF2Element(0b01110,0x25)
1 1 GF2Element(0b11110,0x25) GF2Element(0b11100,0x25)
1 1 GF2Element(0b11010,0x25) GF2Element(0b11101,0x25)
1 1 GF2Element(0b11001,0x25) GF2Element(0b11111,0x25)

Here we have two LFSRs with different characteristic polynomials, and sampling one every 47 samples yields the same output as sampling the other normally (every sample). Note that the LFSR state is different, but the output — remember, this is just the high bit of the LFSR state — is the same.

Trace Parity

Okay, now things get weird.

I mentioned trace maps very briefly in the last article as some hard-to-understand concept. Since then I’ve acquired a copy of Lidl & Niederreiter’s Introduction to Finite Fields and Their Applications and stared at Chapter 2 until I’ve got Greek letters coming out of my ears, and some of it has made sense. So I’ve reached ex-Pralite monk status, just barely.

Some evil genius (probably Gauss or Euler) decided to define a function over finite fields called the trace:

$$ \operatorname{Tr}_{F/K}(u) = u + u^p + u^{p^2} + u^{p^3} + \ldots + u^{p^{d-1}} $$

Here \( u \) is an element of some finite field \( F = GF(p^N) \) with characteristic \( p \), and \( K = GF(p^k) \) is a subfield of \( F \) with \( N = dk \). Then it turns out that \( y = \operatorname{Tr}_{F/K}(u) \) is always an element of the subfield \( K \).

(I’m just listing the definition, it’s a really abstract one to wrap your head around, let alone understand why it works. Essentially the trace is a way of projecting elements in a linear fashion from “\( N \)-dimensional space” to “\( k \)-dimensional space”, sort of like a photograph takes three-dimensional objects and projects them onto a two-dimensional image.)

If \( K \) is the base field \( GF(p) \) then you can drop the \( K \) and just say that

$$ \operatorname{Tr}_F(u) = u + u^p + u^{p^2} + u^{p^3} + \ldots + u^{p^{N-1}} $$

Since all we really care about is \( F=GF(2^N) \) in this series of articles,

$$ \operatorname{Tr}_F(u) = u + u^2 + u^4 + u^8 + \ldots + u^{2^{N-1}} $$

And we’ll get lazy and sloppy and just call it \( \operatorname{Tr}(u) \) rather than \( \operatorname{Tr} _ F(u) \). Let’s calculate this for a couple of examples; in \( H _ {25} \) since \( N=5 \) we have \( \operatorname{Tr}(u) = u + u^2 + u^4 + u^8 + u^{16} \).

def trace1(u):
    """ direct definition of trace for GF(2^n)"""
    result = 0
    for k in xrange(u.field.degree):
        result += u
        u = u * u
    return result.coeffs

for k in xrange(32):
    u = e25 << k
    print '%2d %s %d' % (k, u, trace1(u))
 0 GF2Element(0b00001,0x25) 1
 1 GF2Element(0b00010,0x25) 0
 2 GF2Element(0b00100,0x25) 0
 3 GF2Element(0b01000,0x25) 1
 4 GF2Element(0b10000,0x25) 0
 5 GF2Element(0b00101,0x25) 1
 6 GF2Element(0b01010,0x25) 1
 7 GF2Element(0b10100,0x25) 0
 8 GF2Element(0b01101,0x25) 0
 9 GF2Element(0b11010,0x25) 1
10 GF2Element(0b10001,0x25) 1
11 GF2Element(0b00111,0x25) 1
12 GF2Element(0b01110,0x25) 1
13 GF2Element(0b11100,0x25) 1
14 GF2Element(0b11101,0x25) 0
15 GF2Element(0b11111,0x25) 0
16 GF2Element(0b11011,0x25) 0
17 GF2Element(0b10011,0x25) 1
18 GF2Element(0b00011,0x25) 1
19 GF2Element(0b00110,0x25) 0
20 GF2Element(0b01100,0x25) 1
21 GF2Element(0b11000,0x25) 1
22 GF2Element(0b10101,0x25) 1
23 GF2Element(0b01111,0x25) 0
24 GF2Element(0b11110,0x25) 1
25 GF2Element(0b11001,0x25) 0
26 GF2Element(0b10111,0x25) 1
27 GF2Element(0b01011,0x25) 0
28 GF2Element(0b10110,0x25) 0
29 GF2Element(0b01001,0x25) 0
30 GF2Element(0b10010,0x25) 0
31 GF2Element(0b00001,0x25) 1

See, we always get a 0 or 1; in fact it happens to be the same LFSR sequence you get from \( H_{25} \) itself, but shifted by some number of counts.

We can generalize and compute the trace parity \( \operatorname{Tr}(wu) \) using some constant element \( w \) which is the trace pattern. (I’m using the terminology in Daemen and Rijmen; I haven’t seen it elsewhere, but it’s catchy.) If you remember from the last article, multiplication by a field element is the same as a timeshift: if \( u=x^k \) and \( w=x^d \), then \( wu = x^{k+d} \) and multiplication by \( w \) is identical to shifting by \( d \) samples.

The trace itself has a few other elegant properties.

One is linearity: \( \operatorname{Tr}(ax) + \operatorname{Tr}(by) = \operatorname{Tr}(ax+by) \) as long as \( x, y \in GF(p^N) \) and \( a,b \in GF(p) \). For our use of \( GF(2) \) the only choices for \( a \) and \( b \) are \( \{0,1\} \), so the only nontrivial statement of linearity is \( \operatorname{Tr}(x) + \operatorname{Tr}(y) = \operatorname{Tr}(x+y) \).

Also, the trace itself is invariant under decimation by a power of two:

$$\begin{align} \operatorname{Tr}(u^2) &= u^2 + (u^2)^2 + (u^2)^4 + \ldots + (u^2)^{2^{N-2}} + (u^2)^{2^{N-1}} \cr &= u^2 + u^4 + u^8 + \ldots + u^{2^{N-1}} + u \cr &= \operatorname{Tr}(u) \end{align}$$

with the identity \( (u^2)^{2^{N-1}} = u^{2^N} = u \) coming from the fact that \( u \), like all elements of \( F \), has an order that divides \( 2^N-1 \) so \( u^{{2^N}-1} = 1 \) and \( u^{2^N} = u \). For the trace computation, this essentially gives us a cyclic shift of all the \( u^{2^k} \) elements, leaving their net sum unchanged. So I can create a sequence that remains unchanged under decimation of 2:

for k in xrange(32):
    print '%2d %d %d %d %d' % (
        k, 
        trace1(e25<<k), 
        trace1(e25<<(2*k)),
        trace1(e25<<(4*k)),
        trace1(e25<<(8*k))
    )
 0 1 1 1 1
 1 0 0 0 0
 2 0 0 0 0
 3 1 1 1 1
 4 0 0 0 0
 5 1 1 1 1
 6 1 1 1 1
 7 0 0 0 0
 8 0 0 0 0
 9 1 1 1 1
10 1 1 1 1
11 1 1 1 1
12 1 1 1 1
13 1 1 1 1
14 0 0 0 0
15 0 0 0 0
16 0 0 0 0
17 1 1 1 1
18 1 1 1 1
19 0 0 0 0
20 1 1 1 1
21 1 1 1 1
22 1 1 1 1
23 0 0 0 0
24 1 1 1 1
25 0 0 0 0
26 1 1 1 1
27 0 0 0 0
28 0 0 0 0
29 0 0 0 0
30 0 0 0 0
31 1 1 1 1

Creepy! You get rid of every other value and you’re left with the original sequence. Take only every fourth value and you’re left with the original sequence. It’s like the Lernean Hydra — cut off a bunch of heads and it still remains.

Now, there’s a way to find the trace pattern \( w \) that corresponds with a given sequence.

Let’s look at the LFSR output \( y[k] \) which is the coefficient \( S_{N-1} \) of \( x^k = S_{N-1}x^{N-1} + S_{N-2}x^{N-2} + \ldots + S_2x^2 + S_1x + S_0 \). It’s possible to show that \( y[k] = f(x^k) \) is a linear function of \( x^k \), and it turns out that because it’s a linear function we can write it as \( y[k] = f(x^k) = \operatorname{Tr}(wx^k) \) for some fixed \( w \); we just have to figure out what this value of \( w \) is. (See Lidl & Niederreiter or Daemen & Rijmen for a proof that all linear functions with domain \( GF(p^N) \) and output \( GF(p) \) can be written in this form.)

The values \( y[0], y[1], y[2], \ldots, y[N-1] \) correspond to \( S[0] = x^0 = 1 \). If we take the values \( y[0] = y[d], y[2] = y[d+1], y[4]=y[d+2] \ldots, y[2N-2]=y[N-1+d] \) and use the state recovery technique we mentioned in the last article, then we can figure out some corresponding state \( x^d \), even if we don’t know what \( d \) is. Compute the discrete logarithm if you want to know.

What we are alleging is that \( y[2k] = f(x^{2k}) = y[k+d] = f(x^kx^d) \) for all values of \( k \), not just from \( k=0, 1, \ldots N-1 \).

Because \( \operatorname{Tr}(u^2) = \operatorname{Tr}(u) \),

$$y[k+d] = f(x^kx^d) = \operatorname{Tr}(wx^kx^d) = \operatorname{Tr}(w^2x^{2k}x^{2d}) = f(wx^{2k}x^{2d})$$

If \( wx^{2d} = 1 \) then we have \( y[k+d] = f(x^{2k}) = y[2k] \) for all \( k \), so this means that we can calculate \( w = (x^{2d})^{-1} = ((x^d)^2)^{-1} \).

Essentially here’s all we need to do in order to calculate \( w \):

  • get a collection of the first \( N \) bits if we look at every other output bit of the original LFSR, in other words \( y[0], y[2], y[4], \ldots, y[2N-4], y[2N-2] \)
  • find the state \( v_2 = S[d] = x^d \) from these output bits using state recovery (we don’t need to find the value of \( d \))
  • square it
  • take the multiplicative inverse

Let’s do it!

from libgf2.util import state_from_reversed_output

e25 = H25.wrap(1)
y2 = [highbit(e25 << (2*k)) for k in xrange(H25.degree)]
v2 = H25.wrap(state_from_reversed_output(H25.coeffs, y2[::-1]))
w = (v2*v2).inv
print "w=%s" % (w,)

for k in xrange(31):
    xk = e25 << k
    print '%2d %d %d' % (k, highbit(xk), trace1(w*xk))
w=GF2Element(0b01011,0x25)
 0 0 0
 1 0 0
 2 0 0
 3 0 0
 4 1 1
 5 0 0
 6 0 0
 7 1 1
 8 0 0
 9 1 1
10 1 1
11 0 0
12 0 0
13 1 1
14 1 1
15 1 1
16 1 1
17 1 1
18 0 0
19 0 0
20 0 0
21 1 1
22 1 1
23 0 0
24 1 1
25 1 1
26 1 1
27 0 0
28 1 1
29 0 0
30 1 1

Want to decimate the LFSR sequence by \( j=2 \)? No problem! Just use the high bit of \( v_2x^k \).

What about decimating by \( j=4 \) or some other power \( j=2^b \)? This is a bit trickier; we need some multiplier \( v_j \) such that \( (wv_j)^j = w \), and because the trace is invariant under raising its argument to a power of two (\( \operatorname{Tr}(u^j) = \operatorname{Tr}(u) \)),

$$f(x^{jk}) = \operatorname{Tr}(wx^{jk}) = \operatorname{Tr}((wv_j)^jx^{jk}) = \operatorname{Tr}(wv_jx^k) = f(v_jx^k)$$

which will happen if \( (wv_j)^j = w \implies wv_j = w^{1/j} \implies v_j = w^{-1}w^{1/j} \).

Taking the \( j \)th root may seem tricky, but since \( j \) is a power of two, and raising any field element to the \( 2^N \) order leaves its value unchanged, \( w^{1/j} = w^\frac{2^N}{j} \) and we can compute \( v_j = w^a \) with \( a=\frac{2^N}{j} - 1 \).

v2 = w**(16-1)
v4 = w**(8-1)
v8 = w**(4-1)
v16 = w**(2-1)

for k in xrange(32):
    xk = e25 << k
    print '%2d | %d %d | %d %d | %d %d | %d %d' % (
            k, 
            highbit(e25 << (2*k)), highbit(xk*v2),
            highbit(e25 << (4*k)), highbit(xk*v4),
            highbit(e25 << (8*k)), highbit(xk*v8),
            highbit(e25 << (16*k)), highbit(xk*v16),
    )
 0 | 0 0 | 0 0 | 0 0 | 0 0
 1 | 0 0 | 1 1 | 0 0 | 1 1
 2 | 1 1 | 0 0 | 1 1 | 0 0
 3 | 0 0 | 0 0 | 1 1 | 1 1
 4 | 0 0 | 1 1 | 0 0 | 0 0
 5 | 1 1 | 0 0 | 1 1 | 0 0
 6 | 0 0 | 1 1 | 1 1 | 0 0
 7 | 1 1 | 1 1 | 1 1 | 0 0
 8 | 1 1 | 0 0 | 0 0 | 1 1
 9 | 0 0 | 0 0 | 1 1 | 0 0
10 | 0 0 | 1 1 | 0 0 | 0 0
11 | 1 1 | 1 1 | 1 1 | 1 1
12 | 1 1 | 1 1 | 0 0 | 0 0
13 | 1 1 | 1 1 | 0 0 | 1 1
14 | 1 1 | 1 1 | 0 0 | 1 1
15 | 1 1 | 0 0 | 0 0 | 0 0
16 | 0 0 | 0 0 | 1 1 | 0 0
17 | 0 0 | 0 0 | 0 0 | 1 1
18 | 0 0 | 1 1 | 0 0 | 1 1
19 | 1 1 | 1 1 | 1 1 | 1 1
20 | 1 1 | 0 0 | 0 0 | 1 1
21 | 0 0 | 1 1 | 1 1 | 1 1
22 | 1 1 | 1 1 | 1 1 | 0 0
23 | 1 1 | 1 1 | 0 0 | 0 0
24 | 1 1 | 0 0 | 0 0 | 0 0
25 | 0 0 | 1 1 | 1 1 | 1 1
26 | 1 1 | 0 0 | 1 1 | 1 1
27 | 0 0 | 1 1 | 1 1 | 0 0
28 | 1 1 | 0 0 | 1 1 | 1 1
29 | 0 0 | 0 0 | 1 1 | 1 1
30 | 0 0 | 0 0 | 0 0 | 1 1
31 | 0 0 | 0 0 | 0 0 | 0 0

What about decimation ratios \( j \) that are not powers of two? Then \( \operatorname{Tr} _ F(x^{jk}) = \operatorname{Tr} _ {F’}(wx^k) \) for some other field \( F’ \), but there’s not really a mathematical shortcut for determining \( F’ \) and \( w \) from \( F \) and \( j \), at least not that I’m aware of. The easiest way to determine these is to generate \( 2N \) successive values of the trace from \( k=0 \) to \( k=2N-1 \), then run Berlekamp-Massey to determine the polynomial of the field \( F’ \), then use state recovery to determine \( w \).

There’s one more special case, and that is if you use a decimation ratio of \( -1 \), which essentially means going through the sequence backwards. In this case the polyomials are reflections of each other. With fields isomorphic to \( GF(2^5) \) this corresponds to decimation ratios of \( -1 \equiv 30 \):

for p in [0x25, 0x29, 0x2f, 0x37, 0x3b, 0x3d]:
    e = GF2Element(1, p)
    print ""
    for j in [1,30]:
        poly, deg = berlekamp_massey([highbit(e << (j*k)) for k in xrange(10)])
        print 'minimal polynomial for decimation ratio j={0:2d} is {1:02x}={1:b}'.format(j, poly)
minimal polynomial for decimation ratio j= 1 is 25=100101
minimal polynomial for decimation ratio j=30 is 29=101001

minimal polynomial for decimation ratio j= 1 is 29=101001
minimal polynomial for decimation ratio j=30 is 25=100101

minimal polynomial for decimation ratio j= 1 is 2f=101111
minimal polynomial for decimation ratio j=30 is 3d=111101

minimal polynomial for decimation ratio j= 1 is 37=110111
minimal polynomial for decimation ratio j=30 is 3b=111011

minimal polynomial for decimation ratio j= 1 is 3b=111011
minimal polynomial for decimation ratio j=30 is 37=110111

minimal polynomial for decimation ratio j= 1 is 3d=111101
minimal polynomial for decimation ratio j=30 is 2f=101111

Here each pair of polynomials consists of reflections of each other, for example \( H_{25} = x^5 + x^2 + 1 \) and \( H_{29} = x^5 + x^3 + 1 \).

Trace computation in practice

Mathematically, we have stated that the following three values are equal:

  • the LFSR outputs \( y[k] \) = high bit of LFSR state \( S[k] \)
  • the coefficient of \( x^{N-1} \) in \( x^k \bmod p(x) \), which is really the same thing as LFSR state \( S[k] \) but in polynomial representation
  • the trace calculation \( \operatorname{Tr}(wx^k) \), where we showed how to calculate \( w \) from the LFSR outputs \( y[0], y[2], y[4], \ldots, y[2N-4], y[2N-2] \)

If you really want to compute the trace \( \operatorname{Tr}(u) \) for some reason, then there are a few ways to do it:

  • the mathematical definition \( \operatorname{Tr}(u) = u + u^2 + u^4 + \ldots + u^{2^{N-1}} \)
  • the coefficient of \( x^{N-1} \) (= high bit of LFSR state) in \( w^{-1}u \)

But there’s an easier way, which is to compute the parity of a certain subset of bits in the polynomial bit vector representation of \( u \). Mathematically this is just \( \sum\limits_{k=0}^{N-1}m_ks_k \bmod 2 \) for \( u=\sum\limits_{k=0}^{N-1}s_kx^k \) and some set of coefficients \( m_k \) that represent the mask bits.

We can figure out the mask bits \( m_k \) by using polynomial basis values \( u=x^k \) for \( k=0, 1, 2, \ldots, N-1 \) and computing the trace for each of them. For example, if \( u=x^3 \) then the trace is equal to \( m_3s_3 = m_3 \) since \( s_3 = 1 \) and all the other \( s_k \) coefficients are zero.

Here’s a systematic way to do it:

from libgf2.util import parity

def trace_helper(field):
    mask = 0
    n = field.degree
    e = field.wrap(1)
    y2 = [highbit(e << (2*k)) for k in xrange(n)]
    v2 = field.wrap(state_from_reversed_output(field.coeffs, y2[::-1]))
    winv = (v2*v2)
    u = winv
    for k in xrange(n):
        # loop invariant: u = w^(-1)*x^k
        mask |= (u.coeffs >> (n-1)) << k
        u <<= 1
    return mask, winv

def show_trace_equivalents(field):
    mask, winv = trace_helper(field)
    print "mask=%s, w^(-1)=%s" % (mask, winv)
    for x in xrange(1 << field.degree):
        u = field.wrap(x)
        print '%2d %d %d %d' % (
                x, 
                trace1(u),                              # direct definition 
                (u*winv).coeffs >> (field.degree - 1),  # high bit of u*w^(-1)
                parity(x & mask))                       # parity method
        
show_trace_equivalents(H25)
mask=9, w^(-1)=GF2Element(0b10000,0x25)
 0 0 0 0
 1 1 1 1
 2 0 0 0
 3 1 1 1
 4 0 0 0
 5 1 1 1
 6 0 0 0
 7 1 1 1
 8 1 1 1
 9 0 0 0
10 1 1 1
11 0 0 0
12 1 1 1
13 0 0 0
14 1 1 1
15 0 0 0
16 0 0 0
17 1 1 1
18 0 0 0
19 1 1 1
20 0 0 0
21 1 1 1
22 0 0 0
23 1 1 1
24 1 1 1
25 0 0 0
26 1 1 1
27 0 0 0
28 1 1 1
29 0 0 0
30 1 1 1
31 0 0 0
H73 = GF2QuotientRing(0x73)
show_trace_equivalents(H73)
mask=22, w^(-1)=GF2Element(0b010001,0x73)
 0 0 0 0
 1 0 0 0
 2 1 1 1
 3 1 1 1
 4 1 1 1
 5 1 1 1
 6 0 0 0
 7 0 0 0
 8 0 0 0
 9 0 0 0
10 1 1 1
11 1 1 1
12 1 1 1
13 1 1 1
14 0 0 0
15 0 0 0
16 1 1 1
17 1 1 1
18 0 0 0
19 0 0 0
20 0 0 0
21 0 0 0
22 1 1 1
23 1 1 1
24 1 1 1
25 1 1 1
26 0 0 0
27 0 0 0
28 0 0 0
29 0 0 0
30 1 1 1
31 1 1 1
32 0 0 0
33 0 0 0
34 1 1 1
35 1 1 1
36 1 1 1
37 1 1 1
38 0 0 0
39 0 0 0
40 0 0 0
41 0 0 0
42 1 1 1
43 1 1 1
44 1 1 1
45 1 1 1
46 0 0 0
47 0 0 0
48 1 1 1
49 1 1 1
50 0 0 0
51 0 0 0
52 0 0 0
53 0 0 0
54 1 1 1
55 1 1 1
56 1 1 1
57 1 1 1
58 0 0 0
59 0 0 0
60 0 0 0
61 0 0 0
62 1 1 1
63 1 1 1

Traces and Decimation in libgf2

In libgf2, the Python module for computation in \( GF(2^n) \) which no one uses, trace computations can be performed in a couple of ways.

As far as \( \operatorname{Tr}(u) \) itself, you can calculate it with element.trace:

e = GF2Element(1, 0x25)
for k in xrange(31):
    print "Tr(x^%2d) = %x = %x" % (k, trace1(e << k), (e << k).trace)
Tr(x^ 0) = 1 = 1
Tr(x^ 1) = 0 = 0
Tr(x^ 2) = 0 = 0
Tr(x^ 3) = 1 = 1
Tr(x^ 4) = 0 = 0
Tr(x^ 5) = 1 = 1
Tr(x^ 6) = 1 = 1
Tr(x^ 7) = 0 = 0
Tr(x^ 8) = 0 = 0
Tr(x^ 9) = 1 = 1
Tr(x^10) = 1 = 1
Tr(x^11) = 1 = 1
Tr(x^12) = 1 = 1
Tr(x^13) = 1 = 1
Tr(x^14) = 0 = 0
Tr(x^15) = 0 = 0
Tr(x^16) = 0 = 0
Tr(x^17) = 1 = 1
Tr(x^18) = 1 = 1
Tr(x^19) = 0 = 0
Tr(x^20) = 1 = 1
Tr(x^21) = 1 = 1
Tr(x^22) = 1 = 1
Tr(x^23) = 0 = 0
Tr(x^24) = 1 = 1
Tr(x^25) = 0 = 0
Tr(x^26) = 1 = 1
Tr(x^27) = 0 = 0
Tr(x^28) = 0 = 0
Tr(x^29) = 0 = 0
Tr(x^30) = 0 = 0

Trace parity, namely \( \operatorname{Tr}(wu) \) for some fixed trace pattern \( w \), can be calculated using GF2TracePattern. This implementation uses parity and a bitmask rather than the direct definition of computing trace, for efficiency. Instances of GF2TracePattern can be created in three different ways:

  • from the trace pattern \( w \), by calling GF2TracePattern.from_pattern(field, w)
  • from the bitmask \( m \), by calling GF2TracePattern.from_mask(field, m)
  • from the equivalent initial LFSR state \( s \), by calling GF2TracePattern.from_lfsr_initial_state(field, s); this is kind of odd, but basically the value \( s \) is such that the trace parity \( \operatorname{Tr}(wx^k) = y[k] \) where \( y[k] \) is the high bit of a Galois LFSR with initial state \( s \).

Each GF2TracePattern instance has pattern, mask, and lfsr_initial_state properties which can be accessed, and the trace itself is calculated by using the GF2TracePattern like a function:

from libgf2.gf2 import GF2TracePattern

e = H25.wrap(1)
trp1 = GF2TracePattern.from_lfsr_initial_state(H25, e)
print "pattern w=0x%x, mask=0x%x, lfsr_initial_state=0x%x" % (
    trp1.pattern,
    trp1.mask,
    trp1.lfsr_initial_state)
w = trp1.pattern
for k in xrange(31):
    u = e << k
    print 'Tr(w*x^%2d) = %d = %d' % (k, trace1(w*u), trp1(u))
pattern w=0xb, mask=0x10, lfsr_initial_state=0x1
Tr(w*x^ 0) = 0 = 0
Tr(w*x^ 1) = 0 = 0
Tr(w*x^ 2) = 0 = 0
Tr(w*x^ 3) = 0 = 0
Tr(w*x^ 4) = 1 = 1
Tr(w*x^ 5) = 0 = 0
Tr(w*x^ 6) = 0 = 0
Tr(w*x^ 7) = 1 = 1
Tr(w*x^ 8) = 0 = 0
Tr(w*x^ 9) = 1 = 1
Tr(w*x^10) = 1 = 1
Tr(w*x^11) = 0 = 0
Tr(w*x^12) = 0 = 0
Tr(w*x^13) = 1 = 1
Tr(w*x^14) = 1 = 1
Tr(w*x^15) = 1 = 1
Tr(w*x^16) = 1 = 1
Tr(w*x^17) = 1 = 1
Tr(w*x^18) = 0 = 0
Tr(w*x^19) = 0 = 0
Tr(w*x^20) = 0 = 0
Tr(w*x^21) = 1 = 1
Tr(w*x^22) = 1 = 1
Tr(w*x^23) = 0 = 0
Tr(w*x^24) = 1 = 1
Tr(w*x^25) = 1 = 1
Tr(w*x^26) = 1 = 1
Tr(w*x^27) = 0 = 0
Tr(w*x^28) = 1 = 1
Tr(w*x^29) = 0 = 0
Tr(w*x^30) = 1 = 1

This is also useful for seeing what kind of mask is necessary to achieve various time delays:

def powers_of_x(field):
    e = tr.field.wrap(1)
    while True:
        yield e
        e <<= 1

def get_bits_as_string(tr,m):
    egen = powers_of_x(tr.field)
    return ''.join('1' if tr(egen.next()) else '0' for _ in xrange(m))

for d in xrange(31):
    tr = GF2TracePattern.from_lfsr_initial_state(H25, H25.wrap(1) >> d)
    print "delay={0:2d}, s={1:02x}, mask={2:05b}".format(
        d, tr.lfsr_initial_state, tr.mask), get_bits_as_string(tr,60)
delay= 0, s=01, mask=10000 000010010110011111000110111010100001001011001111100011011101
delay= 1, s=12, mask=00001 100001001011001111100011011101010000100101100111110001101110
delay= 2, s=09, mask=00010 010000100101100111110001101110101000010010110011111000110111
delay= 3, s=16, mask=00101 101000010010110011111000110111010100001001011001111100011011
delay= 4, s=0b, mask=01010 010100001001011001111100011011101010000100101100111110001101
delay= 5, s=17, mask=10101 101010000100101100111110001101110101000010010110011111000110
delay= 6, s=19, mask=01011 110101000010010110011111000110111010100001001011001111100011
delay= 7, s=1e, mask=10111 111010100001001011001111100011011101010000100101100111110001
delay= 8, s=0f, mask=01110 011101010000100101100111110001101110101000010010110011111000
delay= 9, s=15, mask=11101 101110101000010010110011111000110111010100001001011001111100
delay=10, s=18, mask=11011 110111010100001001011001111100011011101010000100101100111110
delay=11, s=0c, mask=10110 011011101010000100101100111110001101110101000010010110011111
delay=12, s=06, mask=01100 001101110101000010010110011111000110111010100001001011001111
delay=13, s=03, mask=11000 000110111010100001001011001111100011011101010000100101100111
delay=14, s=13, mask=10001 100011011101010000100101100111110001101110101000010010110011
delay=15, s=1b, mask=00011 110001101110101000010010110011111000110111010100001001011001
delay=16, s=1f, mask=00111 111000110111010100001001011001111100011011101010000100101100
delay=17, s=1d, mask=01111 111100011011101010000100101100111110001101110101000010010110
delay=18, s=1c, mask=11111 111110001101110101000010010110011111000110111010100001001011
delay=19, s=0e, mask=11110 011111000110111010100001001011001111100011011101010000100101
delay=20, s=07, mask=11100 001111100011011101010000100101100111110001101110101000010010
delay=21, s=11, mask=11001 100111110001101110101000010010110011111000110111010100001001
delay=22, s=1a, mask=10011 110011111000110111010100001001011001111100011011101010000100
delay=23, s=0d, mask=00110 011001111100011011101010000100101100111110001101110101000010
delay=24, s=14, mask=01101 101100111110001101110101000010010110011111000110111010100001
delay=25, s=0a, mask=11010 010110011111000110111010100001001011001111100011011101010000
delay=26, s=05, mask=10100 001011001111100011011101010000100101100111110001101110101000
delay=27, s=10, mask=01001 100101100111110001101110101000010010110011111000110111010100
delay=28, s=08, mask=10010 010010110011111000110111010100001001011001111100011011101010
delay=29, s=04, mask=00100 001001011001111100011011101010000100101100111110001101110101
delay=30, s=02, mask=01000 000100101100111110001101110101000010010110011111000110111010

Under successive time delays, the mask bits appear to evolve as the state of a Fibonacci LFSR — that is, if \( M[d] \) is the bitmask for delay \( d \), then the rightmost bit of \( M[d+1] \) is an XOR of a selected subset of bits \( M[d] \), and the other bits are just a left shift from \( M[d] \). This is not hard to show; let’s take an example:

Suppose we have \( p_{25}(x) = x^5 + x^2 + 1 \) and our mask \( M \) is some collection of mask bits \( m_4, m_3, m_2, m_1, m_0 \), so if the LFSR state \( S[k]= x^k = s_4[k]x^4 + s_3[k]x^3 + s_2[k]x^2 + s_1[k]x + s_0[k] \), then the output sequence is \( y[k] = m_4s_4[k] + m_3s_3[k] + m_2s_2[k] + m_1s_1[k] + m_0s_0[k] \), by definition: we take the bitwise AND of the mask bits with the state bits, and compute the parity of the resulting value.

Now look at a delayed version of the output: \( y[k-1] = m_4s_4[k-1] + m_3s_3[k-1] + m_2s_2[k-1] + m_1s_1[k-1] + m_0s_0[k-1] \); we just use the mask bits \( M \) with the previous state bits \( S[k-1] \). But there’s another way to compute \( y[k-1] \) from the state bits at this state \( S[k] \) and some new mask coefficients \( m_i’ \): \( y[k-1] = m_4’s_4[k] + m_3’s_3[k] + m_2’s_2[k] + m_1’s_1[k] + m_0’s_0[k] \). And we know how the state bits evolve over time:

$$\begin{aligned} s_4[k] &= s_3[k-1] + a_4s_4[k-1]\cr s_3[k] &= s_2[k-1] + a_3s_4[k-1]\cr s_2[k] &= s_1[k-1] + a_2s_4[k-1]\cr s_1[k] &= s_0[k-1] + a_1s_4[k-1]\cr s_0[k] &= \phantom{s_1[k-1] + {} } a_0s_4[k-1]\cr \end{aligned}$$

with \( p(x) = x^5 + a_4x^4 + a_3x^3 + a_2x^2 + a_1x + a_0 \); we can invert the state evolution equations to determine

$$\begin{aligned} s_4[k-1] &= \phantom{s_1[k] + {}} a_0s_0[k]\cr s_3[k-1] &= s_4[k] + a_4s_0[k]\cr s_2[k-1] &= s_3[k] + a_3s_0[k]\cr s_1[k-1] &= s_2[k] + a_2s_0[k]\cr s_0[k-1] &= s_1[k] + a_1s_0[k]\cr \end{aligned}$$

Combining together these equations with the two methods of computing \( y[k-1] \) and matching coefficients of \( s_0[k], s_1[k], \ldots \), we get

$$\begin{aligned} y[k-1] &=m_4’s_4[k] + m_3’s_3[k] + m_2’s_2[k] + m_1’s_1[k] + m_0’s_0[k] \cr &= m_4a_0s_0[k] + m_3s_4[k] + m_3a_4s_0[k] + m_2s_3[k]+m_2a_3s_0[k] \cr &\phantom{=} + m_1s_2[k] + m_1a_2s_0[k] + m_0s_1[k] +m_0a_1s_0[k] \end{aligned}$$

or

$$\begin{aligned} m_4’ &= m_3 \cr m_3’ &= m_2 \cr m_2’ &= m_1 \cr m_1’ &= m_0 \cr m_0’ &= m_4a_0 + m_3a_4 + m_2a_3 + m_1a_2 + m_0a_1 \end{aligned}$$

In other words, the mask \( M’ \) with delay \( d+1 \) is equal to the mask \( M \) with delay \( d \), shifted left and with its rightmost bit computed from the mask bits and polynomial bits. There you go.

Enough about masks and time delays.

Back to decimation: the libgf2 module also has a function libgf2.lfsr.decimate for computing the field representation that corresponds to decimated LFSR output from a given field representation under decimation ratio \( j \) (again, note that the quotient rings are invariant under decimation ratios which are powers of two):

from libgf2.lfsr import decimate

for j in xrange(32):
    print 'j=%2d' % j, decimate(H25, j)
j= 0 GF2QuotientRing(0b1)
j= 1 GF2QuotientRing(0b100101)
j= 2 GF2QuotientRing(0b100101)
j= 3 GF2QuotientRing(0b111101)
j= 4 GF2QuotientRing(0b100101)
j= 5 GF2QuotientRing(0b110111)
j= 6 GF2QuotientRing(0b111101)
j= 7 GF2QuotientRing(0b101111)
j= 8 GF2QuotientRing(0b100101)
j= 9 GF2QuotientRing(0b110111)
j=10 GF2QuotientRing(0b110111)
j=11 GF2QuotientRing(0b111011)
j=12 GF2QuotientRing(0b111101)
j=13 GF2QuotientRing(0b111011)
j=14 GF2QuotientRing(0b101111)
j=15 GF2QuotientRing(0b101001)
j=16 GF2QuotientRing(0b100101)
j=17 GF2QuotientRing(0b111101)
j=18 GF2QuotientRing(0b110111)
j=19 GF2QuotientRing(0b101111)
j=20 GF2QuotientRing(0b110111)
j=21 GF2QuotientRing(0b111011)
j=22 GF2QuotientRing(0b111011)
j=23 GF2QuotientRing(0b101001)
j=24 GF2QuotientRing(0b111101)
j=25 GF2QuotientRing(0b101111)
j=26 GF2QuotientRing(0b111011)
j=27 GF2QuotientRing(0b101001)
j=28 GF2QuotientRing(0b101111)
j=29 GF2QuotientRing(0b101001)
j=30 GF2QuotientRing(0b101001)
j=31 GF2QuotientRing(0b1)

Computing the Decimation Ratio

If you can determine quotient rings from an LFSR that is equivalent to the decimated output from another LFSR, the question arises:

  • Suppose I have two LFSRs LFSR1 and LFSR2 with characteristic polynomials \( p_1(x) \) and \( p_2(x) \) and output bits \( y_1[k] \) and \( y_2[k] \)
  • Suppose furthermore that I know that \( y_2[k] = y_1[jk+b] \) for some unknown values \( j \) and \( b \); that is, LFSR2 output is the same as the decimation of the LFSR1 output for some unknown decimation ratio \( j \) with an unknown time offset.
  • How can I compute the decimation ratio \( j \)?

This was something that stumped me for over a month, and I finally had to ask a question on mathoverflow.net to get some help. Here’s the answer in a nutshell:

  • Field \( F_1 = GF(2)[x]/p_1(x) \) has the property that \( p_1(x) = 0 \), that is, \( x \) is a root of \( p_1(x) \) in \( F_1[x] \).
  • Field \( F_2 = GF(2)[x]/p_2(x) \) has the property that \( p_2(x) = 0 \), that is, \( x \) is a root of \( p_2(x) \) in \( F_2[x] \). But it is also true that \( x^j \) is a root of \( p_2(y) \) in \( F_1[y] \). This is a weird statement and it doesn’t make much sense unless you’re familiar with finite field theory and the implications of minimal polynomials; I’m about at my limit here. What are these \( F_1[y] \) things?

Lidl & Niederreiter has a chapter on factorization of polynomials over finite fields. This book is full of fast-paced theory with proofs, so lots of it goes over my head, but I did learn a number of things from it, some in terminology:

\( F[x] \) represents a polynomial ring with coefficients in \( F \). For example, \( GF(2)[x] \) is the ring of all polynomials \( p(x) \) with coefficients in \( GF(2) \), which includes polynomials like \( x \) and \( x^7 + x^6 + x^3 + 1 \) and \( x^{1234} + x^{44} \). Basically all the coefficients are 0 and 1, and the arithmetic to compute sums and products is in \( GF(2) \), so for example \( (x^3 + x)(x^2 + 1) \) = \( x^5 + x^3 + x^3 + x = x^5 + x \).

\( F_1[y] \) would be a polynomial ring with polynomials \( p(y) \) with coefficients in \( F_1 \), and our \( F_1 \) above is a quotient ring with elements that are themselves polynomials. Let’s look at a simple example, \( F_1 = GF(2)[x]/(x^3 + x + 1) \) which has 8 elements, \( 0, 1, x, x+1, x^2, x^2+1, x^2+x, x^2+x+1 \). Then \( F_1[y] \) would contain any polynomial of \( y \) with coefficients that are one of these 8 elements, for example \( p(y) = xy^3 + (x^2+1)y + (x^2 + x) \) has coefficients \( a_3 = x, a_2 = 0, a_1 = x^2+1, a_0 = x^2 + x \).

There are two variables here, \( x \) and \( y \), and I guess that’s what makes it confusing. Perhaps it becomes less so if you just represent \( F_1 \) by bit vectors, for example \( x \equiv {\tt 010}, x^2+1 \equiv {\tt 101} \), and so on, so the \( p(y) \) in the last paragraph could be written as \( p(y) = {\tt 010}y^3 + {\tt 101}y + {\tt 110} \).

In any case, there are algorithms for determining roots of polynomials over finite fields, and they use some weird mathematical identities to help compute the roots efficiently. By roots, for example, I would mean some value \( y \) in \( F_1 \) that makes \( p(y) = {\tt 010}y^3 + {\tt 101}y + {\tt 110} = 0 \).

One of the methods that works well for fields \( GF(p^N) \) of small characteristic \( p \) is the identity \( x^{p^N} - x = \prod\limits_{c\in GF(p)} (S(x) - c) \) where \( S(x) = x + x^p + x^{p^2} + x^{p^3} + \ldots + x^{p^{N-1}} \). For \( p=2 \) this is just \( x^{2^N} - x = S(x)(S(x)-1) \) with \( S(x) = x + x^2 + x^4 + x^8 + \ldots + x^{2^{N-1}} \).

There’s a way to use this identity in our decimation problem to write \( p(y) = \gcd(p(y),S(x^ky))\gcd(p(y),S(x^ky)-1) \). This splits a polynomial \( p(y) \) into two subfactors \( \gcd(p(y),S(x^ky)) \) and \( \gcd(p(y),S(x^ky)-1) \). You end up computing \( S(x^ky) \) modulo \( p(y) \) for various values of \( k \), first for \( k=0 \) and coming up with two factors of \( p(y) \), which you can then run through the algorithm again with \( k=1 \), and so on, dividing until you end up with a linear factor \( (y-y_1) \) of \( p(y) \) in \( F_1 \). The value of \( y_1 \) is the desired root, which can be written as \( y_1 = x^j \) and then you just take the logarithm of \( x^j \) to get \( j \).

So the solution works like this:

  • Find any root \( y_1 \in F_1 \) of \( p_2(y) \) in the polynomial ring \( F_1[y] \)
  • Take the discrete logarithm to determine \( j \) such that \( x^j = y_1 \)
  • Any decimation ratio in the cyclotomic coset of \( j \) is a valid decimation ratio.

If you like theory, read Lidl and Niederreiter; the appropriate section in the 1994 edition is section 4.3 (Calculation of Roots of Polynomials), equation 4.25, illustrated in example 4.16.

Otherwise you can forget how it works, and just use my libgf2.lfsr.undecimate function:

from libgf2.lfsr import undecimate

H211 = GF2QuotientRing(0x211)
e = H211.wrap(1)
test_j = 109
samples = [(e << (test_j * k)).coeffs >> 8 for k in xrange(18)]
p2,n = berlekamp_massey(samples)
print samples
H211_decimated = GF2QuotientRing(p2)
assert H211_decimated == decimate(H211, test_j)
print 'p2=0x%x' % p2
undecimate(H211_decimated, H211)
[0, 1, 1, 0, 0, 0, 1, 0, 1, 0, 1, 0, 0, 1, 0, 1, 0, 1]
p2=0x37f
109

This works with larger fields too, but be warned it can take quite a while.

import time

H48a = GF2QuotientRing(0x10000000000b7)
H48b = decimate(H48a, 123457)
print "decimated field characteristic polynomial = %x" % H48b.coeffs
t0 = time.time()
print "j=", undecimate(H48b, H48a)
t1 = time.time()
print "elapsed time: %.2f s" % (t1-t0)
decimated field characteristic polynomial = 110dccf2f72ab
j= 123457
elapsed time: 14.53 s

That’s pretty much it for today. There’s some interesting stuff in Lidl & Niederreiter on bases, which are ways of decomposing an element in a finite field into weighted sums of basis elements — much like decomposing waveforms into sums of sine waves — but I haven’t had any direct application of this yet, so I’ll leave it alone for now.

References

Acknowledgements

Thanks to Dr. Max Alekseyev for answering my stumper of a question about computing an unknown decimation ratio, and to Dr. Jyrki Lahtonen for answering several questions on math.stackexchange.com.

Wrapup

Today we looked at decimation of LFSR output by a fixed decimation ratio \( j \).

  • We showed that this always yields another sequence that can be generated from successive samples of another LFSR; the LFSR in question may or may not have the same characteristic polynomial as the original LFSR
  • Decimating by a power of two always yields a shifted version of the original sequence from the same LFSR
  • The trace \( \operatorname{Tr}(u) = u + u^2 + u^4 + \ldots + u^{2^{N-1}} \) can be used to help analyze decimation, as well as calculate arbitrary output samples from an LFSR
  • Decimation ratios that differ by a factor of two, modulo \( 2^N-1 \), form a cyclotomic coset (for example, \( 3, 6, 12, 24, 17 \) for \( N=5 \)), and all yield the same decimated LFSR output (with varying shifts) from the same characteristic polynomial
  • It is possible to determine the unknown decimation ratio \( j \) from two sequences, one of which is a decimation of the other, by using an algorithm to find the root of a polynomial and taking the discrete logarithm
  • The libgf2 module contains classes and functions for assisting with these computations, namely GF2TracePattern, decimate, and undecimate.

Aside from one more theoretical topic of interest — efficient generation of primitive polynomials of a given degree, which we’ll cover later — and some short asides, this ends our discussion of finite field theory. You can now relax.

Next time we will transition into the area of LFSR applications, starting with counters and encoders.


© 2017 Jason M. Sachs, all rights reserved.


Previous post by Jason Sachs:
   Linear Feedback Shift Registers for the Uninitiated, Part VIII: Matrix Methods and State Recovery
Next post by Jason Sachs:
   Linear Feedback Shift Registers for the Uninitiated, Part X: Counters and Encoders

Comments:

To post reply to a comment, click on the 'reply' button attached to each comment. To post a new comment (not a reply to a comment) check out the 'Write a Comment' tab at the top of the comments.

Registering will allow you to participate to the forums on ALL the related sites and give you access to all pdf downloads.

Sign up
or Sign in