A Tour (de Force) of APL in 16 Expressions
Roger K.W. Hui
 




0. Array Logic

(x>0)-(x<0)

Signum of x for real x , assuming the values 1 , 0 , or ¯1 according as x is strictly positive, 0, or strictly negative. The phrase dates from the earliest days of APL, being found in section 1.4 of A Programming Language [a6]. Three decades later, the idea was adopted by Knuth [a7], who wrote that “Iverson’s convention” led to improvements in exposition and technique. Some writers call these “data-driven conditionals” [a8]; during a discussion between Phil Last, John Scholes, and myself it was suggested that they be called “array logic” [a9].

The usefulness of array logic in APL is due to the following:

• Boolean functions have value 0 and 1 rather than true and false [a10, #implementers2].
Functions apply to entire arrays, as in, for example, +/x>100 to compute the number of elements of vector x greater than 100 [a10, #Maple].
A simple function precedence (“right to left”).



1. On Average

(+⌿÷≢) x

The average of x . The construct +⌿÷≢ is a fork, a 3-train. Using +⌿ andinstead of the more familiar +/ and is superior for reasons described in [a11].
 


2. Index-Of Selfie

x⍳x

Index-of with both arguments being the same, an instance of a “selfie” [a13]. x⍳x are like ID numbers; questions of identity on x can often be answered more efficiently on x⍳x than on x itself.

The x⍳x selfie occurs in several important and useful computations:

((⍳≢x)=x⍳x)⌿x nub (unique)
(x⍳x)∘.=(x⍳y) an efficient computation for x∧.=⍉y for matrices x and y [a13]
(x⍳x)∘.=(x⍳y) an efficient computation for x∘.≡y for non-simple vector x [a15]
(⍉↑x⍳¨x)⍳(⍉↑x⍳¨y)  inverted table index-of [a16]
(x⍳x) f⌸ y equivalent to x f⌸ y



3. Inverted Table Index-Of

(⍉↑x⍳¨x) ⍳ (⍉↑x⍳¨y)

A table is a set of values organized into rows and columns. The rows are records. Values in a column have the same type and shape. A table has a specified number of columns but can have any number of rows. The extended index-of on tables finds record indices.

   tx                   ty                    tx ⍳ ty
┌──────┬─┬───┬──┐    ┌──────┬─┬───┬──┐     3 1 5 2 5 5
│John  │M│USA│26│    │Min   │F│CN │17│
├──────┼─┼───┼──┤    ├──────┼─┼───┼──┤        tx ⍳ tx
│Mary  │F│UK │24│    │Mary  │F│UK │24│     0 1 2 3 4
├──────┼─┼───┼──┤    ├──────┼─┼───┼──┤
│Monika│F│DE │31│    │John  │M│UK │26│        ty ⍳ ty
├──────┼─┼───┼──┤    ├──────┼─┼───┼──┤     0 1 2 3 4 4
│Min   │F│CN │17│    │Monika│F│DE │31│
├──────┼─┼───┼──┤    ├──────┼─┼───┼──┤
│Max   │M│IT │29│    │Mesut │M│DE │24│
└──────┴─┴───┴──┘    ├──────┼─┼───┼──┤
                     │Mesut │M│DE │24│
                     └──────┴─┴───┴──┘

An inverted table is a table with the values of a column collected together. Comma-bar each (⍪¨) applied to an inverted table makes it “look” more like a table. And of course the columns have the same tally (). A table can be readily inverted and vice versa.

   x
┌──────┬─────┬───┬──────────────┐
│John  │MFFFM│USA│26 24 31 17 29│
│Mary  │     │UK │              │
│Monika│     │DE │              │
│Min   │     │CN │              │
│Max   │     │IT │              │
└──────┴─────┴───┴──────────────┘
   ⍪¨x                     ≢¨x      
┌──────┬─┬───┬──┐       5 5 5 5
│John  │M│USA│26│
│Mary  │F│UK │24│
│Monika│F│DE │31│
│Min   │F│CN │17│
│Max   │M│IT │29│
└──────┴─┴───┴──┘

   invert ← {↑¨↓⍉⍵}  
   vert   ← {⍉↑⊂⍤¯1¨⍵}
            
   x  ≡ invert tx
1
   tx ≡ vert x
1

A table has array overhead per element. An inverted table has array overhead per column. The difference that this makes becomes apparent when you have a sufficiently large number of rows. The other advantage of an inverted table is that column access is much faster.

An important computation is x index-of y where x and y are compatible inverted tables. Obviously, it can not be just x⍳y . The computation obtains by first “verting” the arguments (un-inverting the tables) and then applying , but often there is not enough space for that.

   ⍪¨x                   ⍪¨y
┌──────┬─┬───┬──┐     ┌──────┬─┬───┬──┐
│John  │M│USA│26│     │Min   │F│CN │17│
│Mary  │F│UK │24│     │Mary  │F│UK │24│
│Monika│F│DE │31│     │John  │M│UK │26│
│Min   │F│CN │17│     │Monika│F│DE │31│
│Max   │M│IT │29│     │Mesut │M│DE │24│
└──────┴─┴───┴──┘     │Mesut │M│DE │24│
                      └──────┴─┴───┴──┘
   x ⍳ y
4 4 4 4

   (vert x) ⍳ (vert y)
3 1 5 2 5 5

We derive a more efficient computation of index-of on inverted tables:

(vert x) ⍳ (vert y)  (a)
({⍉↑⊂⍤¯1¨⍵}x) ⍳ ({⍉↑⊂⍤¯1¨⍵}y)   (b)
(⍉↑⊂⍤¯1¨x) ⍳ (⍉↑⊂⍤¯1¨y)  (c)
(⍉↑x⍳¨x) ⍳ (⍉↑x⍳¨y)  (d)

(a)   The indices obtain by first uninverting the tables, that is, by first applying vert .
(b) Replace vert by its definition.
(c) Replace the D-fn by its definition. We see that ⊂⍤¯1 plays a major role. ⊂⍤¯1 encloses, or alternatively computes a scalar representation.
(d) For purposes of index-of x⍳¨x and x⍳¨y have the same information as ⊂⍤¯1¨x and ⊂⍤¯1¨y , but are much more efficient representations (small integers v the data itself).

Point (d) illustrated on column 0:

   ⊂⍤¯1⊢x0←0⊃x
┌──────┬──────┬──────┬──────┬──────┐
│John  │Mary  │Monika│Min   │Max   │
└──────┴──────┴──────┴──────┴──────┘
   x0 ⍳ x0
0 1 2 3 4

   ⊂⍤¯1⊢y0←0⊃y
┌──────┬──────┬──────┬──────┬──────┬──────┐
│Min   │Mary  │John  │Monika│Mesut │Mesut │
└──────┴──────┴──────┴──────┴──────┴──────┘
   x0 ⍳ y0
3 1 0 2 5 5

That is, the function {(⍉↑⍺⍳¨⍺)⍳(⍉↑⍺⍳¨⍵)} computes index-of on inverted tables.

I believe that in another language a derivation such as the one above would be very long (in part because the program would be very long), possibly impractically long.



Appeared in J in [a94] and in APL in [a95] and [a96].
 



4. Bar Chart

'.⎕'[x∘.>⍳⌈/x]

A bar chart of non-negative integer vector x .

   x← 3 1 4 1 5 9 
   '.⎕'[x∘.>⍳⌈/x]
⎕⎕⎕......
⎕........
⎕⎕⎕⎕.....
⎕........
⎕⎕⎕⎕⎕....
⎕⎕⎕⎕⎕⎕⎕⎕⎕

I encountered this phrase when I first learned APL in the mid 1970s. By then, it was part of the APL milieu, found in STATPACK [a17] and APL in Exposition [a18,p.13]. (I believe I saw it in the former.) The phrase brought home to me that indices can be entire arrays.

Anecdote

One of the best known epiphanies in Western culture is possibly that of Saul on the way to Damascus. My own epiphany in programming languages was much less dramatic than Saul’s. I was not going anywhere, there was no blinding light or voice from above, nor did I fall down. It was in 1965 or 1966, and I was sitting quietly in my office in the Department of Computing Science at the University of Alberta thinking about the commonly occurring statistical problem of classifying a number of observations given the left-hand end of the first frequency class, the class width, and the number of classes. This is a problem I had already handled in other languages but now I was trying to program it in a new language called Iverson’s notation or APL.

— Keith Smillie, Kenneth Iverson, APL and J:
Some Personal Recollections
, November 2004. [a19]
 



5. Permutations

Generate the matrix of all permutations of ⍳⍵ in lexicographic order.

Deriving a Solution

Commonly, in TDD (test-driven development) you start with a very simple case and trying to extend it to successively more general cases. It’s all too easy to be led into a dead-end because the simple case may have characteristics absent in a more general case. For myself, for this problem, I would start “in the middle”: Suppose I have perm 3 , obtained by whatever means:

   p
0 1 2
0 2 1
1 0 2
1 2 0
2 0 1
2 1 0

How do I get perm 4 from that? One way is as follows:

   p1←0,1+p
   (0 1 2 3[p1]) (1 0 2 3[p1]) (2 0 1 3[p1]) (3 0 1 2[p1])
┌───────┬───────┬───────┬───────┐
│0 1 2 3│1 0 2 3│2 0 1 3│3 0 1 2│
│0 1 3 2│1 0 3 2│2 0 3 1│3 0 2 1│
│0 2 1 3│1 2 0 3│2 1 0 3│3 1 0 2│
│0 2 3 1│1 2 3 0│2 1 3 0│3 1 2 0│
│0 3 1 2│1 3 0 2│2 3 0 1│3 2 0 1│
│0 3 2 1│1 3 2 0│2 3 1 0│3 2 1 0│
└───────┴───────┴───────┴───────┘

So it’s indexing each row of a matrix m by 0,1+p. There are various ways of forming the matrix m , one way is:

   ⍒⍤1∘.=⍨0 1 2 3
0 1 2 3
1 0 2 3
2 0 1 3
3 0 1 2

(Some authors [a70] waxed enthusiastic about this “magical matrix”.) In any case, a solution obtains readily from the above description: Form a matrix from the above individual planes; replace the 0 1 2 3 by ⍳⍵ ; and make an appropriate computation for the base case (when 0=⍵). See the 2015-07-12 subsection below.

The Best perm Function

What is the “best” perm function I can write in APL? This “best” is a benchmark not only on my own understanding but also on advances in APL over the years.

“Best” is a subjective and personal measure. Brevity comes into it but is not the only criteria. For example, {(∧/(⍳⍵)∊⍤1⊢t)⌿t←⍉(⍵⍴⍵)⊤⍳⍵*⍵} is short but requires space and time exponential in the size of the result, and that disqualifies it from being “best”. The similarly inefficient {(∧/(⍳⍵)∊⍤1⊢t)⌿t←↑,⍳⍵⍴⍵} is shorter still, but does not work for 1=⍵ .

• 1981, The N Queens Problem [a71]

    p←perm n;i;ind;t
   ⍝ all permutations of ⍳n
    p←(×n,n)⍴⎕io
    →(1≥n)⍴0
    t←perm n-1
    p←(0,n)⍴ind←⍳n
    i←n-~⎕io
   l10:p←(i,((i≠ind)/ind)[t]),[⎕io]p
    →(⎕io≤i←i-1)⍴l10

It was the fashion at the time that function be written to work in either index-origin and therefore have ⎕io sprinkled hither, thither, and yon.

• 1987, Some Uses of { and } [a28]

   perm:  ⍪⌿k,⍤¯1 (⍙⍵-1){⍤¯ 1 k~⍤1 0 k←⍳⍵
       :  1≥⍵
       :  (1,⍵)⍴0

Written in Dictionary APL [a44], wherein: ⍪⌿⍵ ←→ ⊃⍪⌿⊂⍤¯1⊢⍵ and differs from its definition in Dyalog APL;is equivalent toin dfns; ⍺{⍵ ←→ (⊂⍺)⌷⍵ ; and ¯ by itself is infinity.

• 1990-2007

I worked on perm from time to time in this period, but in J rather than in APL. The results are described in a J essay [a72] and in a Vector article [a70]. The lessons translate directly into Dyalog APL.

• 2008 [a73]

   pmat2←{{,[⍳2]↑(⊂⊂⎕io,1+⍵)⌷¨⍒¨↓∘.=⍨⍳1+1↓⍴⍵}⍣⍵⍉⍪⍬}

In retrospect, the power operator is not the best device to use, because the left operand function needs both the previous result (equivalent to perm ⍵-1) and . It is awkward to supply two arguments to that operand function, and the matter is finessed by computing the latter as 1+1↓⍴⍵ .

In this formulation, ⍉⍪⍬ is rather circuitous compared to the equivalent 1 0⍴0 . But the latter would have required aor similar device to separate it from the right operand of the power operator.

• 2015-07-12

   perm←{0=⍵:1 0⍴0 ⋄ ,[⍳2](⊂0,1+∇ ¯1+⍵)⌷⍤1⍒⍤1∘.=⍨⍳⍵}

For a time I thought the base case can be ⍳1 0 instead of 1 0⍴0 , and indeed the function works with that as the base case. Unfortunately (⍳1 0)≢1 0⍴0 , having a different prototype and datatype.

• Future

Where might the improvements come from?

  • We are contemplating an under operator whose monadic case is f⍢g ⍵ ←→ g⍣¯1 f g ⍵ . Therefore 1+∇ ¯1+⍵ ←→ ∇⍢(¯1∘+)⍵
  • Moreover, it is possible to define ≤⍵ as ⍵-1 (decrement) and ≥⍵ as ⍵+1 (increment), as in J; whence 1+∇ ¯1+⍵ ←→ ∇⍢≤⍵
  • Monadic = can be defined as in J, =⍵ ←→ (∪⍳⍨⍵)∘.=⍳≢⍵ (self-classify); whence ∘.=⍨⍳⍵ ←→ =⍳⍵

Putting it all together:

   perm←{0=⍵:1 0⍴0 ⋄ ,[⍳2](⊂0,∇⍢≤⍵)⌷⍤1⍒⍤1=⍳⍵}

We should do something about the ,[⍳2] ☺​.​​



Appeared in [a74].
 


6. Parentheses Nesting

+\1 ¯1 0['()'⍳x]

   x←'⍵((∇<S),=S,(∇>S))⍵⌷⍨?≢⍵'
   x ⍪ ⍉⍪ +\1 ¯1 0['()'⍳x]
⍵ ( ( ∇ < S ) , = S , ( ∇ > S ) ) ⍵ ⌷ ⍨ ? ≢ ⍵
0 1 2 2 2 2 1 1 1 1 1 2 2 2 2 1 0 0 0 0 0 0 0

The phrase has an ancient and storied history. As recounted by Alan Perlis [a20]:

I was at a meeting in Newcastle, England, where I’d been invited to give a talk, as had Don Knuth of Stanford, Ken Iverson from IBM, and a few others as well. I was sitting in the audience sandwiched between two very esteemed people in computer science and computing — Fritz Bauer, who runs computing in Bavaria from his headquarters in Munich, and Edsger Dijkstra, who runs computing all over the world from his headquarters in Holland.

Ken was showing some slides — and one of his slides had something on it that I was later to learn was an APL one-liner. And he tossed this off as an example of the expressiveness of the APL notation. I believe the one-liner was one of the standard ones for indicating the nesting level of the parentheses in an algebraic expression. But the one-liner was very short — ten characters, something like that — and having been involved with programming things like that for a long time and realizing that it took a reasonable amount of code to do, I looked at it and said, “My God, there must be something in this language.” Bauer, on my left, didn’t see that. What he saw or heard was Ken’s remark that APL is an extremely appropriate language for teaching algebra, and he muttered under his breath to me, in words I will never forget, “As long as I am alive, APL will never be used in Munich.” And Dijkstra, who was sitting on my other side, leaned toward Bauer and said, “Nor in Holland.” The three of us were listening to the same lecture, but we obviously heard different things.
 


7. Quicksort [a29, a30, a31]

Q←{1≥≢⍵:⍵ ⋄ (∇ ⍵⌿⍨0>s),(⍵⌿⍨0=s),∇ ⍵⌿⍨0<s←⍵ ⍺⍺ ⍵⌷⍨?≢⍵}

The operand ⍺⍺ is required to be a scalar function which returns ¯1 , 0 , or 1 , according to whetheris less than, equal to, or greater than . In Q , the phrase s←⍵ ⍺⍺ ⍵⌷⍨?≢⍵ compares each cell ofagainst a randomly-chosen “pivot” ⍵⌷⍨?≢⍵ ; the function then applies recursively to the elements of which are less than the pivot (0>s) and those which are greater than the pivot (0<s) .

Example with numbers:

   ⎕←x←?13⍴20
2 2 7 10 10 11 3 10 14 5 9 1 16

   (×-)Q x
1 2 2 3 5 7 9 10 10 10 11 14 16

   3(×-)4 3 ¯17
¯1 0 1

Example with strings:

   strcmp←{⍺≡⍵:0 ⋄ ¯1*</⍋↑⍺ ⍵}

   'syzygy' strcmp 'syzygy'
0
   'syzygy' strcmp 'eleemosynary'
1
   'syzygy' strcmp 'zumba'
¯1

   strcmp¨ Q 'syzygy' 'boustrophedon' 'kakistocracy' 'chthonic'
┌─────────────┬────────┬────────────┬──────┐
│boustrophedon│chthonic│kakistocracy│syzygy│
└─────────────┴────────┴────────────┴──────┘

Triplets

The variant Q1 encloses each of the three parts, resulting in an array with an interesting structure. The middle item of each triplet is the value of the pivot at each recursion. Since the pivot is randomly-chosen the results of Q1 can be different on the same argument.

   Q1←{1≥≢⍵:⍵ ⋄ (⊂∇ ⍵⌿⍨0>s),(⊂⍵⌿⍨0=s),⊂∇ ⍵⌿⍨0<s←⍵ ⍺⍺ ⍵⌷⍨?≢⍵}

   Q1 x
┌─────────┬─┬────────────────────────────────┐
│┌─┬───┬─┐│5│┌┬─┬───────────────────────────┐│
││1│2 2│3││ │││7│┌─┬────────┬──────────────┐││
│└─┴───┴─┘│ │││ ││9│10 10 10│┌┬──┬────────┐│││
│         │ │││ ││ │        │││11│┌┬──┬──┐││││
│         │ │││ ││ │        │││  │││14│16│││││
│         │ │││ ││ │        │││  │└┴──┴──┘││││
│         │ │││ ││ │        │└┴──┴────────┘│││
│         │ │││ │└─┴────────┴──────────────┘││
│         │ │└┴─┴───────────────────────────┘│
└─────────┴─┴────────────────────────────────┘
   Q1 x
┌────────────────────────────────┬──┬────────┐
│┌────────────────────┬────────┬┐│11│┌┬──┬──┐│
││┌──────────────┬─┬─┐│10 10 10│││  │││14│16││
│││┌─┬───┬──────┐│7│9││        │││  │└┴──┴──┘│
││││1│2 2│┌─┬─┬┐││ │ ││        │││  │        │
││││ │   ││3│5││││ │ ││        │││  │        │
││││ │   │└─┴─┴┘││ │ ││        │││  │        │
│││└─┴───┴──────┘│ │ ││        │││  │        │
││└──────────────┴─┴─┘│        │││  │        │
│└────────────────────┴────────┴┘│  │        │
└────────────────────────────────┴──┴────────┘

Order Statistic [a32; a33,§3.7]

S1←{⍵⌷⍨⍺⌷⍋⍵} computes the order statistic: 0 S1 ⍵ is the smallest cell ofand (¯1+≢⍵)S1 ⍵ is the biggest (and (2÷⍨¯1+≢⍵)S1 ⍵ is the median if ≢⍵ is odd).

   x←?(1+1e6)⍴2e9

   (⌊/x), 0 S1 x
87 87
   (⌈/x), (¯1+≢x)S1 x
1999998274 1999998274

   ⎕←m←(2÷⍨¯1+≢x)S1 x
999669657
   +/m<x
500000
   +/m>x
500000

A linear algorithm for the order statistic derives from the Quicksort strategy. The argument is partitioned into groups of cells less than a randomly-chosen pivot, equal to the pivot, and greater than the pivot. As above, the operand ⍺⍺ is required to be a scalar function which returns ¯1 , 0 , or 1 , according to whetheris less than, equal to, or greater than .

S←{
  s←⍵ ⍺⍺ p←⍵⌷⍨?≢⍵
  ⍺<  +/b←0>s:⍺ ∇ b⌿⍵
  ⍺<m←+/b←0≥s:p
  (⍺-m) ∇ (~b)⌿⍵
}

   0 (×-)S x
87
   (¯1+≢x) (×-)S x
1999998274
   (2÷⍨¯1+≢x) (×-)S x
999669657

S1 is less efficient than S because S1 sorts the entire argument, which S does not do. A version Sn for numeric vectors with a “hardwired” comparison is even faster:

Sn←{
  p←⍵⌷⍨?≢⍵
  ⍺<  +/b←⍵<p:⍺ ∇ b⌿⍵
  ⍺<m←+/b←⍵≤p:p
  (⍺-m) ∇ (~b)⌿⍵
}

   i←?≢x

   cmpx 'i S1 x' 'i (×-)S x' 'i Sn x'
  i S1 x    → 7.34E¯2 |   0% ⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕
  i (×-)S x → 3.53E¯2 | -52% ⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕
  i Sn x    → 1.13E¯2 | -85% ⎕⎕⎕⎕⎕


8. Golden Ratio

+∘÷\n⍴1

Convergents to the golden ratio by continued fractions [a45].

   +∘÷\16⍴1
1 2 1.5 1.66667 1.6 1.625 1.61538 1.61905 1.61765 1.61818
      1.61798 1.61806 1.61803 1.61804 1.61803 1.61803
   1 ∧ +∘÷\16⍴1
1 2 3 5 8 13 21 34 55 89 144 233 377 610 987 1597


9. Inner Product

x (+⌿ ×⍤¯1)⍤1 2 ⊢y

Inner product using row-at-a-time operation. For example:

   x←?3 4⍴100
   y←?4 5⍴100
   (x+.×y) ≡ x (+⌿ ×⍤¯1)⍤1 2 ⊢y
1

The main advantages of row-at-a-time over the conventional row-by-column are better cache utilization and ability to exploit sparsity [a47; a48; a49,#IC2013].
 


A. Cayley’s Theorem

r←,⍳⊢ ⋄ (r G) ≡ r ∘.{⍺[⍵]}⍨ ↓r G

The finite case of Cayley’s theorem [a28,§4.4]: every group G is isomorphic to a subgroup of the group of permutations. For example:

   ⎕←T←2 2∘⍴¨(4⍴2)∘⊤¨9 6 7 11 13 14
┌───┬───┬───┬───┬───┬───┐
│1 0│0 1│0 1│1 0│1 1│1 1│
│0 1│1 0│1 1│1 1│0 1│1 0│
└───┴───┴───┴───┴───┴───┘
   ⎕←G←∘.(≠.∧)⍨T
┌───┬───┬───┬───┬───┬───┐
│1 0│0 1│0 1│1 0│1 1│1 1│
│0 1│1 0│1 1│1 1│0 1│1 0│
├───┼───┼───┼───┼───┼───┤
│0 1│1 0│1 1│1 1│0 1│1 0│
│1 0│0 1│0 1│1 0│1 1│1 1│
├───┼───┼───┼───┼───┼───┤
│0 1│1 0│1 1│1 1│0 1│1 0│
│1 1│1 1│1 0│0 1│1 0│0 1│
├───┼───┼───┼───┼───┼───┤
│1 0│0 1│0 1│1 0│1 1│1 1│
│1 1│1 1│1 0│0 1│1 0│0 1│
├───┼───┼───┼───┼───┼───┤
│1 1│1 1│1 0│0 1│1 0│0 1│
│0 1│1 0│1 1│1 1│0 1│1 0│
├───┼───┼───┼───┼───┼───┤
│1 1│1 1│1 0│0 1│1 0│0 1│
│1 0│0 1│0 1│1 0│1 1│1 1│
└───┴───┴───┴───┴───┴───┘

T are the non-singular 2-by-2 boolean matrices, and G is its group table under boolean matrix multiplication. The function r←,⍳⊢ relabels a group table according to its index in its ravel (cf. x⍳x); such relabelling, like displaying in a different font, does not affect group-theoretic properties.

   r ← , ⍳ ⊢
   r G
0 1 2 3 4 5
1 0 4 5 2 3
2 3 5 4 1 0
3 2 1 0 5 4
4 5 3 2 0 1
5 4 0 1 3 2

The rows of r G are permutations, and the relabeled version of its group table as permutations is the same as the relabelling of the original group:

   (r G) ≡ r ∘.{⍺[⍵]}⍨ ↓r G
1


B. Symmetries of the Square

   D8
┌──┬──┬──┬──┬──┬──┬──┬──┐
│⊢ │⍒ │⍒⍒│⍋⌽│⌽ │⍋ │⍋⍒│⍒⌽│
├──┼──┼──┼──┼──┼──┼──┼──┤
│⍒ │⍒⍒│⍋⌽│⊢ │⍒⌽│⌽ │⍋ │⍋⍒│
├──┼──┼──┼──┼──┼──┼──┼──┤
│⍒⍒│⍋⌽│⊢ │⍒ │⍋⍒│⍒⌽│⌽ │⍋ │
├──┼──┼──┼──┼──┼──┼──┼──┤
│⍋⌽│⊢ │⍒ │⍒⍒│⍋ │⍋⍒│⍒⌽│⌽ │
├──┼──┼──┼──┼──┼──┼──┼──┤
│⌽ │⍋ │⍋⍒│⍒⌽│⊢ │⍒ │⍒⍒│⍋⌽│
├──┼──┼──┼──┼──┼──┼──┼──┤
│⍋ │⍋⍒│⍒⌽│⌽ │⍋⌽│⊢ │⍒ │⍒⍒│
├──┼──┼──┼──┼──┼──┼──┼──┤
│⍋⍒│⍒⌽│⌽ │⍋ │⍒⍒│⍋⌽│⊢ │⍒ │
├──┼──┼──┼──┼──┼──┼──┼──┤
│⍒⌽│⌽ │⍋ │⍋⍒│⍒ │⍒⍒│⍋⌽│⊢ │
└──┴──┴──┴──┴──┴──┴──┴──┘

A permutation can be represented as an integer vector or as a square boolean matrix with exactly one 1 in each row and each column. Functions pm←{⍵∘.=⍳≢⍵} and mp←⍳∘1⍤1 transform from one to the other.

   p                     x
6 3 2 1 5 4 0         96 84 59 5 19 47 36
   pm p      
0 0 0 0 0 0 1            x[p]
0 0 0 1 0 0 0         36 5 59 84 47 19 96
0 0 1 0 0 0 0
0 1 0 0 0 0 0            (pm p)+.×x
0 0 0 0 0 1 0         36 5 59 84 47 19 96
0 0 0 0 1 0 0
1 0 0 0 0 0 0
   mp pm p
6 3 2 1 5 4 0

⊢ ⍋ ⍒ ⌽ on permutations produce permutation results. They can be identified with ⊢ ⍉ ⊖⍉ ⊖ on square matrices.

   (⊢ pm p) ≡ pm ⊢p
1
   (⍉ pm p) ≡ pm ⍋p
1
   (⊖⍉pm p) ≡ pm ⍒p
1
   (⊖ pm p) ≡ pm ⌽p
1

Since ⊢ ⍉ ⊖⍉ ⊖ on matrices are transpositions of the square, then so are ⊢ ⍋ ⍒ ⌽ on permutations. The group table D8 is a compact presentation of numerous identities involving ⊢ ⍋ ⍒ ⌽ on permutations — D8[i;0] composed with D8[0;j] is D8[i;j] . For example:

i j D8[i;0] D8[0;j] D8[i;j]
5 5
2 2 ⍒⍒ ⍒⍒
1 2 ⍒⍒ ⍋⌽
1 5
1 6 ⍋⍒

That is, the 2 2 entry asserts that ⍒⍒⍒⍒p ←→ ⊢p ; the 1 5 entry asserts that ⍒⍋p ←→ ⌽p ; and so on. The veracity of these assertions can be checked by (⍎¨D8,¨'p') ≡ ⍎¨(∘.,⍨0⌷D8),¨'p' ⊣ p←?⍨n .

Finally, per the previous chapter where r←,⍳⊢ relabels the group elements, (r D8) ≡ r ∘.{⍺[⍵]}⍨ ↓r D8 .



Earlier versions of these ideas appeared in [a28, §4.4] and [a50].
 


C. Interval Index

I←{¯1+(r↓⍴⍵)⍴(≢⍺)↓i⊣i[i]←+\(≢⍺)>i←⍋⍺⍪,[r↓⍳⍴⍴⍵]⍵⊣r←1-⍴⍴⍺}

is a sorted (strictly ascending) vector; ⍺ I ⍵[i] is the least j such that ⍺[j] follows ⍵[i] in the ordering, or ≢⍺ if ⍵[i] follows ¯1↑⍺ in the ordering or ifis empty. The function is named classify in [a28,§1.2] and is similar to the I. primitive in J [a51, a52].

   x ← ¯1 2 3 7 8.5
   y ← 0 4 7 6 11 2 1 3 ¯5
   y ,[¯0.5] x I y
0 4 7 6 11 2 1 3 ¯5
0 2 3 2  4 1 0 2 ¯1

   x ← ↑ 'Fi' 'Jay' 'John' 'Morten' 'Roger'
   y ← ↑ 'JD' 'Jd' 'Geoff' 'Anna' 'Scott' 'Zeus  '

   x I y
0 1 0 ¯1 4 4

   x I x
1 2 3 4 5

¯∞       Fi       Jay      John     Morten   Roger   ∞
↓        ↓        ↓        ↓        ↓        ↓       ↓
(       )[       )[       )[       )[       )[       )
    ¯1       0        1        2        3        4    
   Anna    Fi        Jay      John    Morten   Roger
           JD        Jd                        Scott
           Geoff                               Zeus

D. 50847534

I was re-reading A Mathematician’s Apology [a1] before recommending it to Suki Tekverk, our summer intern, and came across a statement that the number of primes less than 1e9 is 50847478 (§14, page 23). The function pco [a89] from the dfns workspace does computations on primes; ¯1 pco n is the number of primes less than n :

   )copy dfns pco

   ¯1 pco 1e9
50847534

The two numbers 50847478 and 50847534 can not both be correct. A search of 50847478 on the internet reveals that it is incorrect and that 50847534 is correct. 50847478 is erroneously cited in various textbooks and even has a name, Bertelsen’s number, memorably described by MathWorld [a127] as “an erroneous name erroneously given the erroneous value of π(109) = 50847478”.

Although several internet sources give 50847534 as the number of primes less than 1e9 , they don’t provide a proof. Relying on them would be doing the same thing — arguing from authority — that led other authors astray, even if it is the authority of the mighty internet. Besides, I’d already given Suki, an aspiring mathematician, the standard spiel about the importance of proof.

How do you prove the 50847534 number? One way is to prove correct a program that produces it. Proving pco correct seems daunting — it has 103 lines and two large tables. Therefore, I wrote a function from scratch to compute the number of primes less than , a function written in a way that facilitates proof.

sieve←{
  b←⍵⍴{∧⌿↑(×/⍵)⍴¨~⍵↑¨1}q←2 3 5
  b[⍳6⌊⍵]←(6⌊⍵)⍴0 0 1 1 0 1 
  49≥⍵:b
  p←(≢q)↓⍸∇⌈⍵*0.5
  m←1+⌊(⍵-1+p×p)÷2×p
  b ⊣ p {b[⍺×⍺+2×⍳⍵]←0}¨ m
}

   +/ sieve 1e9
50847534

sieve ⍵ produces a boolean vector b with length such that b/⍳⍴b are all the primes less than . It implements the sieve of Eratosthenes: mark as not-prime multiples of each prime less than ⍵*0.5 ; this latter list of primes obtains by applying the function recursively to ⌈⍵*0.5 . Some obvious optimizations are implemented:

  Multiples of 2 3 5 are marked by initializing b with ⍺⍴{∧⌿↑(×/⍵)⍴¨~⍵↑¨1}2 3 5 rather than with ⍺⍴1 .
Subsequently, only odd multiples of primes > 5 need to be marked.
Multiples of a prime to be marked can start at its square.

Further examples:

   +/∘sieve¨ ⍳10
0 0 0 1 2 2 3 3 4 4

   +/∘sieve¨ 10*⍳10
0 4 25 168 1229 9592 78498 664579 5761455 50847534

There are other functions which are much easier to prove correct, for example, siev←{2=+⌿0=(⍳3⌈⍵)∘.|⍳⍵} . However, siev requires O(⍵*2) space and siev 1e9 can not produce a result with current technology. (Hmm, a provably correct program that can not produce the desired result …)

We can test that sieve is consistent with pco and that pco is self-consistent. pco is a model of the p: primitive function in J [a52]. Its cases are:

   pco n the n-th prime
¯1 pco n the number of primes less than n
 1 pco n 1 iff n is prime
 2 pco n the prime factors and exponents of n
 3 pco n the prime factorization of n
¯4 pco n the next prime < n
 4 pco n the next prime > n
10 pco r,s   a boolean vector b such that r+b/⍳⍴b are the primes in the half-open interval [r,s) .
   ¯1 pco 10*⍳10      ⍝ the number of primes < 1e0 1e1 ... 1e9
0 4 25 168 1229 9592 78498 664579 5761455 50847534

   +/ 10 pco 0 1e9    ⍝ sum of the sieve between 0 and 1e9
50847534
                      ⍝ sum of sums of 10 sieves 
                      ⍝ each of size 1e8 from 0 to 1e9
   +/ t← {+/10 pco ⍵+0 1e8}¨ 1e8×⍳10
50847534
                      ⍝ sum of sums of 1000 sieves 
                      ⍝ each of size 1e6 from 0 to 1e9
   +/ s← {+/10 pco ⍵+0 1e6}¨ 1e6×⍳1000
50847534

   t ≡ +/ 10 100 ⍴ s
1
                      ⍝ sum of sums of sieves with 1000 randomly
                      ⍝ chosen end-points, from 0 to 1e9
   +/ 2 {+/10 pco ⍺,⍵}/ 0,({⍵[⍋⍵]}1000?1e9),1e9
50847534

   ¯1 pco 1e9         ⍝ the number of primes < 1e9
50847534
   pco 50847534       ⍝ the 50847534-th prime
1000000007
   ¯4 pco 1000000007  ⍝ the next prime < 1000000007
999999937
   4 pco 999999937    ⍝ the next prime >  999999937
1000000007
   4 pco 1e9          ⍝ the next prime > 1e9
1000000007

                      ⍝ are 999999937 1000000007 primes?
   1 pco 999999937 1000000007
1 1
                      ⍝ which of 999999930 ... 1000000009 
                      ⍝ are prime?
   1 pco 999999930+4 20⍴⍳80
0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1

   +/∘sieve¨ 999999937 1000000007 ∘.+ ¯1 0 1 
50847533 50847533 50847534
50847534 50847534 50847535



Appeared in [a128].
 


E. A Memo Operator

Jay sent out an e-mail asking who would be interested to go see the film The Man Who Knew Infinity [a54]. Fiona’s response: Ah, Mr. 1729; I’m in. John was game, and so was I.

 

Partition Function

Nik found an essay by Stephen Wolfram [a55] on Ramanujan written on occasion of the film, wherein Wolfram said:

A notable paper [Ramanujan] wrote with Hardy [a56] concerns the partition function (PartitionsP in the Wolfram Language) that counts the number of ways an integer can be written as a sum of positive integers. The paper is a classic example of mixing the approximate with the exact. The paper begins with the result for large n:

But then, using ideas Ramanujan developed back in India, it progressively improves the estimate, to the point where the exact integer result can be obtained. In Ramanujan’s day, computing the exact value of PartitionsP[200] was a big deal — and the climax of his paper. But now, thanks to Ramanujan’s method, it’s instantaneous:

In[1]:= PartitionsP[200]
Out[1]= 3 972 999 029 388

[In Wolfram Language AKA Mathematica, user input is indicated by In[.]:= and the corresponding system output by Out[.]= . For readability, groups of 3 decimal digits are delimited by a blank instead of the comma commonly used in the UK or USA.]

A partition [a57] of a non-negative integer n is a vector v of positive integers such that n = +/v, where the order in v is not significant. For example, the partitions of the first seven non-negative integers are:

┌┐
││
└┘
┌─┐
│1│
└─┘
┌─┬───┐
│2│1 1│
└─┴───┘
┌─┬───┬─────┐
│3│2 1│1 1 1│
└─┴───┴─────┘
┌─┬───┬───┬─────┬───────┐
│4│3 1│2 2│2 1 1│1 1 1 1│
└─┴───┴───┴─────┴───────┘
┌─┬───┬───┬─────┬─────┬───────┬─────────┐
│5│4 1│3 2│3 1 1│2 2 1│2 1 1 1│1 1 1 1 1│
└─┴───┴───┴─────┴─────┴───────┴─────────┘

Can we compute the partition counting function “instantaneously” in APL?

We will use a recurrence relation derived from the pentagonal number theorem, proved by Euler more than 160 years before Hardy and Ramanujan:

equation (11) in [a59]. In APL:

   pn  ← {1≥⍵:0≤⍵ ⋄ -/+⌿∇¨rec ⍵}
   rec ← {⍵ - (÷∘2 (×⍤1) ¯1 1 ∘.+ 3∘×) 1+⍳⌈0.5*⍨⍵×2÷3}

   pn¨ 0 1 2 3 4 5
1 1 2 3 5 7

   pn 30
5604

Warning: don’t try pn 200 because that would take a long time. Why? pn 200 would engender pn being applied to each element of rec 200 and:

   rec 200
199 195 188 178 165 149 130 108 83 55 24 ¯10
198 193 185 174 160 143 123 100 74 45 13 ¯22

Each of the non-negative values would itself engender further recursive applications.

A Memo Operator

I recalled from J the memo adverb (monadic operator) M. Paraphrasing the J dictionary [a60; a52]: f M is the same as f but may keep a record of the arguments and results for reuse. It is commonly used for multiply-recursive functions. Sounds like just the ticket!

The Functionistas had previously implemented a (dyadic) memo operator [a61]. The version here is more restrictive but simpler. The restrictions are as follows:

  • The cache is constructed “on the fly” and is discarded once the operand finishes execution.
  • The operand is a dfn {condition:basis ⋄ expression} where none of condition, base and expression contain aand recursion is denoted byin expression.
  • The arguments are scalar integers.
  • The operand function does not have ¯1 as a result.
  • A recursive call is on a smaller integer.

With these restrictions, our memo operator can be defined as follows:

M←{
 f←⍺⍺
 i←2+'⋄'⍳⍨t←3↓,⎕CR'f'
 0=⎕NC'⍺':⍎' {T←(1+  ⍵)⍴¯1 ⋄  {',(i↑t),'¯1≢T[⍵]  :⊃T[⍵]   ⋄ 
                                        ⊃T[⍵]  ←⊂',(i↓t),'⍵}⍵'
          ⍎'⍺{T←(1+⍺,⍵)⍴¯1 ⋄ ⍺{',(i↑t),'¯1≢T[⍺;⍵]:⊃T[⍺;⍵] ⋄ 
                                        ⊃T[⍺;⍵]←⊂',(i↓t),'⍵}⍵'
}

   pn M 10
42
   pn 10
42

Obviously, thelines in M are crucial. When the operand is pn, the expression which is executed is as follows; the characters inserted by M are shown in red.

{T←(1+⍵)⍴¯1 ⋄ {1≥⍵:0≤⍵ ⋄ ¯1≢T[⍵]:⊃T[⍵] ⋄ ⊃T[⍵]←⊂-/+⌿∇¨rec ⍵}⍵}

The machinations produce a worthwhile performance improvement:

   pn M 30
5604
   cmpx 'pn M 30'
3.94e¯4
   cmpx 'pn 30'
4.49e1 
   4.49e1 ÷ 3.94e¯4
113959

That is, pn M 30 is faster than pn 30 by a factor of 114 thousand.

Can APL compute pn M 200 “instantaneously”? Yes it can, for a suitable definition of “instantaneous”.

   cmpx 'm←pn M 200'
4.47e¯3
   m
3.973e12
   0⍕m
 3972999029388

How long would pn 200 have taken without a memo operator? A lower bound derives by counting the number of function calls:

   pn   ← {1≥⍵:0≤⍵ ⋄ -/+⌿∇¨rec ⍵}

   pnc  ← {1≥⍵:1 ⋄ 1++/∇¨,rec ⍵}    ⍝ # of function calls
   pnc1 ← {C←0 ⋄ C⊣{C+←1 ⋄ 1≥⍵:0≤⍵ ⋄ -/+⌿∇¨rec ⍵}⍵}

   pnc¨ ⍳15
1 1 5 9 17 29 49 87 149 261 451 781 1351 2333 4041
   pnc M¨ ⍳15
1 1 5 9 17 29 49 87 149 261 451 781 1351 2333 4041
   pnc1¨ ⍳15
1 1 5 9 17 29 49 87 149 261 451 781 1351 2333 4041

pnc and pnc1 both count the number of function calls required by pn . pnc is in the form {condition:basis ⋄ expression} that M can handle; pnc1 is not. Extrapolating that pn 30 took 45 seconds to compute:

   pnc M¨ 30 200
2.59912e7 7.5715e47
   ⊢ s ← 45 × ÷⍨/ pnc M¨ 30 200    ⍝ # of seconds
1.3109e42
   s ÷ ×/ 365.2425 24 60 60        ⍝ # of years
4.15407e34

The algorithm pn M can compute the number of partitions for large n rapidly and exactly if extended precision arithmetic were available. From J:

   timer 'pn M. 1000'
0.0524486

   pn M. 1000
24061467864032622473692149727991

M also works on operands which are anonymous, dyadic, or have non-scalar results.



Presented in [a62].
 


F. Euler’s Identity

 
Euler’s identity

   0 = 1 + * ○ 0j1
 

is the most beautiful equation in all of mathematics, relating in one short phrase the fundamental quantities 0, 1, e, π, and 0j1 and the basic operations plus, times, and exponentiation. It is a particular case of Euler’s formula

   (*0j1×⍵) = (cos ⍵) + 0j1 × (sin ⍵)

Herewith a proof of Euler’s formula and Euler’s identity.

*0j1×⍵    (!⍳∞) ÷⍨ e ← ∞ ⍴ 1 0j1  ¯1 0j¯1  1 0j1  ¯1 0j¯1
cos ⍵     (!⍳∞) ÷⍨ c ← ∞ ⍴ 1   0  ¯1    0  1   0  ¯1    0
sin ⍵     (!⍳∞) ÷⍨ s ← ∞ ⍴ 0   1   0   ¯1  0   1   0   ¯1

The table above lists three functions and their Taylor series coefficients. The symboldenotes infinity and is not (yet) implemented in APL. Per [a125], ⍳∞ is a vector where (⍳∞)[i] is i ; and ∞⍴x with x a finite vector, is a vector y such that y[i]=x[i|⍨≢x] .

For terms 0 2 4 6 8 … e is the same as for c , and for terms 1 3 5 7 … e is the same as for 0j1×s :

   e ⍪ c ⍪ ⍉⍪ 0j1×s
1 0j1  ¯1 0j¯1  1 0j1  ¯1 0j¯1  1 0j1
1   0  ¯1    0  1   0  ¯1    0  1   0 …
0 0j1   0 0j¯1  0 0j1   0 0j¯1  0 0j1
That is:
   e ⍪ ⍉⍪ c + 0j1×s
1 0j1  ¯1 0j¯1  1 0j1  ¯1 0j¯1  1 0j1 …
1 0j1  ¯1 0j¯1  1 0j1  ¯1 0j¯1  1 0j1

   e = c + 0j1×s
1 1 1 1 1 1 1 1 1 1 …

Whence

   (*0j1×⍵) = (cos ⍵) + 0j1 × (sin ⍵)

It is not too far a leap to plug in an extremal value of the functions:

   cos ○1
¯1
   sin ○1
1.22461e¯16
   (cos ○1) + 0j1 × (sin ○1)
¯1j1.22465e¯16

   * ○ 0j1
¯1

I used to think math was no fun,
’Cause I couldn’t see how it was done.
    Now Euler’s my hero,
    For I now see why 0
Equals e to the π i plus 1.

e raised to the π times i,
And plus 1 leaves you nought but a sigh,
    This fact amazed Euler,
    That genius toiler,
And still gives us pause, bye the bye.

              — Anon y Mous
 



Appeared in J in [a126].
 



created:  2016-07-25 07:25
updated: