Doc/Articles/Play083
How the Roll Function Works (in APL\360 and its Descendants)
Eugene McDonnell Eemcd AT aol DOT com
If you look up roll in the J Introduction and Dictionary [1], you read
? y yields a uniform random selection from the population i. y The random seed used begins at 7^5 (16807).
If you had access to the inside of the computer, you would see that the following occurs when the roll function is executed:
seed =: 2147483647 | 16807 seed (1)
Following this, ? y is determined as <. y * seed % 2147483647
This is about all you need to know to reproduce the results that the interpreter gives for the roll function, but knowing all this you still might not understand the roll function. In particular, you wouldn't know the significance of the two constants 2147483647 and 16807. To explain them (which is the purpose of this article) requires that a little groundwork be laid.
The number 2,147,483,647 (which assembly language programmers know as 7FFFFFFF) is famous in number theory. It was the largest prime known for over a hundred years (1772-1876). It is in fact (2^31)-1, and one of the Mersenne primes, which are of the form (2^P)-1, for P a prime. Not all primes P give a prime (2^P)-1, but that is another story. In 1644, 2,147,483,647 had been alleged to be prime by Father Marin Mersenne (one of the most famous allegators in all of mathematics), but it was not proven to be prime until Leonhard Euler (65 years old and completely blind) did so in 1772.
Euler's prime is interesting because, quite conveniently for 32-bit word computers, the largest integer that can be represented in 31 bits and a sign bit is also a prime number, namely Euler's prime.
The number 16807 is not as celebrated as Euler's prime. It is interesting to us because 2147483647|16807^2147483646 is equal to 1, and because there is no smaller exponent than the one shown that makes the expression equal to 1. In order to make this meaningful to you, study the results of the following two expressions:
8|^/~>:i.7 1 1 1 1 1 1 1 2 4 0 0 0 0 0 3 1 3 1 3 1 3 4 0 0 0 0 0 0 5 1 5 1 5 1 5 6 4 0 0 0 0 0 7 1 7 1 7 1 7 7|^/~>:i.6 1 1 1 1 1 1 2 4 1 2 4 1 3 2 6 4 5 1 4 2 1 4 2 1 5 4 6 2 3 1 6 1 6 1 6 1
Looking at the rows of these two results, we see that for the modulus 8, all of the rows contain repeated elements. For the modulus 7, on the other hand, two of the rows contain elements that are all distinct. A 1 appears as the last element of each row, but only for rows 3 and 5 is there only a single 1. It is generally true that when the modulus is a prime, some rows will have all distinct elements. You'll never guess how many such rows there are, in general. For a prime P there are always as many rows of this kind as there are numbers less than P-1 which have 1 as greatest common divisor (gcd) with P-1 (that is, which are relatively prime to P-1). As we can see above, for the prime 7, two rows have distinct values. Trying the rule given above (and remember J uses +. to signify gcd), 7-1 gives 6, and 6 +. 1 2 3 4 5 is 1 2 3 2 1; since there are two 1's, there should be two rows that have distinct values - - and there are. The elements of i. P-1 that produce all P-1 distinct values when their first P-1 powers are taken mod P are called primitive roots of P. So 3 and 5 are primitive roots of 7.
Let's put it another way, one that relates it to our purpose. If we start with 3 or 5 and take successive powers mod 7, then we will get all the values between 1 and 6, in some order, ending with 1 and with no repetitions. When we reach 1, the cycle will start again.
Perhaps you see now the relevance of primes and their primitive roots to the problem of random number generation. In the case of the roll function, we start with 16807 and take successive powers mod 2,147,483,647, generating all the values between 1 and 2,147,483,646, in some order, ending with 1 and with no repetitions. When we reach 1, the cycle starts again. The useful property of the permutation so generated is that with the proper primitive root, it satisfies many of the more important tests for a random series. When used in conjunction with (1) very satisfactory random numbers can be generated. The man who suggested this technique is D. H. Lehmer, the eminent second- generation number theorist from Berkeley, California. He made the suggestion in 1948. It had been described several times in print before APL\360 was implemented, first by B. W. Holz and C. E. Clark [2], and then again by D. W. Hutchinson [3]. Both of these articles referred to a 36-bit word machine, and both used the prime (2*35)-31 and the primitive roots 5*5 and 5*13. Both credited Lehmer for suggesting the method and the specific values to be used for the prime and the primitive root.
History of the Roll Function
The roll function was introduced to describe the indeterminacy of certain operations of the IBM System/360 computers which the computer architects decided not to define. This allowed the designers of the individual members of the System/360 family to have some freedom in designing their machines. The formal definition of System/360 written in APL [4] employed the query symbol (?) to describe this freedom.
When the symbol was introduced, its definition was different from what it is now. It came in three forms: niladic, monadic, and dyadic. In the niladic use, ? meant that the result was either 0 or 1, and one couldn't tell which. The monadic form ?(B) yielded a logical vector of length B containing an arbitrary pattern of 1's and 0's. In J this would be ?B#2. The dyadic form ?A(B) gave a logical vector of length B containing an arbitrary pattern of A 1's and B-A 0's. In terms of J's deal dyad, this would be written (B?B) { B {. A#1.
As you can see, the random functions were originally specific to the Boolean domain, and strongly oriented to machine description.
PROBLEM: Write a function that prepares a table of all patterns of A 1's and B-A 0's.
By the time the 7090 time sharing version of APL came into being in late 1965, the definitions of the random functions had changed almost to what they are today. However, only the roll function was implemented on the 7090. The algorithm implemented, which was one that had been used at the IBM Research Center in Yorktown Heights, New York, was crude in the extreme. A "random" number was generated by performing the exclusive-or of 20 sequential words in memory, beginning with words 0-19, and progressing through memory with each new random number request. Between the time of this implementation and the APL\360 of late 1966, Hutchinson's article, giving Lehmer's method for the 7090 (a 36-bit word machine, remember), appeared in CACM. The article gives (2*35)-31 as the largest prime less than 2*35, and gives the primitive roots 5*5 and 5*13, but does not give any criteria for choosing a primitive root, nor does it tell how to find one.
We come now to a murky period in our history. We do know that the implementers had decided to replace the dubious algorithm they had with Lehmer's superior one. It is easy to see how the prime (2*31)-1 was chosen as the modulus. The question that I have not been able to answer is, how was the primitive root 16807 (or 7*5) chosen? I have checked with Larry Breed and Dick Lathwell, who were the only people involved, and neither of them can recall this detail.
I do know that several years later, while I was working in the APL group at the IBM Research Center, an article appeared [5] that explored the question of determining a good random number generator for the System/360. Before the article appeared, one of the authors (I forget which) called me to ask about the APL random number generator. When he asked what method we used, I told him; and then he asked what primitive root we used. I said, "16807". There was a long silence at the other end of the telephone line; then, "How did you know?" I wasn't able to tell him, any more than I am able to tell you. The article gave the report of a group at the IBM Research Center which had just conducted an intensive study on the design of a good random number generator for 32-bit machines such as the IBM System/360. They undertook the study because some articles that had appeared in the computing literature had questioned the possibility of doing so. Their article said, yes, "a pseudo- random number generator for a 31-bit (sic) machine has to be chosen carefully. In particular only two values . . . of the many investigated gave test results as good as those obtained for . . . 7^5) The article contained no reference to the fact that the method and the particular constants recommended had been in use for over two years on APL\360 before their results were found.
The authors also wrote that they "are indebted to Dr. B. Tuckerman of the IBM Research Center for his invaluable help in finding the positive primitive root of the prime" (2^31)-1. By the way, Tuckerman at one time was the record holder in the large prime stakes. In 1971, using a 360/91 at IBM's Research Center, he found that (2*19937)-1 (a 6002-digit number) is prime. The search took 39 minutes and 26.4 seconds of computer execution time (Guinness Book of Records, you could look it up). I called Dr. Tuckerman just recently to see if he remembered the circumstances around his choice of the primitive root 16807. Although he didn't recall the precise circumstances surrounding his choice, he was very helpful in discussing the general question of how one goes about finding a primitive root for a given prime.
Two last pieces of history. The first APL\360 implementation of the random number generator was the same as it is today, with one small exception. The assembly language code had a test following the generation of the next random link to see whether the result was zero; if it was, 16807 was supplied in its place -- three instructions to look for a situation that can never occur! As we have seen, number theory guarantees that the residue on division by (2*31)-1 is always a positive integer less than (2*31)-1 and can never be zero! Three machine instructions were involved in the test, two that were always executed for each use of the random number generator, and one that could never be executed. I recently argued the case for removing this test after reviewing the code, and it has (or should have) been removed.
Next, the implementers of APL\360 (Breed, Lathwell, and Roger Moore, principally) carried on a low-key battle with Adin Falkoff and Iverson in this way: while the implementation was going on at breakneck speed, Falkoff and Iverson were writing a manual for using the system, a pamphlet they called "The APL Terminal System: Instructions for Operation". The first version of this pamphlet was dated 4 November 1966. A later version, dated 30 November 1966, had a dark brown cover bearing the cryptic legend "APL\360". The low-key battle I refer to concerned the claims made in the manual versus what was actually being implemented. One of these disparities had to do with the random functions. The roll function was not mentioned, and the deal function was described, although not by that name, in two forms. The first form corresponds to the current deal function. The second form permitted a vector right argument, so that N?V gave a result of shape (#V),N. Thus the following was permitted:
3?3 4 5 3 1 2 1 4 3 5 2 4
This last feature was never implemented. Episodes like this were surely behind Iverson's somewhat plaintive remarks in "The March on Armonk" [6]:
Along this line I would like to address a word to implementors, God bless them; (I say "God bless them" with the same sort of mixture of admiration, appreciation and exasperation as we say to the ladies, "God bless them"). It is really, in a way, the implementors that have control of a lot of this. If implementors decided they are going to provide a system which becomes popular with variations, then there is nothing any of us could do to control this.
All was not sweetness and light in the design of APL.
Finding a Primitive Root
How exactly does one find a primitive root of a prime? Daniel Shanks [7] writes:
Given a prime p, it is always possible to compute a primitive root by trial and error... but no general, explicit, nontentative method has been devised, and this, like a good criterion for primality, remains an important unsolved problem.
The way to proceed then is the following: for A and P positive integers, A: P - 1 , where P is a prime and R is one of its primitive roots, we get a vector called the indices of the powers of R mod P. I won't explain these further than to say that they are of great value in number theory: they are to the powers of the primitive root as logarithms are to exponentials.
There are three things about indices that are not well-known. The indices corresponding to a prime P form a permutation of length P-1. This permutation has these properties:
It is equal to its double downgrade (that is, it is a self-double-downgrading permutation, an sddp [9]: P -: /:/: P).
Its first difference, mod P-1, is a permutation.
This first difference permutation is also an sddp.
For example, the indices for the primitive root 3 mod 7 are 0 2 1 4 5 3. The double downgrade of this is also 0 2 1 4 5 3. Its first difference (mod 6) is 2 5 3 1 4, also a permutation. Its double downgrade (1-origin) is 2 5 3 1 4. The total number of permutations of length 6 is 720; there are only 48 sddp among these. The total number of permutations of length 5 is 120; there are only 8 sddp among these.
References
1. Iverson, K. E. J Introduction and Dictionary, p. 203. 2. Clark, C. E., and Holz, B. W. Tests of Randomness of the Bits of a Set of Pseudo-Random Numbers. Operations Research Office, (ASTIA Publication AB207553), Dec. 1958. 3. Hutchinson, D. W. A new uniform pseudo-random number generator. Communications of the ACM (June 1966). 4. Falkoff, A. D., Iverson, K. E., and Sussenguth, E. H. A formal description of System/360. IBM Systems Journal 3, 3 (1964), 198-261. 5. Lewis, T. A., Goodman, A. S., and Miller, J. M. A pseudo-random number generator for the System/360. IBM Systems Journal 8, 2 (1969) 136-45. 6. Iverson, K. The march on Armonk. Proceedings of the APL User's Conference, Binghamton, N. Y., July 1969, p. 50. 7. Shanks, D. Solved and Unsolved Problems in Number Theory. Spartan Books, Washington, D. C., 1962, P. 79. 8. Knuth, D. E. Seminumerical Algorithms. Addison-Wesley, Reading, Mass., 1969, Sec. 4.6-3. 9. McDonnell, E. E. Magic squares and permutations. APL Quote Quad 7, 3 (Fall 1976), 25-28.
This article originally appeared in APL Quote Quad Volume 8 Number 3 (March 1978) Page 42 and has been rewritten to use J notation in place of APL.
Notes
- "2147483647|16807^2147483646 is equal to 1" can be proved by Fermat's little theorem. It can also be verified using J with 16807 (2147483647&|@^) 2147483646.