Systematically Random Numbers

by Ralph L. Fox

STSC is often asked how APL's random number generator works. The following function emulates the primitive "roll" function in APL*PLUS (and many other APL) Systems.

    ’ R„OneRoll N;A;P
[1]   © Emulate ?N for scalar integer N. Global ‘RL
[2]   ©   is like ŒRL, the random link or "seed".
[3]
[4]   P„2147483647   © P„¯1+2*31 is a prime number
[5]                  © convenient because it is also the largest
[6]                  © integer on 32-bit machines.
[7]
[8]   A„16807        © A„7*5 is a "primitive root" of P
[9]                  © since its first P-1 powers (modulo P)
[10]                 © generate all integers from 1 to P-1.
[11]
[12]  ‘RL„P|‘RL×A    © First generate the next random link:
[13]
[14]  R„ŒIO+˜NבRL÷P © Then base the result on ‘RL's ratio to P:
    ’

g←⎕NEW Solitaire ⋄ g.Visible←0


Œ„¼…½’‘œ¯–•—

Funny font! Œ„¼…½’‘œ¯–•—

Note that OneRoll and ? gave the same value, and ‘RL is keeping pace with ŒRL. Run them again:

? 6789 ª OneRoll 6789
5130
5130
      ‘RL ª ŒRL
1622650073
1622650073

For more details on this approach, refer to the excellent article by Eugene McDonnell on the historical development of APL's random number generator in APL Quote-Quad, Volume 8, Number 3, March, 1978.

The obvious advantage of ? over OneRoll is that the primitive extends to higher dimensional arrays, whereas multi-element processing has to be programmed into a general emulator function. Here is a first attempt to program this emulator for a vector, without looping:

    ’ R„VRoll N;A;P;RLS
[1]   © (Poorly) emulate ?N for a vector N.
[2]
[3]   P„2147483647  © P„¯1+2*31
[4]   A„16807       © A„7*5
[5]   ‘RL„¯1†RLS„P|‘RL×A*(~ŒIO)+¼½N„,N
[6]   R„ŒIO+˜N×RLS÷P
    ’

This Function doesn't perform satisfactorily because raising A to large powers quickly produces round-off errors and eventually a LIMIT ERROR. For example,

      ‘RL„ŒRL„7*5
      VRoll 8½100
14 76 46 54 1 1 1 1
      ‘RL
0
      ? 8½100
14 76 46 54 22 5 68 68
      ŒRL
1458777923

Note that the VRoll result starts out correctly but quickly degenerates. On systems that support nested arrays and user-defined functions with operators, an elegant solution to the above quandary is available. The goal is to calculate the residue of pairwise-products sequentially, rather than obtain all residues at the end. To accomplish this, define an auxiliary subroutine just to perform the residue of the product:

    ’ R„A ResProd B
[1]   R„2147483647|A×B © R(¯1+2*31)|A×B
    ’

Use this function with the scan operator to produce:

    ’ R„VecRoll N;A;P;RLS
[1]   © Emulate ?N for a vector N. A subfn ResProd
[2]   ©   and a nested array system is required.
[3]
[4]   P„2147483647 © P„¯12*31
[5]   A„16807 © A„7*5
[6]   ‘RL„¯1†RLS„ResProd\(‘RL×A),(¯1+½N„,N)½A
[7]   R„ŒIO+˜N×RLS÷P
    ’

We can now verify that ? is properly emulated:

‘RL„7*5 ª VecRoll 8½100
14 76 46 54 22 5 68 68
      ‘RL
1458777923

Now that we know how the roll primitive is implemented, the next question that comes to mind is: how truly random can we consider its results? The answer depends on the type of simulation you have in mind. APL's roll function is a standard multiplicative congruential pseudo-random number generator that is perfectly adequate for general-purpose use. For specialize use, you should first apply relevant statistical tests before relying on its outpu. (For example, using roll to generate points in higher-dimensional space has been shown to have some geometrical deficiencies.) For more infformation on the randomness of roll, see the article by Thomas Herzog in APL Quote-Quad, Volume 12, Number 2, December, 1981. There he points out that the following algorithm offers a significantly improved random number generator at relatively modest performance cost:

    ’ R„ROLL N;K
[1]   © Variant of ?N for vector N; "randomness" is
[2]   ©   improved by adding a random shuffling.
[3]
[4]    K„½N„,N
[5]    R„?N
[6]    R„R[K?K]
    ’

If you would like to explore further the concept of highly-random number generators, statistical software can provide various goodness-of-fit, auto-correlation, and nonparametric tests for the uniformity of "randomly-generated" numbers. You can compare data samples generated by ? and the ROLL function. You might even like to scale down (the illustrated or a looping) VecRoll to use the prime P„31 and a primitive root A„11, and then generate and analyze the corresponding small univers (before repetition) of random numbers with

‘RL„11 ª DATA„VecRoll 30½30 ’ R„OneRoll N;A;P [1] © Emulate ?N for scalar integer N. Global ‘RL [2] © is like ŒRL, the random link or "seed". [3] [4] P„2147483647 © P„¯1+2*31 is a prime number [5] © convenient because it is also the largest [6] © integer on 32-bit machines. [7] [8] A„16807 © A„7*5 is a "primitive root" of P [9] © since its first P-1 powers (modulo P) [10] © generate all integers from 1 to P-1. [11] [12] ‘RL„P|‘RL×A © First generate the next random link: [13] [14] R„ŒIO+˜NבRL÷P © Then base the result on ‘RL's ratio to P: ’ To illustrate this emulator, we set ‘RL and ŒRL to the default value of ŒRL in a clear workspace, and compare executions: ‘RL„ŒRL„7*5 OneRoll 12345 1624 ?12345 1624 ‘RL 282475249 ŒRL 282475249 Note that OneRoll and ? gave the same value, and ‘RL is keeping pace with ŒRL. Run them again: ? 6789 ª OneRoll 6789 5130 5130 ‘RL ª ŒRL 1622650073 1622650073 For more details on this approach, refer to the excellent article by Eugene McDonnell on the historical development of APL's random number generator in APL Quote-Quad, Volume 8, Number 3, March, 1978. The obvious advantage of ? over OneRoll is that the primitive extends to higher dimensional arrays, whereas multi-element processing has to be programmed into a general emulator function. Here is a first attempt to program this emulator for a vector, without looping: ’ R„VRoll N;A;P;RLS [1] © (Poorly) emulate ?N for a vector N. [2] [3] P„2147483647 © P„¯1+2*31 [4] A„16807 © A„7*5 [5] ‘RL„¯1†RLS„P|‘RL×A*(~ŒIO)+¼½N„,N [6] R„ŒIO+˜N×RLS÷P ’ This Function doesn't perform satisfactorily because raising A to large powers quickly produces round-off errors and eventually a LIMIT ERROR. For example, ‘RL„ŒRL„7*5 VRoll 8½100 14 76 46 54 1 1 1 1 ‘RL 0 ? 8½100 14 76 46 54 22 5 68 68 ŒRL 1458777923 Note that the VRoll result starts out correctly but quickly degenerates. On systems that support nested arrays and user-defined functions with operators, an elegant solution to the above quandary is available. The goal is to calculate the residue of pairwise-products sequentially, rather than obtain all residues at the end. To accomplish this, define an auxiliary subroutine just to perform the residue of the product: ’ R„A ResProd B [1] R„2147483647|A×B © R(¯1+2*31)|A×B ’ Use this function with the scan operator to produce: ’ R„VecRoll N;A;P;RLS [1] © Emulate ?N for a vector N. A subfn ResProd [2] © and a nested array system is required. [3] [4] P„2147483647 © P„¯12*31 [5] A„16807 © A„7*5 [6] ‘RL„¯1†RLS„ResProd\(‘RL×A),(¯1+½N„,N)½A [7] R„ŒIO+˜N×RLS÷P ’ We can now verify that ? is properly emulated: ‘RL„7*5 ª VecRoll 8½100 14 76 46 54 22 5 68 68 ‘RL 1458777923 Now that we know how the roll primitive is implemented, the next question that comes to mind is: how truly random can we consider its results? The answer depends on the type of simulation you have in mind. APL's roll function is a standard multiplicative congruential pseudo-random number generator that is perfectly adequate for general-purpose use. For specialize use, you should first apply relevant statistical tests before relying on its outpu. (For example, using roll to generate points in higher-dimensional space has been shown to have some geometrical deficiencies.) For more infformation on the randomness of roll, see the article by Thomas Herzog in APL Quote-Quad, Volume 12, Number 2, December, 1981. There he points out that the following algorithm offers a significantly improved random number generator at relatively modest performance cost: ’ R„ROLL N;K [1] © Variant of ?N for vector N; "randomness" is [2] © improved by adding a random shuffling. [3] [4] K„½N„,N [5] R„?N [6] R„R[K?K] ’ If you would like to explore further the concept of highly-random number generators, statistical software can provide various goodness-of-fit, auto-correlation, and nonparametric tests for the uniformity of "randomly-generated" numbers. You can compare data samples generated by ? and the ROLL function. You might even like to scale down (the illustrated or a looping) VecRoll to use the prime P„31 and a primitive root A„11, and then generate and analyze the corresponding small univers (before repetition) of random numbers with ‘RL„11 ª DATA„VecRoll 30½30

Back to the articles index.