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 1Looking 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.
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 4This 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.
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.