Systematically Random Numbers

The following article is from the Treasue Chest column in the Summer 1988 issue of APL*PLUS News. 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: ’ 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