Let's Twist Again

| 0 Comments

In the last few posts[1][2] we have been looking at linear feedback shift registers, or LFSRs, which generate random bits from a buffer by taking the rightmost bit as the output, shifting the buffer one bit to the right and, if the output is one, applying an exclusive or operation with a set of bits known as taps. Since the leftmost bit is always one of the taps, they effectively shift the output bit back into the leftmost bit of the buffer.
We saw that despite looking completely different to the congruential generators that we covered earlier[3], they are surprisingly similar, not least in that we must be extremely careful when setting them up if we want to generate sequences of bits that take as long as possible to enter into a cycle.
Finally, we noted that it was possible to generalise LFSRs so that they generate numbers rather than bits, which we shall do in this post.

Integer LFSRs

Figure 1 shows a graphical representation of an eight bit LFSR with the filled circles representing ones, the empty circles representing zeros, the \(\oplus\) symbols representing exclusive or operations and the arrows indicating the flow of bits through the register at each step.

Figure 1: A Bitwise Linear Feedback Shift Register

If we note that the exclusive or operation of ones and zeros is equivalent to the remainder of their addition after integer division by two
\[ x \oplus y = (x+y)\;\%\;2 \]
where \(\%\) is the remainder operation, then we have a hint as to how we might generalise LFSRs to positive integers; fill the buffer with positive integers instead of bits and use the taps to calculate the remainders of sums after division by an integer greater than two, known as the base of the LFSR.
In bitwise LFSRs we can consider the presence of a tap as multiplying the output bit by one and the absence of a tap as multiplying it by zero since in the latter case the subsequent exclusive or operation wouldn't have any effect and so we might as well remove the tap. In generalising them to integers, each tap will therefore multiply the output by some constant before being added to a number in the register.

Program 1 illustrates the idea with a buffer of three integers with the taps multiplying the output by zero, one and five, adding the results to the first, second and third elements of the buffer respectively and calculating the remainder of the sums after integer division by seven.

Program 1: An Integer LFSR

Just like bitwise LFSRs, the sequences generated by integer LFSRs will be the maximum possible length if the buffer cycles through every possible arrangement of three integers greater than or equal to zero and less than its base, with the exception of all zeros. In program 1 the base is seven and so we cannot do better than
\[ 7^3 - 1 = 343 - 1 = 342 \]
and program 2 demonstrates that this is indeed the length of the generated sequence.

Program 2: The Length Of The Sequence

If you try out a few different values for the taps' multipliers you will find that it is no easier to set up integer LFSRs so that they generate sequences of maximal length than it was to set up congruential generators so that they do so. Which is a pity, really.

Note that integer LFSRs that have buffers of just one element are identical to multiplicative congruential generators since that element is set to zero when the output is shifted out, then the product of the tap and the output is added to it before applying the remainder operation which is trivially equivalent to
\[ x_{i+1} = (a \times x_i)\;\%\;m \]
for a base, or modulus in the language of congruential generators, of \(m\) and a tap having a multiplier of \(a\).

We can generalise integer feedback registers still further by using arbitrary operators to manipulate the numbers in their buffers; in particular we might prefer to use bitwise operators since they are faster than arithmetic operators.
Before take a look at one that does so, however, it's worth reminding ourselves about how JavaScript's bitwise operators work.

JavaScript's Bitwise Operators

Since bitwise operators are typically applicable to integers rather than floating point numbers, it isn't particularly surprising that JavaScript converts its floating point numbers to integers before applying them. What is rather surprising, however, is that it converts them to signed integers rather than unsigned integers.
I rather glossed over the details when implementing bitwise LFSRs, which I could get away with since we only used bitwise operations which, provided we remember to use the >>> operator to perform right shifts so that the sign bit is also shifted, behave just as they would with unsigned integers.
Things get a little trickier when we combine both bitwise and arithmetic operations, and before doing so we'll need to understand exactly how JavaScript converts its numbers to integers.

Firstly, all digits to the right of the point are set to zero, effectively rounding the number towards zero. For a number x we might do this with

x = x<0 ? -Math.floor(-x) : Math.floor(x);

Note that NaN, plus or minus zeros and infinities are all converted to zero with a positive sign bit.
Next, x is set to its positive remainder after division by \(2^{32}\)

x %= Math.pow(2, 32);
if(x<0) x += Math.pow(2, 32);

Finally, if x is greater than or equal to \(2^{31}\) then \(2^{32}\) is subtracted from it

if(x>=Math.pow(2, 31)) x -= Math.pow(2, 32);

These rules mean that bitwise operations on positive numbers can have negative results. Fortunately they also mean that we can recover the unsigned integer that we might have expected by simply adding \(2^{32}\) to them.

The Mersenne Twister

The Mersenne Twister[4], arguably the most popular pseudo random number generator at present, uses bitwise operators to update the 624 elements in its buffer after shifting them, albeit in a single sweep after all of them have been drawn as outputs rather than after each of them has been.
This update is implemented by the fillMT function given in listing 1.

Listing 1: fillMT
var N = 624;
var M = 397;

var MA = 0x9908b0df;

function fillMT(mt) {
 var k, y;

 for(k=0;k<N-M;++k) {
  y = (mt[k+1]&0x7fffffff) | (mt[k]&0x80000000);
  mt[k] = y&0x01 ? mt[k+M] ^ (y>>>1) ^ MA : mt[k+M] ^ (y>>>1);
 }
 for(;k<N-1;++k) {
  y = (mt[k+1]&0x7fffffff) | (mt[k]&0x80000000);
  mt[k] = y&0x01 ? mt[k+M-N] ^ (y>>>1) ^ MA : mt[k+M-N] ^ (y>>>1);
 }
 y = (mt[0]&0x7fffffff) | (mt[N-1]&0x80000000);
 mt[N-1] = y&0x01 ? mt[M-1] ^ (y>>>1) ^ MA : mt[M-1] ^ (y>>>1);
}

Whilst I typically try to give proofs of why an algorithm is effective, in this case I must admit defeat; I have no idea why this is optimal and have simply translated the reference implementation from C to JavaScript.
I'm willing to trust its inventors' claim that the buffer cycles through \(2^{19,937}-1\) of its \(2^{19,968}-1\) possible non-zero states. Note that prime numbers, being those positive integers that cannot be wholly divided by any but themselves and one, of the form \(2^p-1\) are known as a Mersenne primes and this algorithm gets its name from the fact that its cycle length is just such a number.

The ak.mtRnd function given in listing 2 returns a function that uses fillMT to fill the buffer after every element has been drawn as an output.

Listing 2: ak.mtRnd
var TWO_M32 = Math.pow(2, -32);

ak.mtRnd = function(s) {
 var mt = new Array(N);
 var i = N;
 var y;

 seedMT[ak.nativeType(s)](s, mt);

 return function() {
  if(i===N) {
   fillMT(mt);
   i = 0;
  }

  y = mt[i++];
  y ^= y>>>11;
  y ^= (y<<7) & 0x9d2c5680;
  y ^= (y<<15) & 0xefc60000;
  y ^= y>>>18;

  return y<0 ? y*TWO_M32+1 : y*TWO_M32;
 };
};

Note that each output is further manipulated with a series of bitwise operations, known as tempering, before being mapped to a number greater than or equal to zero and strictly less than one by dividing it by \(2^{32}\) and adding one if it's negative to correct for the fact that buffer's elements are signed, rather than unsigned, integers.

And no, I don't know why the tempering is optimal either.

Seeding The Mersenne Twister

The reference implementation of the Mersenne Twister supports seeding with both integers and arrays of integers and so our implementation shall do the same. We shall use our usual scheme of dispatching the call to seedMT to methods of an object depending on the type of the argument, as illustrated by listing 3.

Listing 3: seedMT
var TWO_P16 = Math.pow(2, 16);

var A  = 1812433253;
var AL = A&0xffff;
var AH = A>>>16;

var seedMT = {};

seedMT[ak.NUMBER_T] = function(s, mt) {
 var i, x, xl, xh;

 if(!isFinite(s) || s<0 || s!=ak.floor(s)) {
  throw new Error('invalid seed in ak.mtRnd');
 }

 mt[0]= s%TWO_P32;
 for (i=1;i<N;++i) {
  x  = mt[i-1] ^ (mt[i-1]>>>30);
  xl = x&0xffff;
  xh = x>>>16;

  mt[i] = (AH*xl+AL*xh)*TWO_P16 + AL*xl + i;
 }
};

seedMT[ak.UNDEFINED_T] = function(s, mt) {
 seedMT[ak.NUMBER_T](ak.floor(Math.random()*TWO_P32), mt);
};

If the seed is undefined then we randomly generate one with Math.random and use the seedMT function for numbers to fill the buffer with the results of something like a congruential generator.

The first difference between the latter and the reference implementation are that, since JavaScript numbers are floating point, we must explicitly check that the seed is a finite non-negative integer. Next, the product of its multiplier A and the largest 32 bit integer has more bits than can be exactly represented with a floating point number and so we must take care not to lose the least significant bits to rounding.
Specifically, we must split both A and x into low and high parts and use long multiplication to compute the product. Note that all of these values will be positive since the low parts don't include the sign bit and the high parts have had it shifted to the right.
Since the generator computes the remainder after integer division by a modulus of \(2^{32}\) we can drop the product of the two high parts, which must be a multiple of it, so that the results of the congruence will fit within the 53 bits of precision of a floating point number. Finally, we don't actually need to calculate the remainder after division by the modulus since the bitwise operations of the Mersenne Twister will convert them to 32 bit integers and so will do it for us.

And no, I haven't a clue as to why this algorithm is optimal. Would you please stop asking awkward questions!

The function that seeds the Mersenne Twister with an array is given in listing 4.

Listing 4: Seeding With An Array
var TWO_P32 = Math.pow(2, 32);

var B  = 1566083941;
var BL = B&0xffff;
var BH = B>>>16;

var C = 1664525;

var seedMT = {};

seedMT[ak.ARRAY_T] = function(s, mt) {
 var n = s.length;
 var i = 1;
 var j = 0;
 var k, x, xl, xh;

 if(n===0) throw new Error('empty seed in ak.mtRnd');

 for(k=0;k<n;++k) {
  if(ak.nativeType(s[k])!==ak.NUMBER_T) {
   throw new Error('non-numeric seed in ak.mtRnd');
  }

  if(!isFinite(s[k]) || s[k]<0 || s[k]!=ak.floor(s[k])) {
   throw new Error('invalid seed in ak.mtRnd');
  }
  s[k] %= TWO_P32;
 }

 seedMT[ak.NUMBER_T](19650218, mt);

 for(k=Math.max(n, N);k>0;--k) {
  x = mt[i-1] ^ (mt[i-1]>>>30);
  mt[i] = (x>=0 ? (mt[i] ^ (x*C)) : (mt[i] ^ ((x+TWO_P32)*C))) + s[j] + j;

  if(++i===N) {mt[0] = mt[N-1]; i=1;}
  if(++j===n) j=0;
 }

 for(k=N-1;k>0;--k) {
  x  = mt[i-1] ^ (mt[i-1]>>>30);
  xl = x&0xffff;
  xh = x>>>16;

  mt[i] = (mt[i] ^ ((BH*xl+BL*xh)*TWO_P16 + BL*xl)) + TWO_P32-i;

  if(++i===N) {mt[0] = mt[N-1]; i=1;}
 }
 mt[0] = 0x80000000;
};

Note that since this uses both bitwise and arithmetic operations, we must be careful to add \(2^{32}\) to the results of the former to convert them to unsigned integers before we apply the latter, relying upon the fact that we have 53 bits to work with to be sure that we can exactly represent the results of arithmetic expressions involving them whether they're positive or negative.
Furthermore, just as we did when seeding with numbers, we must check that the elements of the seed are finite non-negative integers and use long multiplication if JavaScript's numbers don't have enough bits to store the results, as is the case in the second loop in which we can again ignore the products of the high terms since they would have been lost if we were using unsigned 32 bit integers. We can get away with simply subtracting i there since we've ensured that the result of the exclusive or operation is a positive number by adding \(2^{32}\) to it and so if the result of the subtraction is negative, it will be small enough that subsequent bitwise operations will behave as if we were using unsigned integers.
Finally, the first element of the buffer is set to a non-zero value to guarantee that they don't all equal zero.

And again, no, I don't know why this particular scheme was chosen; enough already with the blasted questions!

Using ak.mtRnd

Now I am painfully aware of how little justification that I have given for the effectiveness of the Mersenne Twister, but I can at least give you a demonstration of it in action with program 3.

Program 3: Using ak.mtRnd

Hopefully, this will be evidence enough for you that the Mersenne Twister is a pretty good approximation of a random variable, because it's all I have to offer!

Finally, the implementation of ak.mtRnd can be found in MTRnd.js.

References

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

[2] Shifting Up A Gear, www.thusspakeak.com, 2016.

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

[4] Matsumoto, M. & Nishimura, T., Mersenne Twister: A 623-dimensionally equidistributed uniform pseudorandom number generator, ACM Transactions on Modeling and Computer Simulation, Vol. 8, No. 1, 1998.

Leave a comment