Recreational APL

How the Roll Function Works

(in APL\360 and its Descendants)

If you look up roll in the IBM manual APL Language [1], you read The roll function employs an algorithm due to D. H. Lehmer: the result for each scalar argument X is a function of X and the random link variable ŒRL. The result of the roll function is system dependent but typically for X<2*31 is equal to ŒIO plus the integer part of X׌RL÷¯1+2*31. If you look up ŒRL in the same publication, you find that its value in a clear workspace is 7*5 , that its meaningful range is ¼¯2+2*31 , and that it is used in the roll and deal functions. If you had access to the inside of the computer, you would see that the following occurs when the roll function is executed: ŒRL„2147483647 | 16807׌RL Following this, the result is computed in accordance with the algorithm quoted from APL Language. 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|(¼7)°.*¼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|(¼6)°.*¼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. 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 I use "Ÿ" 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 ¼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 the algorithm quoted earlier from [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 current APL this would be ?B½2 (0-origin). 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 the current dyadic definition, this would be written (B†A½1)[B?B]. 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 is the current 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 and P an odd prime, A is a primitive root of P if and only if 1¬P|A*(P-1)÷Q , for any of the prime divisors Q of P-1. For example, the prime divisors of (2*31)-2 are 2, 3, 7, 11, 31, 151, and 331. (By the way, you should be able to determine that 2, 3, and 11 are divisors of 2,147,483,646 by applying the divisibility tests of grade school arithmetic, but if you knew Fermat's theorem, you would be able to tell immediately that 31 should also be a divisor. Fermat's theorem states that, for any prime N and any integer K, 0=N|(K*N)-K. In our case, we have N=31 and K=2, and thus 0=31|(2*31)-2) .) To determine if an integer A is a primitive root of (2*31)-1, all we need to do is test 1^.¬((2*31)-1)|A*((2*31)-2)÷Q . The powers are rather large, however, and it is clear we do not want to compute all the intermediate powers to obtain, say, A*6487866, where the exponent is ((2*31)-2)÷331. (This is the smallest exponent that has to be computed.) There are ways to compute powers of numbers efficiently, but I won't discuss them here, other than to say that in the case in point, we can raise a number to the 6,487,866th power with only 22 squarings and 14 multiplications by 6,487,866 (all of these modulo 2,147,483,647) [8]. Once we have found one primitive root, the remainder can be found fairly easily. To explain how this is so requires one more bit of theory. Let us go back to our elementary example concerning the primitive roots of 7. Suppose we had found that 3 is a primitive root of 7. To find the other primitive roots, obtain the numbers relatively prime to 6, namely 1 and 5. Then 7|3*1 5 gives all the primitive roots of 7, namely 3 and 5. If we take (P-1)|“P|R*¼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. Now here is a scoop for APL Quote Quad: there are three things about indices that I have not seen reported before. 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]. 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 (in 0-origin) 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. Some Problems What is the inverse of 16807, mod 2147483647? That is, what number multiplied by 16807, mod 2147483647, gives 1? This number is the penultimate random link. The ultimate link, of course, is 1. APL Language says, "for X<2*31", and describes the algorithm. What would you do if X‰2*31? Many people are introduced to the roll function with the example ?6 6 to show how the roll of a pair of dice can be simulated. The September-October 1977 issue of Personal Computing had an interesting dice-rolling game that lends itself nicely to APL modeling and also to a clever APL solution. The game is called Petals Around the Rose, and that name is significant. Every answer is an even number. That, and the result of each throw of the dice, are all that can be told to those trying to guess the result of the roll. The roller of the dice (one who knows the game) is called the Potentate of the Rose. Five dice are rolled, and the answer is given. The object of the game is for the other players to give the answer before the Potentate. They are supposed to give just the answer, not the theory behind the answer. If any player (other than the Potentate) gets the answer right on six successive rolls, this indicates that the player understands the game. The writeup in Personal Computing is quite amusing, and I suggest you locate a copy of the magazine to read about it. All that I will do here is to give 16 rolls and the corresponding answers, and leave to you the job of determining how to find the answer for a given throw. References APL Language. IBM Publication GC26-3847, Jan. 1976, p. 22. 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. Hutchinson, D. W. A new uniform pseudo-random number generator. Communications of the ACM (June 1966). Falkoff, A. D., Iverson, K. E., and Sussenguth, E. H. A formal description of System/360. IBM Systems Journal 3, 3 (1964), 198-261. 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. Iverson, K. The march on Armonk. Proceedings of the APL User's Conference, Binghamton, N. Y., July 1969, p. 50. Shanks, D. Solved and Unsolved Problems in Number Theory. Spartan Books, Washington, D. C., 1962, P. 79. Knuth, D. E. Seminumerical Algorithms. Addison-Wesley, Reading, Mass., 1969, Sec. 4.6-3. McDonnell, E. E. Magic squares and permutations. APL Quote Quad 7, 3 (Fall 1976), 25-28.