J-ottings 13
Vector 14:1, page 33
Transcribed from the printed issue; not yet reviewed. Read it in the PDF of the issue, page 33
Howard Peelle recently gave alternative algorithms for the verbs which appeared in J-ottings 9 to obtain best rational approximations to any given fraction. I propose to show here that there is still lots of mileage in rational numbers! This article, following a section in the encyclopaedic volume "Concrete Mathematics" by Graham, Knuth and Patashnik [1], describes how all positive rational fractions can be generated by a scheme called the Stern-Brocot tree, named after the two mathematicians (French and German respectively) who discovered it. Stern-Brocot trees can be described in a remarkably concise way which nevertheless embraces a large number of J features.
Rational fractions are represented by a two-item vector of co-prime integers, representing numerator and denominator in that order, so that, for example, 1.5 is represented by 3 2. The algorithm starts with the two-item vector of fractions 0 1;1 0, which represent zero and infinity. At every step a new fraction is inserted between every pair of fractions currently in the list, the new fraction consisting of the sums of the neighbouring numerators and denominators. Thus after one step, the result sequence is 0 1; 1 1; 1 0, and after two steps it is 0 1; 1 2; 1 1; 2 1; 1 0.
This is an algorithm which has a succinct realisation in J. The first subproblem is to compute the inserted vector of fractions at any stage, that is to add the values of all neighbours — the familiar 2+/ of APL2 which becomes 2&(+/\) in a tacit J expression. The new layer so generated is
lyr=.<"1@(2&(+/\)@:>) NB. generate next layer of Stern-Brocot tree
lyr v=.0 1;1 1;1 0
┌───┬───┐
│1 2│2 1│
└───┴───┘
The @: is worthy of comment. If @ were used, the verb 2&(+\) would be applied to each of the three rank 0 items of v individually to produce the three-item vector 1 2 3, with each item boxed. The leftmost <"1 is to ensure that the result is a list of fractions and not a matrix.
The new layer is then merged with the old:
malt=./:@,&(i.@#){, NB. merge x. and y., alternate elements from each
v malt lyr v
┌───┬───┬───┬───┬───┐
│1 2│1 2│1 1│2 1│1 0│
└───┴───┴───┴───┴───┘
A single step of the algorithm is then
sb=.malt lyr NB. single Stern-Brocot step
m0
┌───┬───┐
│0 1│1 0│
└───┴───┘
sb m0
┌───┬───┬───┐
│0 1│1 1│1 0│
└───┴───┴───┘
At this point the power of power comes into play. A tree n layers deep is generated by sb^:n, for example
sb^:4 m0
┌───┬───┬───┬───┬───┬───┬───┬───┬───┬───┬───┬───┬───┬───┬───┬───┬───┐
│0 1│1 4│1 3│2 5│1 2│3 5│2 3│3 4│1 1│4 3│3 2│5 3│2 1│5 2│3 1│4 1│1 0│
└───┴───┴───┴───┴───┴───┴───┴───┴───┴───┴───┴───┴───┴───┴───┴───┴───┘
There are never any fractions requiring cancellation, and by generating sufficient layers, any predetermined rational fraction will ultimately be generated.
The above may remind readers of a somewhat better known pattern of the rational fractions called the Farey series. The principle of this is to generate rational fractions as the lower triangular part of the table t%/t=.1+i.n which, in rational arithmetic, looks like
Again all the rational fractions will eventually appear, but this time cancellation has to be performed, and duplicates eliminated, for example, there is only one distinct item on the diagonal.
Remarkably, Farey series are obtained by two small adjustments to the Stern-Brocot tree, one being to use a new start value which ensures that only one half of the tree, namely the fractional half, is used, and the other a step to eliminate denominators which exceed the chosen size for the Farey square. Cancellation and removal of duplicates has already been achieved at the Stern-Brocot stage. Another fork is called for, this time with copy as the central prong. All the numerators in a Stern-Brocot tree are found by 1&{&.>, and the copy vector is found by the hook >: >.
m1=.0 1;1 1 NB. start vector for fractions half of S-B tree
fy=.(>: >@(1&{&.>))#] NB. n fy sb^:(n-1) m1 is Farey nos, n=max denom
6 fy sb^:5 m1
┌───┬───┬───┬───┬───┬───┬───┬───┬───┬───┬───┬───┬───┐
│0 1│1 6│1 5│1 4│1 3│2 5│1 2│3 5│2 3│3 4│4 5│5 6│1 1│
└───┴───┴───┴───┴───┴───┴───┴───┴───┴───┴───┴───┴───┘
A Stern-Brocot sequence, thought of as a tree, is a strict binary tree, and so is traversed from its root by alternate left and right branches. The value at any node can be found very simply by associating left and right branching with the two matrices l and r given below, and applying the ordinary matrix multiplication rule:
ip=.+/ .* NB. inner product for matrices
l=.2 2$1 1 0 1 NB. l and r = left and right search matrices for ..
r=.2 2$1 0 1 1 NB. .. a S-B tree; use e.g. +/>(ip&.>)/l;r;r;l
The fraction to be found by following the branches left, right, right, left is:
Given a number, rational or not, it is interesting to obtain its two closest bounds within any given Stern-Brocot tree. This is easily achieved by first converting the tree to its numeric values:
nvs=.>@(%/&>) NB. Numerical Values of items in S-Brocot tree
nvs sb^:3 m0
0 0.3333333 0.5 0.6666667 1 1.5 2 3 _
and then finding the indexes of the lowest fraction in the tree which exceeds it in value, and the greatest fraction which does not.
nbrs=.+&_1 0@(+/@(>nvs)){] NB. x. nbrs y. is best approx to x. in tree y.
k=.^1
k nbrs sb^:3 m0
┌───┬───┐
│2 1│3 1│
└───┴───┘
It is not a requirement that a Stern-Brocot tree start with the "natural" vector pair 0 1;1 0 — any pair of nodes which exist in the tree can be taken as argument. This means that successive rational approximations can be made using nodes further and further down the tree. For example, a succession of approximations to pi can be obtained by observing that the values 3 1 and 7 2 occur in the primary tree sb^:5 m0.
Execute
and observe that the familiar approximation 22 7 appears in p.
Now try repeatedly
and obtain at every step closer and closer rational approximations to pi. Alternatively the left and right branch technique described above can be used to write a verb which obtains the first rational approximation to achieve a given degree of accuracy.
approx=.3 :0
:
p=.y.nbrs sb^:5 m0
while. x.<<./|y.-nvs p do.
p=.y. nbrs sb^:5 p end.
p=:(t=<./t=.|y.-nvs p)#p
)
0.00001 approx o.1
355 113
Here in conclusion is the complete set of the nouns and verbs used:
malt=.(/:@,&(i.@#)){, NB. merge x. and y., alternate elements from each
lyr=.<"1@(2&(+/\)@:>) NB. generate next layer of Stern-Brocot tree
m0=.0 1;1 0 NB. start vector for full S-Brocot tree
sb=.malt lyr NB. sb^:n m0 is Stern-Brocot tree n deep
m1=.0 1;1 1 NB. start vector for fractions half of S-B tree
fy=.(>: >@(1&{&.>))#] NB. n fy sb^:(n-1) m1 is Farey nos, n=max denom
nvs=.>@(%/&>) NB. Numerical Values of items in S-Brocot tree
nbrs=.+&_1 0@(+/@(>nvs)){] NB. x. nbrs y.=best approxns to x. in tree y.
l=.2 2$1 1 0 1 NB. l and r =left and right search matrices for..
r=.2 2$1 0 1 1 NB. .. a S-B tree; use e.g. +/>(ip&.>)/l;r;r;l
ip=.+/ .* NB. Matrix Inner Product
Quite a small amount of J to achieve all possible rational numbers!
Reference¶
[1] Concrete Mathematics - A Foundation for Computer Science: R.L. Graham, D.E. Knuth, O. Patashnik. Addison-Wesley 1989.
Technical Note¶
From: Norman Thomson, April 1997
Chris Burke has asked me to point out that my attribution to him of the J Phrase book is misleading (see Vector 13 No.3 p.117). The principal author of the J Phrase Book is Roger Hui; Ken Iverson and Eugene McDonnell were also major contributors, and other contributions were made by Chris Burke and Donald McIntyre. Also Roger is revising the Phrase Book to correct known errors, and to accommodate language changes in J3.03.
Chris would also like me to point out that à propos the article on windows programming in J-ottings 12, the windows driver was completely revised in J3, and is now much simpler. In particular wd 'wait' is unnecessary, and not recommended.