Systematically Random Numbers
by Ralph L. Fox
this is an example of code:
g←⎕NEW Solitaire ⋄ g.Visible←0Funny 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 1622650073For 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 1458777923Note 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 1458777923Now 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 illiustrate 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