Talkin' 'Bout My Generators

| 0 Comments

Last time[1] we took a look at pseudo random number generators, being entirely deterministic sequences of numbers that pass many of the tests that we might use to try to identify truly random sequences.
Specifically, we considered congruential generators which, starting from an initial value of \(x_0\), known as the seed, follow the rule
\[ x_{i+1} = \left(a \times x_i + c\right)\;\%\;m \]
where \(x\;\%\;m\) is the remainder of integer division of \(x\) by \(m\), \(a\) is known as the multiplier, \(c\) as the increment and \(m\) as the modulus.
We saw that we needed to be careful with our choices of the values of \(a\), \(c\), \(m\) and \(x_0\) if we wished to construct sequences that approximate randomness, not least that we should choose values for \(a\), \(c\) and \(x_0\) that are strictly less than \(m\), and the easiest way to do so is to simply copy those that have stood the test of time.

A Minimum Standard

In the late eighties, Park and Miller wrote a review of such generators[2] and suggested that a reasonable choice might be
\[ x_{i+1} = \left(16,807 \times x_i\right)\;\%\;2,147,483,647 \]
To implement this on a system with 32 bit integers requires some care, however, since the largest possible result of the product is
\[ 16,807 \times 2,147,483,646 = 36,092,757,621,515 > 4,294,967,295 = 2^{32}-1 \]
Derivation 1: Schrage's Algorithm
Firstly, we have
\[ x \;\%\; m = x - m \times \bigg\lfloor\frac{x}{m}\bigg\rfloor \]
so that
\[ \begin{align*} a \times (x \;\%\; q) &= a \times \left(x - q \times \bigg\lfloor\frac{x}{q}\bigg\rfloor\right)\\ &= a \times x - a \times q \times \bigg\lfloor\frac{a \times x}{a \times q}\bigg\rfloor\\ &= a \times x - (m-r) \times \bigg\lfloor\frac{a \times x}{m-r}\bigg\rfloor\\ &= a \times x - m \times \bigg\lfloor\frac{a \times x}{m-r}\bigg\rfloor + r \times \bigg\lfloor\frac{x}{q}\bigg\rfloor \end{align*} \]
Noting that
\[ \begin{align*} \frac{a \times x}{m-r} - \frac{a \times x}{m} & = \frac{a \times x \times m - a \times x \times (m-r)}{m \times (m-r)}\\ & = \frac{x \times a \times r}{m \times a \times q} = \frac{x}{m} \times \frac{r}{q} < 1 \end{align*} \]
it must be the case that
\[ \bigg\lfloor\frac{a \times x}{m-r}\bigg\rfloor - \bigg\lfloor\frac{a \times x}{m}\bigg\rfloor \]
equals either zero or one.

Consequently, either
\[ a \times (x \;\%\; q) - r \times \bigg\lfloor\frac{x}{q}\bigg\rfloor = (a \times x ) \;\%\; m \]
or
\[ a \times (x \;\%\; q) - r \times \bigg\lfloor\frac{x}{q}\bigg\rfloor = (a \times x ) \;\%\; m - m \]
Fortunately, Schrage discovered a method[3] for computing the remainder of a product upon integer division that only requires the addition and subtraction of terms having magnitudes that are no greater than the divisor by noting that
\[ m = a \times q + r \]
is satisfied by
\[ \begin{align*} q &= \bigg\lfloor\frac{m}{a}\bigg\rfloor\\ r &= m \;\%\; a \end{align*} \]
where the odd looking brackets stand for the largest integer less than or equal to the term that they surround, or in other words the floor, and if \(r\) is no greater than \(q\) we have
\[ \begin{align*} y &= a \times (x\;\%\; q)\\ z &= r \times \bigg\lfloor\frac{x}{q}\bigg\rfloor\\ (a \times x)\;\%\;m &= \begin{cases} y-z &\mathrm{if}\,z \leqslant y\\ m-(z-y) & \mathrm{otherwise} \end{cases} \end{align*} \]
as proven in derivation 1.

A Recursive Schrage's Algorithm

Even if \(r\) exceeds \(q\), we have
\[ \begin{align*} y &= a \times (x\;\%\; q)\\ &\leqslant a \times (q-1)\\ &= a \times q - a\\ &< a \times q + r\\ &= m \end{align*} \]
and so only \(z\) can exceed \(m\). To calculate it, note that if we define
\[ \begin{align*} a^\prime &= r\\ x^\prime &= \bigg\lfloor\frac{x}{q}\bigg\rfloor \end{align*} \]
then its remainder upon integer division by \(m\) is given by
\[ z \;\%\; m = \left(a^\prime \times x^\prime\right) \;\%\; m \]
and we can use Schrage's algorithm again!

This time, we have
\[ \begin{align*} q^\prime &= \bigg\lfloor\frac{m}{a^\prime}\bigg\rfloor\\ r^\prime &= m \;\%\; a^\prime \end{align*} \]
and since \(a^\prime\) equals \(r\), which must be strictly less than \(a\)
\[ \begin{align*} q^\prime &\geqslant q\\ r^\prime &< r \end{align*} \]
Now, if \(r^\prime\) is greater than \(q^\prime\) we can play the same trick again and, by doing so recursively, we'll get a sequence of non-decreasing \(q^\prime\) terms and strictly decreasing \(r^\prime\) terms and so must eventually get to a point where the latter is no greater than the former and we can terminate the recursion.
The schrageMul function given in listing 1 exploits this trick to build a function that calculates the remainder for \(x\) given fixed \(a\) and \(m\).

Listing 1: schrageMul
function schrageMul(a, m) {
 var q = ak.floor(m/a);
 var r = m%a;
 var f, g;

 if(r<=q) {
  f = function(x) {
   var lhs = a*(x%q);
   var rhs = r*ak.floor(x/q);
   return lhs>=rhs ? lhs-rhs : m-(rhs-lhs);
  };
 }
 else {
  g = schrageMul(r, m);
  f = function(x) {
   var lhs = a*(x%q);
   var rhs = g(ak.floor(x/q));
   return lhs>=rhs ? lhs-rhs : m-(rhs-lhs);
  };
 }
 return f;
}

The important thing to note about schrageMul is that the function that it builds doesn't itself test for the termination of the recursion, but rather that it is recursively constructed to behave appropriately whether or not the termination condition has been reached by creating a new recursive schrageMul function g to calculate \(z \;\%\; m\) if necessary.
The reason for this is that we shall typically call the resulting function many, many times for a particular choice of \(a\) and \(m\) and so it is much more efficient to perform the test whilst building the function than it would be to do so when calling it.

We can improve upon this by noting that if
\[ r \times \bigg\lfloor\frac{x}{q}\bigg\rfloor \]
is smaller than the largest integer that JavaScript can exactly represent then we can calculate \(z \;\%\; m\) without further recursion. The improved schrageMul function given in listing 2 does this by keeping track of the largest possible value that \(x^\prime\) can have at any point in the recursion in sup, starting with
\[ \bigg\lfloor\frac{m-1}{q}\bigg\rfloor \]
and successively setting it to the integer division of the previous maximum value by the current \(q^\prime\) term.

Listing 2: An Improved schrageMul
function schrageMul(a, m, sup) {
 var q = ak.floor(m/a);
 var r = m%a;
 var f, g;

 sup = isNaN(sup) ? ak.floor((m-1)/q) : ak.floor(sup/q);

 if(r<=q) {
  f = function(x) {
   var lhs = a*(x%q);
   var rhs = r*ak.floor(x/q);
   return lhs>=rhs ? lhs-rhs : m-(rhs-lhs);
  };
 }
 else if(r*sup<=ak.INT_MAX+1) {
  f = function(x) {
   var lhs = a*(x%q);
   var rhs = (r*ak.floor(x/q))%m;
   return lhs>=rhs ? lhs-rhs : m-(rhs-lhs);
  };
 }
 else {
  g = schrageMul(r, m, sup);
  f = function(x) {
   var lhs = a*(x%q);
   var rhs = g(ak.floor(x/q));
   return lhs>=rhs ? lhs-rhs : m-(rhs-lhs);
  };
 }
 return f;
}

You may be forgiven if you think that there's something fishy about using ak.INT_MAX+1 as the largest integer, rather than ak.INT_MAX. The reason for this is that JavaScript numbers are IEEE double precision floating point numbers, having 53 bits of precision[4]. When using them as integers it is therefore convenient to limit them to 53 bits plus the sign and ak.INT_MAX is the largest positive value that such integers can have.
However, all powers of two within the range of double precision numbers are also exactly representable and since ak.INT_MAX+1 is just such a number we might as well use it instead!
Note that, for the sake of simplicity when calling schrageMul, if sup is undefined then we treat it as if it were equal to the maximum possible value of \(m-1\), using isNaN to test for this case.

Multiplicative Congruential Generators

The simplest congruential generators are those in which \(c\) equals zero, known as multiplicative congruential generators.
If it's the case that \(a\) times the largest possible value of \(x\), which must be smaller than \(m\), is no greater than ak.INT_MAX+1 then we don't need Schrage's algorithm, as illustrated by smallMCG in listing 3.

Listing 3: A Small Multiplicative Congruential Generator
function isSmall(a, m) {
 return a*(m-1)<=ak.INT_MAX+1;
}

function smallMCG(a, m, state) {
 return function() {
  state.x = (a*state.x)%m;
  return (state.x-1)/(m-1);
 };
}

Derivation 2: The Positivity of x
If \(x\) can ever equal zero, then it must be the case that, for non-negative integers \(n\) and \(k\)

\(\displaystyle{ \begin{align*} a^n \times x_0 &= k \times m\\ \frac{a^n \times x_0}{m} &= k \end{align*} }\)

which is impossible if neither \(a\) nor \(x_0\) have any factors in common with \(m\), which we have already proven necessary if we want to generate a sequence that returns to the seed as slowly as possible.
Here, the state object contains the current value of \(x\) which is updated by the multiplicative congruence each time that the function that it returns is called.
Note that, by returning the updated value minus one divided by \(m\) minus one, that function approximates a standard uniform random variable, returning pseudo random numbers greater than or equal to zero and less than one, just like Math.random. This is because, if we choose the seed and the multiplier sensibly, \(x\) can never equal zero and so \(x\) minus one cannot be negative, as proven in derivation 2.
If it is not such a case, we simply use schrageMul to create a function that does the same thing, as illustrated by the bigMCG function given in listing 4.

Listing 4: A Big Multiplicative Congruential Generator
function bigMCG(a, m, state) {
 var mul = schrageMul(a, m);

 return function() {
  state.x = mul(state.x);
  return (state.x-1)/(m-1);
 };
}

Linear Congruential Generators

For non-zero \(c\) we must also be careful that the addition can't exceed the maximum integer. Noting again that \(x\) must be smaller than \(m\) then if \(a\) times \(m-1\) plus \(c\) doesn't exceed the largest integer we can perform the calculation directly with JavaScript numbers, as illustrated by tinyLCG in listing 5.

Listing 5: A Tiny Linear Congruential Generator
function isTiny(a, c, m) {
 return a*(m-1)+c<=ak.INT_MAX+1;
}

function tinyLCG(a, c, m, state) {
 return function() {
  state.x = (a*state.x+c)%m;
  return state.x/m;
 };
}

If the whole formula won't fit, but the product alone will, we can split the formula into two steps; first compute the remainder of the product and then the remainder of the sum. If it so happens that this sum can also exceed the maximum, then we can exploit the fact that \(c\) must also be less than \(m\) with
\[ -m < x-(m-c) < m \]
as is done by smallLCG in listing 6.

Listing 6: A Small Linear Congruential Generator
function smallLCG(a, c, m, state) {
 var f;

 if(m-1+c<=ak.INT_MAX+1) {
  f = function() {
   state.x = ((a*state.x)%m + c)%m;
   return state.x/m;
  };
 }
 else {
  f = function() {
   state.x = (a*state.x)%m - (m-c);
   if(state.x<0) state.x += m;
   return state.x/m;
  };
 }
 return f;
}

Finally, if the product can also be too big then we simply calculate it with Schrage's algorithm again, as shown by bigLCG in listing 7.

Listing 7: A Big Linear Congruential Generator
function bigLCG(a, c, m, state) {
 var mul = schrageMul(a, m);
 var f;

 if(m-1+c<=ak.INT_MAX+1) {
  f = function() {
   state.x = (mul(state.x) + c)%m;
   return state.x/m;
  };
 }
 else {
  f = function() {
   state.x = mul(state.x) - (m-c);
   if(state.x<0) state.x += m;
   return state.x/m;
  };
 }
 return f;
}

Seeding The Generators

Seeding linear congruential generators is slightly simpler than seeding multiplicative congruential generators since, provided \(c\) has been chosen appropriately, we won't run into trouble if the seed is zero or has a common factor with \(m\) and so we only need to check that it is a finite non-negative integer before assigning its remainder upon division by \(m\) to \(x\).
The seedLCG function given in listing 8 returns a function that does precisely this if it's passed an argument and uses Math.random() to randomly choose a valid seed otherwise.

Listing 8: Seeding A Linear Congruential Generator
function seedLCG(m, x0, state) {
 if(ak.nativeType(x0)===ak.UNDEFINED_T) {
  state.x = ak.floor(Math.random()*m);
 }
 else if(ak.nativeType(x0)!==ak.NUMBER_T) {
   throw new Error('non-numeric seed in ak.congruentialRnd');
 }
 else if(!isFinite(x0) || x0!==ak.floor(x0) || x0<0) {
   throw new Error('invalid seed in ak.congruentialRnd');
 }
 else state.x = x0%m;
 }
}

To seed a multiplicative congruential generator we must replace values of zero with one and increment values greater than one until their highest common factor with \(m\) is one. The function returned by the seedMCG function given in listing 9 uses the ak.hcf function, which we implemented to remove common factors from the numerators and denominators of our ak.rational type[5], to determine when to stop doing the latter after having used a function created by seedLCG to initialise state.x.

Listing 9: Seeding A Multiplicative Congruential Generator
function seedMCG(m, x0, state) {
 seedLCG(m, x0, state);
 while(ak.hcf(state.x, m)!==1) ++state.x;
}

ak.congruentialRnd

Armed with these helper functions, it is a relatively simple matter to implement a function that creates congruential generators that efficiently compute pseudo random numbers for any given choices of \(a\), \(c\) and \(m\), as is done by ak.congruentialRnd in listing 10.

Listing 10: ak.congruentialRnd
ak.congruentialRnd = function(a, m, c, x0) {
 var state = {};
 var rnd;

 if(ak.nativeType(a)!==ak.NUMBER_T) {
  throw new Error('non-numeric multiplier in ak.congruentialRnd');
 }
 a = Number(a);

 if(ak.nativeType(m)!==ak.NUMBER_T) {
  throw new Error('non-numeric modulus in ak.congruentialRnd');
 }
 m = Number(m);

 if(ak.nativeType(c)===ak.UNDEFINED_T) c = 0;
 else if(ak.nativeType(c)!==ak.NUMBER_T) {
  throw new Error('non-numeric offset in ak.congruentialRnd');
 }
 c = Number(c);

 if(m!==ak.floor(m) || m<3 || m>ak.INT_MAX+1) {
  throw new Error('invalid modulus in ak.congruentialRnd');
 }

 if(a!==ak.floor(a) || a<2 || a>=m) {
  throw new Error('invalid multiplier in ak.congruentialRnd');
 }

 if(c!==ak.floor(c) || c<0 || c>=m) {
  throw new Error('invalid offset in ak.congruentialRnd');
 }

 if(c>0) {
  if(isTiny(a, c, m)) rnd = tinyLCG(a, c, m, state);
  else if(isSmall(a, m)) rnd = smallLCG(a, c, m, state);
  else rnd = bigLCG(a, c, m, state);
  seedLCG(m, x0, state);
 }
 else {
  if(ak.hcf(a, m)!==1) throw new Error('invalid multiplier in ak.congruentialRnd');
  if(isSmall(a, m)) rnd = smallMCG(a, m, state);
  else rnd = bigMCG(a, m, state);
  seedMCG(m, x0, state);
 }
 return rnd;
};

This function, which can be found in CongruentialRnd.js, simply ensures that the multiplier \(a\), the increment \(c\) and the modulus \(m\) have valid values before using the helper functions to create an appropriate generator and randomly seeding it.
An important thing to note is that, unlike Math.random, it creates random number generators rather than being one.

Program 1 tests ak.congruentialRnd by comparing the sequences it generates with those generated using arbitrary precision arithmetic[6] with our ak.bignum type and its overloaded arithmetic operators.

Program 1: Testing ak.congruentialRnd

Note that we're using a seed of one since it satisfies the requirements for a good seed for any valid choices of \(a\), \(c\) and \(m\). Incidentally, the program's default choice for the latter are taken from Cray's RANF generator and requires that schrageMul is applied recursively. I'd encourage you to try a variety of choices for them to convince yourself that two sequences are always equal.

A Fistful Of Generators

Given just how difficult it is to pick good values for \(a\), \(c\) and \(m\), we shall typically use ak.congruentialRnd to create generators that others have devised. The ak library does so for

amcGeneratorSource
1,664,5254,294,967,2961,013,904,223ak.knuthRndKnuthRnd.js
16,8072,147,483,6470ak.minstd1RndMINSTDRnd.js
48,2712,147,483,6470ak.minstd2RndMINSTDRnd.js
65,5392,147,483,6480ak.randuRndRANDURnd.js
44,485,709,377,909281,474,976,710,6560ak.ranfRndRANFRnd.js

Note that ak.minstdRnd is also provided as a synonym for ak.minstd2Rnd.

We have taken some trouble to allow the user to supply their own random number generators whenever we have needed a source of randomness. For example, we can use ak.ranfRnd instead of defaulting to JavaScript's Math.random to generate normally distributed random numbers with our ak.normalRnd[7], as demonstrated by program 2 which plots a scaled line graph of the normal distribution's PDF over a bar chart of them.

Program 2: Replacing Math.random

Whilst congruential generators are extremely computationally efficient they're by no means the best source of pseudo random numbers, they must start repeating after no more than \(m\) numbers have been generated for example, and so in the next post we shall take a look at how we might improve them.

References

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

[2] Park, S. & Miller, K., Random Number Generators: Good Ones are Hard to Find, Communications of the ACM, Vol. 31, No. 10, 1988.

[3] Schrage, L., A more portable FORTRAN random number generator, ACM Transactions on Mathematical Software, Vol. 5, No. 2, 1979.

[4] You're Going To Have To Think! Why [Insert Technique Here] Won't Cure Your Floating Point Blues , www.thusspakeak.com, 2013.

[5] You're Going To Have To Think! Why Rational Arithmetic Won't Cure Your Floating Point Blues, www.thusspakeak.com, 2013.

[6] You're Going To Have To Think! Why Arbitrary Precision Arithmetic Won't Cure Your Floating Point Blues, www.thusspakeak.com, 2013.

[7] Define Normal, www.thusspakeak.com, 2014.

Leave a comment