At Play with J: Extended Integers
Vector 13:3, page 127
Transcribed from the printed issue; not yet reviewed. Read it in the PDF of the issue, page 127
Extended Integers¶
J has recently had added to it a new class of number, called extended integer. An extended integer can be produced by applying the extend verb x: to an ordinary integer. An extended integer is displayed as an ordinary integer terminated with an x, like this:
An extended integer may also be written directly by putting an x at the end of an ordinary integer:
If one or more of the integers in a list is extended, they are all extended:
Various primitives produce extended integer results if the argument is extended. For example, very large exact factorials are possible:
Some verbs f signal domain errors on some extended arguments because the result is not integral; however, <.@f and >.@f will work on extended arguments. I think you’ll get the idea from a few examples:
1234 % 5x
|domain error
| 1234 %5x
1234 <.@% 5x
246x
1234 % 2x
617x
1234 <.@% 2x
617x
] r =. <.@%: 2*10^38x
14142135623730950488x
*: r + ,. _1 0 1
199999999999999999971238085416445537169x
199999999999999999999522356663907438144x
200000000000000000027806627911369339121x
What if you want to turn an extended integer into an ordinary integer? Those of you conversant with J will guess that the way to do this is to apply the inverse of x:, like this:
The Application¶
I was excited when I heard about extended integers, because I had an immediate use for them. Jeffrey Shallit gave a paper at APL83 in which he discussed the problem of determining how many times the random number generator had been used, given a value of the random link. In order to give an APL solution he had to include in the paper a portion of an extended arithmetic package. This was because the numbers needed to solve the problem were very large integers. Stating the problem as simply as possible, given the equation:
the problem is to find x, where y, M, and g are known. This is the basis of the random number generator we shall be discussing. The value g is also known as the generator, and this can lead to confusion when we talk about the generator of the generator. I offer my apology in advance. The problem is sometimes known as the logarithm problem, since (forgetting the M-residue for a moment), if we have:
and know y and g, we can find x by taking the base-g logarithm of y:
In solving (A) for the particular problem of APL\360 and its descendants (including J), x can be as large as 2,147,483,646. For this value of x, g^x has over 9,000,000,000 decimal digits, and would take several hundred large volumes to print out. There are tactics one can employ to cut the size of the problem down, but extended-precision arithmetic will still be required. In J’s implementation, the phrase M&|@^ is recognized and can be computed efficiently both in time and space.
The technique discussed by Shallit to solve the problem is due to Pohlig and Hellman. You’ll have to look up the references if you are interested in the mathematical background, since I shall focus on the problem’s algorithmic aspects.
The Gory Details¶
Several global constants are needed. The modulus used in random number generators of the APL\360 kind must be a prime. The largest prime that can be stored as a 4-byte integer is in fact also the largest integer that can be stored, that is, one less than 2^31. This prime was discovered by Euler and for over a hundred years was the largest prime known. It is the Mersenne prime M31, too, for those of you interested in the Euclidian perfect numbers.
It is convenient to have the value of the integer one less than M handy:
For the random number generator to have maximum period, the generator g must be a primitive root of the modulus. A primitive root of a prime has the property that its powers, mod the prime, are distinct. For example, the prime 7 has 3 and 5 as primitive roots, because their powers, mod 7, are distinct:
but the other positive integers less than 7 have repeated elements:
Dr. Bryant Tuckerman, of the IBM Watson Research Laboratory in Yorktown Heights, New York, gave the APL\360 implementors the primitive root 7^5, or 16807. A decade or so later people began exploring random number generators extensively, and were surprised to find that there was no better generator than this for the modulus M:
A prime-power factor is an integer all of whose factors are the same. For example, 32, 9, 125, 49, and 11 are all prime-power factors. The prime-power factors of L play a key role in the algorithm. The primitive q: in J yields the prime factors of a number, but these may be repeated. For example, q: 12 is 2 2 3. The algorithm requires that repeated primes be replaced by their product. The verb h, to be defined later, factors numbers and replaces repeated items by their product:
Certain powers of the generator g are needed. Those needed are the quotients of dividing L by f, and multiplying this quotient by the integers less than the items of f. For example, for the factor 7 we get:
and similarly for the other factors.
The verb p, to be defined later, raises the generator g to to any integer power, mod M. We use it to raise the generator to the powers B:
Such a list is made for each prime-power factor. These are boxed and joined together, forming q, a list of lists, containing 542 numbers altogether (+/2 9 7 11 31 151 331), and too large to display here.
You should be warned that the formation of q takes a minute or so to execute, depending on the speed of your computer. I find it convenient to comment out this line in the script, and insert the value of q directly.
We need the quotient of L with its factors as a separate global noun:
Those are all of the global nouns. Now we have to deploy a number of utility verbs. The verb w:
raises each item of its argument’s nub to its tally:
The verb h:
factors its argument and produces the prime power factors from it:
The verb s:
raises its left argument x to the power of its right argument, mod M:
The verb p:
raises g to the power of its argument:
The verb j:
divides its argument into L, and multiplies this quotient by the non-negative integers less than it.
The Main Problem¶
With all this behind us, we’re ready to discuss the main problem. Suppose we find that the value of y, the random link, is 1209311799:
We define a verb t:
The phrase
can be replaced by
,.D =. y s e NB. let this be (D)
2147483646x
473297587x
1537170743x
353622995x
4096x
709324280x
667991092x
The heart of the matter is that each distinct value that y may take yields a different list D. We look for the index of each item of D in q, finding:
These values are the residues we seek. We form a table F by stitching f and E:
and this is also the result of applying verb t to y:
If you refer to Hui’s article, you’ll see that F is similar to the table shown on page 64, beneath the expression c mr q. That is, F is a table of moduli and residues. For example, the items in list D correspond to the factors in f. In particular, the value 1537170743x corresponds to the factor 7. We gave all the possible values that may be taken on in this position in list C back a bit. In C we find that the value 1537170743x is in position 5, and so in F we find the value 5 next to the 7 in the first column. It is a remarkable fact that the expression y s e produces values which must occur also in the corresponding item of q, and that its index in the list in q is the residue we want for the next step.
The verb r:
inserts Hui’s verb cr1 between each of the items of its argument, and yields a two-item list, with the first item necessarily equal to L, and the second item the power of g yielding y, mod M.
We take the tail of z as our desired x:
We can verify that x is indeed the desired value by applying p to it:
Are we happy? We shouldn’t be yet, because this wasn’t precisely the problem we wanted to solve, which was, how many times had the random number generator been used to arrive at the given random link. This part is easy, because we know that the initial value of the random link is 16807, which corresponds to exponent 1. All we have to do to get the value we want is to decrease x by 1. This gives us at last the verb ner:
References¶
Hui, Roger K. W., The Ball Clock Problem, Vector 12 2, October 1995, 55-66
Pohlig, Stephen C., and Martin E. Hellman, An Improved Algorithm for Computing Logarithms Over GF(p) and Its Cryptographic Significance, IEEE Trans. Info. Theory IT-24 (1978) 106-110
Shallit, J. O., Merrily We Roll Along: Some Aspects of ?, APL83 Conference Proceedings, APL Quote-Quad 13 3, March 1983, 243-248
Appendix¶
NB. Script for finding x in y=m|g^x, knowing y, M, and g,
NB. particularized for use in analyzing the behavior of
NB. the linear congruential random number generator found
NB. in APL\360 and its descendants.
NB. GLOBAL NOUNS
M =. x: <: 2^31
L =. <: M
g =. x: 7^5
f =. h L NB. f is 2 9 7 11 31 151 331
q =. <@p@j"0 f NB. comment this line out
NB. replace line above by value of q directly
e =. L % f
NB. UTILITY VERBS
w =. ~. ^ #/.
h =.w@q:
s=. M&|@^
p =. g&s
j =. L&% * i.
NB. HUI's CHINESE REMAINDER VERBS FROM VECTOR 12 2, P 66
NB. INCLUDING GCD AS A LINEAR COMBINATION
NB. Chinese Remainder
ab =. |.@(gcd/ * [ % +./)@(,&{.)
cr1 =. [: |/\ *.&{. , ,&{: +/ .* ab
chkc =. [: assert ,&{: -: ,&{. | {:@cr1
cr =. cr1 [ chkc
NB. GCD as a Linear Combination
g0 =. , ,. =@i.@2:
it =. {: ,: {. - {: * <.@%&{./
gcd =. (}.@{.)@(it^:(*@{.@{:)^:_)@g0
assert=. 13!:8@(12"_)^:-.
NB. MAIN VERBS
t =. f"_ ,. q"_ i.&> ] s e"_
r =. {: @ (cr1/ @ t)
ner =. <:@r NB. number of executions of roll