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
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.
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
This time, we have
The
The important thing to note about
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
You may be forgiven if you think that there's something fishy about using as the largest integer, rather than
However, all powers of two within the range of double precision numbers are also exactly representable and since is just such a number we might as well use it instead!
Note that, for the sake of simplicity when calling
If it's the case that \(a\) times the largest possible value of \(x\), which must be smaller than \(m\), is no greater than then we don't need Schrage's algorithm, as illustrated by
Here, the
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
If it is not such a case, we simply use
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
Finally, if the product can also be too big then we simply calculate it with Schrage's algorithm again, as shown by
The
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
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
Program 1 tests
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
Note that
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
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.
[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.
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
Consequently, either
\[
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  (mr) \times \bigg\lfloor\frac{a \times x}{mr}\bigg\rfloor\\
&= a \times x  m \times \bigg\lfloor\frac{a \times x}{mr}\bigg\rfloor + r \times \bigg\lfloor\frac{x}{q}\bigg\rfloor
\end{align*}
\]
Noting that
\[
\begin{align*}
\frac{a \times x}{mr}  \frac{a \times x}{m}
& = \frac{a \times x \times m  a \times x \times (mr)}{m \times (mr)}\\
& = \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}{mr}\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
\]
\[
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}
yz &\mathrm{if}\,z \leqslant y\\
m(zy) & \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 (q1)\\
&= 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 nondecreasing \(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 ? lhsrhs : m(rhslhs); }; } else { g = schrageMul(r, m); f = function(x) { var lhs = a*(x%q); var rhs = g(ak.floor(x/q)); return lhs>=rhs ? lhsrhs : m(rhslhs); }; } 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{m1}{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((m1)/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 ? lhsrhs : m(rhslhs); }; } 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 ? lhsrhs : m(rhslhs); }; } else { g = schrageMul(r, m, sup); f = function(x) { var lhs = a*(x%q); var rhs = g(ak.floor(x/q)); return lhs>=rhs ? lhsrhs : m(rhslhs); }; } return f; }
You may be forgiven if you think that there's something fishy about using
ak.INT_MAX+1
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
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 \(m1\), 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
smallMCG
in listing 3.Listing 3: A Small Multiplicative Congruential Generator
function isSmall(a, m) { return a*(m1)<=ak.INT_MAX+1; } function smallMCG(a, m, state) { return function() { state.x = (a*state.x)%m; return (state.x1)/(m1); }; }
Derivation 2: The Positivity of x
If \(x\) can ever equal zero, then it must be the case that, for nonnegative integers \(n\) and \(k\)
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.
\(\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.
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.x1)/(m1); }; }
Linear Congruential Generators
For nonzero \(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 \(m1\) plus \(c\) doesn't exceed the largest integer we can perform the calculation directly with JavaScript numbers, as illustrated bytinyLCG
in listing 5.
Listing 5: A Tiny Linear Congruential Generator
function isTiny(a, c, m) { return a*(m1)+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(mc) < 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(m1+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  (mc); 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(m1+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)  (mc); 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 nonnegative 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('nonnumeric 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 byak.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('nonnumeric multiplier in ak.congruentialRnd'); } a = Number(a); if(ak.nativeType(m)!==ak.NUMBER_T) { throw new Error('nonnumeric 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('nonnumeric 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 useak.congruentialRnd
to create generators that others have devised. The ak
library does so fora  m  c  Generator  Source 

1,664,525  4,294,967,296  1,013,904,223  ak.knuthRnd  KnuthRnd.js 
16,807  2,147,483,647  0  ak.minstd1Rnd  MINSTDRnd.js 
48,271  2,147,483,647  0  ak.minstd2Rnd  MINSTDRnd.js 
65,539  2,147,483,648  0  ak.randuRnd  RANDURnd.js 
44,485,709,377,909  281,474,976,710,656  0  ak.ranfRnd  RANFRnd.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