Shifting Up A Gear

| 0 Comments

Last time[1] we took a look at linear feedback shift registers, or LFSRs, which produce pseudo random bits by shifting a buffer of them to the right, with the rightmost becoming both the output and the leftmost bit which being combined with a subset of the bits, known as taps, using an exclusive or operation as illustrated by figure 1.

Figure 1: A Linear Feedback Shift Register

Here the empty circles represent zeros, the filled circles ones, the \(\oplus\) symbols exclusive or operations and the arrows indicate the flow of the bits through the register and, by convention, the rightmost bit is labelled as bit one rather than bit zero.

We saw that despite looking completely different to the linear congruential generators that we covered a few posts ago[2] they are in fact rather similar to them, albeit involving multiplication and remainder operations involving polynomials rather than integers.
Unfortunately, another thing that LFSRs have in common with linear congruential generators is that it is extremely difficult to set them up so that they generate sequences of bits that are as close to random as possible. Thankfully, there are those who've gone to the trouble of finding decent choices for the taps for us[3] and so all that we need to concern ourselves with is implementing them, which we shall do in this post.

A Small LFSR

A third thing that LFSRs have in common with congruential generators is that the efficiency of their implementations depends upon the specific values that define them, in this case the number of bits and the positions of the taps.
Most trivially, if the number of bits is less than or equal to 32 then we can implement them with a single number, as is done by the smallRnd function given in listing 1.

Listing 1: smallRnd
function smallRnd(bits, taps, seed) {
 var state = ak.floor((Math.pow(2, bits)-1)*ak.ranfRnd(seed)())+1;
 var mask = smallMask(bits, taps);

 return function() {
  var b = state&0x01;
  state >>>= 1;
  if(b) state ^= mask;
  return b;
 };
}

First, we seed the bits of the generator's state using an ak.ranfRnd generator which, since it is randomly initialised if given an undefined seed, means that if seed is undefined the bits are seeded randomly.
Note that because a register full of zero bits will stay that way forever we must task care to avoid initialising the state with zero. Here we do this by exploiting the fact that, like Math.random, the generators created by ak.ranfRnd return a value greater than or equal to zero and strictly less than one and so if we multiply them by the largest unsigned integer with the given number of bits and add one then we'll have an unsigned integer greater than or equal to one and less than or equal to that largest integer.

Next we use the bitwise exclusive or operation to combine the array of taps into the mask value with the smallMask function given in listing 2.

Listing 2: smallMask
function smallMask(bits, taps) {
 var mask = 0;
 var n = taps.length;
 var i;

 for(i=0;i<n;++i) mask |= 1<<(taps[i]-1);
 mask |= 1<<(bits-1);
 return mask;
}

Note that the leftmost bit is always set with the final bitwise or.

Finally, the function that smallRnd returns extracts the rightmost bit with a bitwise and, shifts the bits to the right, taking care to shift the sign bit too by using the >>> operator, and then applies all of the taps at once by replacing the state with the result of the bitwise exclusive or of it and the mask.

A Big LFSR

If the number of bits in the register is greater than 32 then we'll need an array of numbers rather than just one, each holding up to 32 bits of the register. Now it may be the case that not every one of those numbers contains a tapped bit and so we can improve the performance of such LFSRs by only creating masks for those that do, as is done by bigMask in listing 3.

Listing 3: bigMask
function bigMask(bits, taps) {
 var words = ak.ceil(bits/32);
 var mask = [];
 var curr = {};
 var n = taps.length;
 var tap = taps[0]-1;
 var i, pos;

 curr.pos = ak.floor(tap/32);
 curr.bits = 1<<(tap%32);
 mask.push(curr);

 for(i=1;i<n;++i) {
  tap = taps[i]-1;
  pos = ak.floor(tap/32);

  if(pos === curr.pos) {
   curr.bits |= 1<<(tap%32);
  }
  else {
   curr = {};
   curr.pos = pos;
   curr.bits = 1<<(tap%32);
   mask.push(curr);
  }
 }

 if(curr.pos<words-1) {
  curr = {};
  curr.pos = words-1;
  curr.bits = 0;
  mask.push(curr);
 }
 curr.bits |= 1<<((bits-1)%32);

 return mask.length > 1 ? mask : mask[0].bits;
}

This function assumes that the taps are arranged in ascending order, that the first element of the register's array represents its rightmost 32 bits and that the last element represents its leftmost, in keeping with the rightmost to leftmost bits of the register being numbered in increasing order.
The latter assumption means that the \(n^\mathrm{th}\) bit of the register lies within the \(\lfloor n/32\rfloor^\mathrm{th}\) element of the array, where the odd-looking brackets stand for the greatest integer that is less than or equal to the number that they enclose, or in other words its floor. Since this is trivially equal to the result of the integer division of \(n\) by 32 we can determine which bit it is within that element with the remainder operation.
The former means that we can simply iterate through the taps updating the current element of the mask whenever a tap applies to the same element as it does and adding a new element to the mask whenever it doesn't, ensuring that we need only apply exclusive or operations to those elements that contain any taps. Consequently the elements of the mask record both the position within the register's array to which they should be applied and those bits of that element that they should tap in the pos and bits members respectively.

Note that the final if statement ensures that the element containing the leftmost bit is always included in the mask and the following exclusive or operation makes sure that it is always set to one.
This also means that if the only tapped bits are in the element containing the leftmost bits then we can simply apply the mask to that element alone, an efficiency improvement that we indicate is possible by returning the mask as a single number rather than as an array and which is exploited by the smallMaskRnd function given in listing 4.

Listing 4: smallMaskRnd
function smallMaskRnd(state, mask) {
 var words = state.length;

 return function() {
  var b = state[0]&0x01;
  var i;

  state[0] >>>= 1;
  for(i=1;i<words;++i) {
   state[i-1] |= (state[i]&0x01)<<31;
   state[i] >>>= 1;
  }

  if(b) state[words-1] ^= mask;
  return b;
 };
}

The function that it returns must still iterate through the elements of the register to shift their bits to the right, making sure that their rightmost bits are shifted into the leftmost bits of the elements that precede them and hence contain the bits that lie immediately to their right in the register.
After this we only need to apply an exclusive or operation to the last element if the output bit is one, as opposed to the general case implemented by bigMaskRnd in listing 5.

Listing 5: bigMaskRnd
function bigMaskRnd(state, mask) {
 var words = state.length;
 var masks = mask.length;

 return function() {
  var b = state[0]&0x01;
  var i;

  state[0] >>>= 1;
  for(i=1;i<words;++i) {
   state[i-1] |= (state[i]&0x01)<<31;
   state[i] >>>= 1;
  }

  if(b) for(i=0;i<masks;++i) state[mask[i].pos] ^= mask[i].bits;
  return b;
 };
}

The bigRnd function given in listing 6 chooses between using smallMaskRnd and bigMaskRnd to create the LFSR after seeding the elements of the register using ak.ranfRnd much as we did in smallRnd, albeit only ensuring that the first element, containing the rightmost bits, cannot equal zero.

Listing 6: bigRnd
var TWO_P32 = Math.pow(2, 32);

function bigRnd(bits, taps, seed) {
 var words = ak.ceil(bits/32);
 var rnd = ak.ranfRnd(seed);
 var state = new Array(words);
 var mask = bigMask(bits, taps);
 var i;

 state[0] = ak.floor((TWO_P32-1) * rnd())+1;
 for(i=1;i<words-1;++i) state[i] = ak.floor(TWO_P32 * rnd());
 state[words-1] = ak.floor(Math.pow(2, bits%32) * rnd());

 return ak.nativeType(mask)===ak.ARRAY_T ? bigMaskRnd(state, mask)
                                         : smallMaskRnd(state, mask);
}

ak.galoisRnd

As it happens there are two different, functionally equivalent, types of LFSR and those that we have been discussing are known as Galois LFSRs. We shall therefore call the function that glues all of this together ak.galoisRnd, as given in listing 7.

Listing 7: ak.galoisRnd
ak.galoisRnd = function(bits, taps, seed) {
 var n, i;

 if(ak.nativeType(bits)!==ak.NUMBER_T) {
  throw new Error('non-numeric bit count in ak.galoisRnd');
 }

 if(!isFinite(bits) || bits<2 || bits!==ak.floor(bits)) {
  throw new Error('invalid bit count in ak.galoisRnd');
 }

 if(ak.nativeType(taps)===ak.NUMBER_T) taps = [taps];
 if(ak.nativeType(taps)!==ak.ARRAY_T || taps.length===0) {
  throw new Error('invalid taps in ak.galoisRnd');
 }

 n = taps.length;
 taps = taps.slice(0);
 taps.sort(function(a,b){return a-b;});

 if(taps[0]<1 || taps[n-1]>bits) {
  throw new Error('invalid tap in ak.galoisRnd');
 }

 for(i=0;i<n;++i) {
  if(ak.nativeType(taps[i])!==ak.NUMBER_T) {
   throw new Error('non-numeric tap in ak.galoisRnd');
  }

  if(taps[i]!==ak.floor(taps[i])) {
   throw new Error('invalid tap in ak.galoisRnd');
  }

  if(i>0 && taps[i]===taps[i-1]) {
   throw new Error('repeated tap in ak.galoisRnd');
  }
 }

 return bits<=32 ? smallRnd(bits, taps, seed) : bigRnd(bits, taps, seed);
};

Most of this function is simply ensuring that the arguments are valid; firstly that the number of bits is a finite integer greater than one and secondly that the positions of the taps is either a single finite integer, or an array of them, that are greater than zero and less than or equal to the number of bits.
By replacing an integer tap with an array containing just that one tap we need only consider the array case when checking that the taps are valid. We do this by first creating a sorted copy of the array, which means that we only need to check that the first and the last of them are within the valid range. We then check that each of them is a number, then an integer and finally that they are unique. Note that we don't need to check that they are finite since in that case they will fail either the range test, in the case of infinities, or the integer test, in the case of NaNs.
Finally, we return an LFSR created by either smallRnd or bigRnd depending upon the number of bits in the register.

Program 1 demonstrates ak.galoisRnd by using it to implement the eight bit LFSR shown in figure 1.

Program 1: Drawing From The Eight Bit LFSR... Again

As we saw last time this LFSR has the maximum cycle length for an eight bit register of 255 bits, reflected by the fact that the pattern repeats after every five rows of 31 bits.

Since ak.galoisRnd can create LSFRs having many, many more than 32 bits we can create sequences with truly enormous cycle lengths if we choose the taps carefully, as demonstrated by program 2.

Program 2: The Practically Never Ending Cycle

As you will see if you run it several times, there is no discernible pattern in the bits that the generator outputs.

You can find the implementation of ak.galoisRnd in GaloisRnd.js.

ak.wardMoltenoRnd

The LFSR in program 2 is one of those that have cycles of maximal length for particular sizes of registers that Ward and Molteno went to the trouble of finding for us. Specifically, it cycles after \(2^{4096}-1\) bits have been drawn which in decimal is slightly larger than two followed by 1,233 zeros, and so it should come as no surprise that there weren't any discernible patterns!
It rather makes sense then to specifically implement some of their LFSRs with ak.galoisRnd generators, as is done in listing 8.

Listing 8: WardMoltenoRnd.js
ak.wardMolteno8Rnd = function(seed) {
 return ak.galoisRnd(8, [4, 5, 6], seed);
};

ak.wardMolteno16Rnd = function(seed) {
 return ak.galoisRnd(16, [11, 13, 14], seed);
};

ak.wardMolteno32Rnd = function(seed) {
 return ak.galoisRnd(32, [25, 26, 30], seed);
};

ak.wardMolteno64Rnd = function(seed) {
 return ak.galoisRnd(64, [60, 61, 63], seed);
};

ak.wardMolteno128Rnd = function(seed) {
 return ak.galoisRnd(128, [121, 126, 127], seed);
};

ak.wardMolteno256Rnd = function(seed) {
 return ak.galoisRnd(256, [246, 251, 254], seed);
};

ak.wardMolteno512Rnd = function(seed) {
 return ak.galoisRnd(512, [504, 507, 510], seed);
};

ak.wardMolteno1024Rnd = function(seed) {
 return ak.galoisRnd(1024, [1001, 1002, 1015], seed);
};

ak.wardMolteno2048Rnd = function(seed) {
 return ak.galoisRnd(2048, [2029, 2034, 2035], seed);
};

ak.wardMolteno4096Rnd = function(seed) {
 return ak.galoisRnd(4096, [4069, 4081, 4095], seed);
};

ak.wardMoltenoRnd = ak.wardMolteno4096Rnd;

These optimal generators can be found in WardMoltenoRnd.js.

As a final note, the idea of combining elements of a feedback buffer can be generalised to numbers as well as bits and it is to such generators that we shall turn our attention in the next post.

References

[1] Get A Shift On, www.thusspakeak.com, 2016.

[2] Pseudo Algorithmical, www.thusspakeak.com, 2016.

[3] Ward, R. & Molteno, T., Table of Linear Feedback Shift Registers, 2007.

[4] Talkin' 'Bout My Generators.html, www.thusspakeak.com, 2016.

Leave a comment