Unforgettable Numbers

Professor Kip Murray asked on the J Forum for the “unforgettable” times seen on a 24-hour digital clock.

The problem is that every number has something notable about it, so that each number is “unforgettable” and consequently it’s hard to remember any single one of them.

0000 all zeros
0001 first counting number
0002 first prime number
0003 first odd prime number
0004 first composite number
...

      ⎕rl←7*5 ⋄ 24 60⊤¯1+?×/24 60
6 59

0659 the time I wake up if the alarm was set at 0700 ☺

You may have heard of the following story about Hardy and Ramanujan. One day Hardy took a taxi to visit Ramanujan. On arriving, Hardy told Ramanujan that the taxi had the thoroughly unremarkable 4-digit number n on its licence plate. Ramanujan immediately remarked that n is the first number that … . I forget what n or the property was, something like n is the first number that can be written as the sum of two perfect cubes in two different ways, something typically Ramanujanian.

Yes, that was it:

      c←3*⍨⍳200           ⍝ cubes of numbers 1..200
      t←(∘.≤⍨⍳≢c) × ∘.+⍨c  
           ⍝ upper triangle of the addition table of these cubes

      e←(,t){⍺⍵}⌸,⍳⍴t     ⍝ unique sums and their indices
      e←(2=≢¨e[;2])⌿e     ⍝ sums that derive two different ways
      e
┌───────┬─────────────────┐
│1729   │┌────┬────┐      │
│       ││1 12│9 10│      │
│       │└────┴────┘      │
├───────┼─────────────────┤
│1092728│┌─────┬─────┐    │
│       ││1 103│64 94│    │
│       │└─────┴─────┘    │
├───────┼─────────────────┤
     ...
├───────┼─────────────────┤
│5472152│┌───────┬───────┐│
│       ││102 164│128 150││
│       │└───────┴───────┘│
└───────┴─────────────────┘

      e⌷⍨(⊢⍳⌊/)e[;1]      ⍝ the row with the smallest sum
┌────┬───────────┐
│1729│┌────┬────┐│
│    ││1 12│9 10││
│    │└────┴────┘│
└────┴───────────┘
      +/ 1 12 * 3         ⍝ check the sums
1729
      +/ 9 10 * 3
1729

Now that I have worked out the number I can find the story on the net. On hearing the story, J.E. Littlewood remarked that “every positive integer is one of Ramanujan’s personal friends”.

P.S. In my youth, when I needed to remember a (5-digit) number for a time, I computed its largest prime factor by mental calculation. Try it and you’ll see why it works.

Translated and updated from a J Wiki essay which first appeared on 2009-08-22.

Changes of Heart

Karen Shaw started the ball rolling (hearts afluttering?) by asking Jay Foad to come up with a one-liner for St. Valentine’s Day; he then solicited contributions from the language development group. Nick Nickolov responded with the following, with no explanation other than that there is room for improvement:

⎕io←0 ⋄ (⊢,⌽)' X'[{(.5×n*2)>+/(⍺-.6×⍵)⍵*2}/¨0↓n-⍳(2×n),n←20]

           XXXXX        XXXXX           
         XXXXXXXXXX  XXXXXXXXXX         
        XXXXXXXXXXXXXXXXXXXXXXXX        
       XXXXXXXXXXXXXXXXXXXXXXXXXX       
       XXXXXXXXXXXXXXXXXXXXXXXXXX       
       XXXXXXXXXXXXXXXXXXXXXXXXXX       
      XXXXXXXXXXXXXXXXXXXXXXXXXXXX      
      XXXXXXXXXXXXXXXXXXXXXXXXXXXX      
      XXXXXXXXXXXXXXXXXXXXXXXXXXXX      
      XXXXXXXXXXXXXXXXXXXXXXXXXXXX      
       XXXXXXXXXXXXXXXXXXXXXXXXXX       
       XXXXXXXXXXXXXXXXXXXXXXXXXX       
       XXXXXXXXXXXXXXXXXXXXXXXXXX       
       XXXXXXXXXXXXXXXXXXXXXXXXXX       
        XXXXXXXXXXXXXXXXXXXXXXXX        
        XXXXXXXXXXXXXXXXXXXXXXXX        
        XXXXXXXXXXXXXXXXXXXXXXXX        
         XXXXXXXXXXXXXXXXXXXXXX         
         XXXXXXXXXXXXXXXXXXXXXX         
          XXXXXXXXXXXXXXXXXXXX          
           XXXXXXXXXXXXXXXXXX           
           XXXXXXXXXXXXXXXXXX           
            XXXXXXXXXXXXXXXX            
             XXXXXXXXXXXXXX             
             XXXXXXXXXXXXXX             
              XXXXXXXXXXXX              
               XXXXXXXXXX               
                XXXXXXXX                
                 XXXXXX                 
                   XX                   

Subsequently, Nick described how he derived the expression:

Between two meetings I saw Jay’s request for ideas. I sent him a link to the Heart Curve on MathWorld, but he thought we’d need graphics for that. To prove that ASCII art is good enough, I tried to plot the set of points (x,y) where (x2 + y2 – 1)3x2 y3 ≈ 0 but the interpreter crashed with a syserror for some reason. (I don’t recall what I last did and can’t reproduce it.) Since I had to type the whole thing again anyway and I hadn’t been too successful in visualising the solutions of the above equation, I decided to write it in a less beautiful but simpler way: draw half a circle, tilt it, and concatenate a mirror image to form a heart. The moment I got something that resembled a heart, I sent the email and went downstairs for coffee.

Meanwhile, I changed Nick’s expression in steps in the process of trying to understand it. ⎕io←0 throughout.

a.  (⊢,⌽)' X'[{(.5×n*2)>+/(⍺-.6×⍵)⍵*2}/¨0↓n-⍳(2×n),n←20]
Nick’s original expression.

b. ,∘⌽⍨' X'[{(.5×n*2)>+/(⍺-.6×⍵)⍵*2}/¨0↓n-⍳(2×n),n←20]
The fork (⊢,⌽) computes x,⌽x without use of a temporary variable. It is equivalent to ,∘⌽⍨.

c. ,∘⌽⍨' X'[{(.5×n*2)>+/(⍺-.6×⍵)⍵*2}/¨n-⍳(2×n),n←20]
The 0↓ is a no-op and is likely scaffolding left from the development process.

d. ,∘⌽⍨' X'[{(.5×n*2)>+/(⍺-.6×⍵)⍵*2}/¨n-⍳2 1×n←20]
(2×n),n←20 is equivalent to 2 1×n←20.

e. ,∘⌽⍨' X'[(.5×n*2)>{+/(⍺-.6×⍵)⍵*2}/¨n-⍳2 1×n←20]

The comparison with .5×n*2 is on the result of the {...} and is likely more efficient “outside” rather than explicitly applied to each scalar “inside”. More importantly, moving it outside avoids the use of a lexical global, an aesthetic infelicity (your opinion may differ).

At this step, it dawned on me that is applied to a 2-element vector, resulting in a 2×n by n matrix of 2-element integer vectors:

   ⍳2 1×n←20
┌───┬───┬───┬───┬───┬
│0 0│0 1│0 2│0 3│0 4│
├───┼───┼───┼───┼───┼
│1 0│1 1│1 2│1 3│1 4│ ...
├───┼───┼───┼───┼───┼
│2 0│2 1│2 2│2 3│2 4│
├───┼───┼───┼───┼───┼
│3 0│3 1│3 2│3 3│3 4│
├───┼───┼───┼───┼───┼
     ...

f. ,∘⌽⍨' X'[(.5×n*2)>⊃∘.{+/(⍺-.6×⍵)⍵*2}/n-⍳¨2 1×n←20]
The reduction on each 2-element vector is equivalent to an outer product; that is, f/¨⍳a,b is equivalent to ⊃∘.f/⍳¨a,b. Not much of a step forward but facilitates what comes next.

g. ,∘⌽⍨' X'[(.5×n*2)>(n-⍳2×n)∘.{+/(⍺-.6×⍵)⍵*2}n-⍳n←20]
Applying twice removes the reduction and makes the outer product more straightforward.

h. ,∘⌽⍨' X'[(n×*⍨.5)>(n-⍳2×n)∘.{|⍵+0j1×⍺-.6×⍵}n-⍳n←20]
The function {+/(⍺-.6×⍵)⍵*2} computes the sum-of-squares on ⍺-.6×⍵ and . It is equally the square of the magnitude of a complex number whose real and imaginary parts are the two expressions. The subsequent comparison to .5×n*2 has the same outcome if instead we compare the square root of both sides, that is, compare n×*⍨.5 and the magitude (not the square of the magnitude) of the complex number.

i. (⎕ucs 32 9829)[(n×*⍨.5)>i∘.{|⍵+0j1×⍺-.6×⍵}|i←n-⍳2×n←20]

A couple of ideas here: (0) Instead of using X in the picture, use the heart symbol, Unicode codepoint 9829. (1) Instead of creating half a heart and then stitching the two halves together with ,∘⌽⍨, the whole heart can be created using an appropriate right argument to the outer product.

This result differs slightly from the previous ones, but makes a prettier picture.

j. ' ♥'[(n×*⍨.5)>∘.(|1 ¯.6j1+.×,)∘|⍨n-⍳2×n←20]

I remembered that Unicode codepoints can be included directly in an APL expression, so ⎕ucs 32 9829 can be replaced by ' ♥'. (Much to my relief. Talk about lack of aesthetics!)

The arguments to the outer product are i and |i, and i can be eliminated altogether by doing f∘|⍨ instead of i f |i.

The function {|⍵+0j1×⍺-.6×⍵}, which computes the magnitude of a complex number, remains the same if the real and imaginary parts are switched, and in this case it is advantageous to switch them for brevity. Moreover, the computation can be coded as a train involving an inner product.

k. ' ♥'[(n×*⍨.5)>|∘.{⍺-⍵×.6j1}∘|⍨n-⍳2×n←20]
Sometimes, trains are better coded as dfns. The magnitude | is moved “outside”, and is unchanged if the complex coefficient is replaced by the negation of its conjugate and the operation changed from + to -.

l. ' ♥'[(n÷2*÷2)>|i∘.-.6j1×|i←n-⍳2×n←20]
The expression can be shortened by using the variable i after all (last seen in step i). Replace n×*⍨.5 by n÷2*÷2 and voilà, a retro expression without dfns, trains, or any of them new-fangled operators.

In summary

a.     (⊢,⌽)' X'[{(.5×n*2)>+/(⍺-.6×⍵)⍵*2}/¨0↓n-⍳(2×n),n←20]
b.     ,∘⌽⍨' X'[{(.5×n*2)>+/(⍺-.6×⍵)⍵*2}/¨0↓n-⍳(2×n),n←20]
c.     ,∘⌽⍨' X'[{(.5×n*2)>+/(⍺-.6×⍵)⍵*2}/¨n-⍳(2×n),n←20]
d.     ,∘⌽⍨' X'[{(.5×n*2)>+/(⍺-.6×⍵)⍵*2}/¨n-⍳2 1×n←20]
e.     ,∘⌽⍨' X'[(.5×n*2)>{+/(⍺-.6×⍵)⍵*2}/¨n-⍳2 1×n←20]
f.     ,∘⌽⍨' X'[(.5×n*2)>⊃∘.{+/(⍺-.6×⍵)⍵*2}/n-⍳¨2 1×n←20]
g.     ,∘⌽⍨' X'[(.5×n*2)>(n-⍳2×n)∘.{+/(⍺-.6×⍵)⍵*2}n-⍳n←20]
h.     ,∘⌽⍨' X'[(n×*⍨.5)>(n-⍳2×n)∘.{|⍵+0j1×⍺-.6×⍵}n-⍳n←20]
i.     (⎕ucs 32 9829)[(n×*⍨.5)>i∘.{|⍵+0j1×⍺-.6×⍵}|i←n-⍳2×n←20]
j.     ' ♥'[(n×*⍨.5)>∘.(|1 ¯.6j1+.×,)∘|⍨n-⍳2×n←20]
k.     ' ♥'[(n×*⍨.5)>|∘.{⍺-⍵×.6j1}∘|⍨n-⍳2×n←20]
l.     ' ♥'[(n÷2*÷2)>|i∘.-.6j1×|i←n-⍳2×n←20]

The Diamond Kata

Acknowledgments

Morten Kromberg is the other co-author of this text but the blogging software prevents his being listed as such.

We are indebted to Jay Foad, Nick Nikolov, John Scholes, and Fiona Smith for comments on successive drafts of the MS.

The Problem

The diamond kata is a programming exercise used in the agile development, XP and TDD communities. [0, 1, 2, 3]. It first came to our attention in 2012 in discussions with Gianfranco Alongi of Ericcson AB, Sweden, and again more recently in [3]. The following problem description is from [2] (the description varies slightly in each of the above references).

We’re going to write a program that builds a diamond picture like this, for some set of letters from A to whatever:

--A--
-B-B-
C---C
-B-B-
--A--

Solutions in Dyalog APL

We are actually going to solve a slightly different problem: the function will have two arguments. The left argument is a scalar of the background element; the right argument is a vector of the “letters” to be used. The result described in the opening section can be produced as:

   '-' dia 'ABCD'
---A---
--B-B--
-C---C-
D-----D
-C---C-
--B-B--
---A---

(In an APL session, what you enter is indented and the resultant output is aligned at the left margin. As well, in this text, we sometimes present the indented-input/output pairs across the page.)

Making the argument the actual letters rather than the number of letters or the last letter, facilitates working with “alphabets” other than A-Z. For example:

   0 dia 3 1 4 2
0 0 0 3 0 0 0
0 0 1 0 1 0 0
0 4 0 0 0 4 0
2 0 0 0 0 0 2
0 4 0 0 0 4 0
0 0 1 0 1 0 0
0 0 0 3 0 0 0

The diamond result is readily seen to be composed of reflections of one of the quadrants. For example, '-' dia 'ABCD' (shown above) is composed of reflections of

A---
-B--
--C-
---D

The functions (reverse last axis) and (reverse first axis) provide the requisite reflections. If q is the quadrant shown above, then:

   ⌽q       q
---A     A---
--B-     -B--
-C--     --C-
D---     ---D

   ⌽⊖q      ⊖q
D---     ---D
-C--     --C-
--B-     -B--
---A     A---

There are various ways to stitch the quadrants together (and to drop the common middle row or column) to get the required result. The following is a terse way:

   f←{⍉⍵⍪1↓⊖⍵}

   f ⌽q          f f ⌽q        f⍣2 ⌽q
---D---       ---A---       ---A---
--C-C--       --B-B--       --B-B--
-B---B-       -C---C-       -C---C-
A-----A       D-----D       D-----D
              -C---C-       -C---C-
              --B-B--       --B-B--
              ---A---       ---A---

The same pattern (metakata?) of f⍣2 involving can be used to compute the minors of a matrix [4].

It remains to produce the quadrant q. We note that it is a diagonal matrix and resembles the identity matrix for which creation many techniques are known. We use two of the 34 different methods in [5].

The first method: For a vector with n letters, if each vector is followed by n copies of the background element, then the 2⍴n reshape of that rotates each row an additional step to the right, and the result is the required diagonal matrix.

   n←≢v←'ABCD'

   v,'-'⍴⍨s←2⍴n           s ⍴ v,'-'⍴⍨s←2⍴n
A----                  A---
B----                  -B--
C----                  --C-
D----                  ---D

The second method: for the proposed merge operator [6] ⍺ f merge ⍵ is an array like with at positions selected by f. As such, it is a functional version of selective assignment. Merge with 1 1∘⍉ as the operand does the trick.

   4 4⍴⍳16                1 1∘⍉ 4 4⍴⍳16
 1  2  3  4            1 6 11 16
 5  6  7  8
 9 10 11 12
13 14 15 16

   '-'⍴⍨2⍴n              'ABCD' (1 1∘⍉merge) '-'⍴⍨2⍴n
----                   A---
----                   -B--
----                   --C-
----                   ---D

Putting it together:

dia  ← {{⍉⍵⍪1↓⊖⍵}⍣2⌽s⍴⍵,⍺⍴⍨s←2⍴≢⍵}
diam ← {{⍉⍵⍪1↓⊖⍵}⍣2⌽⍵(1 1∘⍉merge)⍺⍴⍨2⍴≢⍵}


   '-' dia 'ABCD'         0 dia 3 1 4 2
---A---                0 0 0 3 0 0 0
--B-B--                0 0 1 0 1 0 0
-C---C-                0 4 0 0 0 4 0
D-----D                2 0 0 0 0 0 2
-C---C-                0 4 0 0 0 4 0
--B-B--                0 0 1 0 1 0 0
---A---                0 0 0 3 0 0 0

   '-' diam 'ABCD'        0 diam 3 1 4 2
---A---                0 0 0 3 0 0 0
--B-B--                0 0 1 0 1 0 0
-C---C-                0 4 0 0 0 4 0
D-----D                2 0 0 0 0 0 2
-C---C-                0 4 0 0 0 4 0
--B-B--                0 0 1 0 1 0 0
---A---                0 0 0 3 0 0 0

Further Examples

   ↑∘⎕a¨ 0 1 5 3           ⍝ four prefixes of the alphabet
┌┬─┬─────┬───┐
││A│ABCDE│ABC│
└┴─┴─────┴───┘
   '.' dia¨ ↑∘⎕a¨ 0 1 5 3  ⍝ apply dia to each prefix
┌┬─┬─────────┬─────┐
││A│....A....│..A..│
││ │...B.B...│.B.B.│
││ │..C...C..│C...C│
││ │.D.....D.│.B.B.│
││ │E.......E│..A..│
││ │.D.....D.│     │
││ │..C...C..│     │
││ │...B.B...│     │
││ │....A....│     │
└┴─┴─────────┴─────┘
   '.' {{⍉⍵⍪1↓⊖⍵}⍣2⌽s⍴⍵,⍺⍴⍨s←2⍴≢⍵}¨ ↑∘⎕a¨ 0 1 5 3
┌┬─┬─────────┬─────┐
││A│....A....│..A..│
││ │...B.B...│.B.B.│
││ │..C...C..│C...C│
││ │.D.....D.│.B.B.│
││ │E.......E│..A..│
││ │.D.....D.│     │
││ │..C...C..│     │
││ │...B.B...│     │
││ │....A....│     │
└┴─┴─────────┴─────┘

The last input expression is standalone (replacing dia in the penultimate input expression with its definition) and you can try it in a http://tryapl.org session.

Testing

Tests are often descriptive rather than prescriptive, asserting what a result should be rather than how to compute it. For that reason they tend to be more robust and easier to develop than the actual function. For example, it is much easier to test that r is a root of a polynomial than to compute it, or to check that S is a solution to an NP-complete problem than to derive it.

In this instance, we did not write the tests before writing the diamond functions, but we easily could have. The same expressive power in APL that enabled terse solutions of the diamond kata also enabled their concise and comprehensive testing.

testd d checks that d is a correct result of a diamond function, and returns a 1 if it is. Like the diamond functions, testd treats d in its totality, as an object (matrix) rather than as individual lines. This treatment allows straightforward statement of properties that would be more laborious otherwise. (For example, d≡⊖d and d≡⌽d assert that d is symmetric in the horizontal and vertical axes.)

assert←{⍺←'assertion failure' ⋄ 0∊⍵:⍺ ⎕signal 8 ⋄ shy←0}

testd←{                   ⍝ test the result of a diamond function
 assert (2=⍴⍴⍵)∧=/⍴⍵:     ⍝ square matrix
 assert (0=≢⍵)∨1=2|≢⍵:    ⍝ dimensions are 0 or odd
 assert ⍵≡⌽⍵:             ⍝ symmetric in vertical   axis
 assert ⍵≡⊖⍵:             ⍝ symmetric in horizontal axis
 n←⌈2÷⍨≢⍵                 ⍝ number of input "letters"
 q←(n,-n)↑⍵               ⍝ upper right quadrant
 assert (q=⍬⍴⍵)∨∘.=⍨⍳n:   ⍝ background or diagonal
 1                        ⍝ 1 iff everything OK
}

   testd '-' dia 'ABCD'
1

   testd 0 dia 3 1 4 2
1

If a check fails (or if there is an outright programming error), execution suspends at the first line with the error. While suspended an interactive debugging environment is provided for investigating the problem.

   testd 'ABCD'
assertion failure
testd[1] assert(2=⍴⍴⍵)∧=/⍴⍵: ⍝ square matrix
     ∧
   ⍵
ABCD
   ⍴⍵
4
   ⍴⍴⍵
1

testd is not a complete test, because it does not “know” what the specified background element or what the “letters” were. The operator T provides a complete test:

T←{                       ⍝ test a diamond function ⍺⍺
 d←⍺ ⍺⍺ ⍵                 ⍝ result of diamond function
 assert testd d:          ⍝ checks not knowing arguments
 n←≢⍵                     ⍝ number of input "letters"
 assert (⍴d)≡0 0⌈¯1+2×n:  ⍝ dimensions are 0 or ¯1+2×n
 assert (1≥n)∨⍺=⍬⍴d:      ⍝ background element
 assert ⍵≡1 1⍉(n,-n)↑d:   ⍝ "letters" in diagonal of quadrant
 1                        ⍝ 1 iff everything OK
}

⍺ dia T ⍵ checks that ⍺ dia ⍵ produces a correct result, and returns a 1 if it does (and suspends at the first line of error if it does not).

   '-' dia T 'ABCD'         '-' diam T 'ABCD'
1                        1

   0 dia T 3 1 4 2          0 diam T 3 1 4 1 5 9
1                        1

The following expression tests '-' dia x where x is a prefix of the alphabet ⎕a, from 0 to 26 letters:

   '-' (dia T)∘(⍴∘⎕a)¨ ¯1+3 9⍴⍳27
1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1

The following expression tests 0 dia n?n for n from 0 to 99. The “letters” are n?n, a random permutation of order n:

   0 (dia T)∘(?⍨)¨ 10 10⍴¯1+⍳100
1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1

References

  1. Rose, Seb, Recycling Tests in TDD, Claysnow Blog, 2014-11-23.
  2. Cockburn, Alistair, Thinking Before Programming, Alistair Cockburn Blog, 2014-11-29.
  3. Jeffries, Ron, TDD on the Diamond Problem, RonJeffries.com Blog, 2014-11-29.
  4. Rusk, John, An Experiment in Think-First Development, AgileWiki Blog, 2014-12-01.
  5. Hui, Roger, and Kenneth E. Iverson, J Introduction and Dictionary,1990-2014; \. entry
  6. Hui, Roger, Identity Matrix, Jwiki Essay, 2005-11-18.
  7. Scholes, John, Merge, Dfns, 2012-03-30.
  8. Legrand, Bernard, Mastering Dyalog APL, 2009-11.

NOTE: A full description and tutorial of the language can be found in [7].

Cholesky Decomposition

Morten was visiting Dyalog clients and forwarded a request: Can we have the Cholesky decomposition?

If A is a Hermitian, positive-definite matrix, its Cholesky decomposition [0] is a lower-triangular matrix L such that A ≡ L +.× +⍉L. The matrix L is a sort of “square root” of the matrix A.

For example:

   ⎕io←0 ⋄ ⎕rl←7*5 ⋄ ⎕pp←6

   A←t+.×⍉t←¯10+?5 5⍴20
   A
231   42  ¯63  16  26
 42  199 ¯127 ¯68  53
¯63 ¯127  245  66 ¯59
 16  ¯68   66 112 ¯75
 26   53  ¯59 ¯75  75

   L←Cholesky A
   L
15.1987   0        0        0       0
 2.7634  13.8334   0        0       0
¯4.1451  ¯8.35263 12.5719   0       0
 1.05272 ¯5.12592  2.1913   8.93392 0
 1.71067  3.48957 ¯1.81055 ¯6.15028 4.33502

   A ≡ L +.× +⍉L
1

For real matrices, “Hermitian” reduces to symmetric and the conjugate transpose +⍉ to transpose . The symmetry arises in solving least-squares problems.

Some writers asserted that an algorithm for the Cholesky decomposition “cannot be expressed without a loop” [1] and that “a Pascal program is a natural way of expressing the essentially iterative algorithm” [2]. You can judge for yourself whether the algorithm presented here belies these assertions.

The Algorithm [3]

A recursive solution for the Cholesky decomposition obtains by considering A as a 2-by-2 matrix of matrices. It is algorithmically interesting but not necessarily the best with respect to numerical stability.

Cholesky←{
 ⍝ Cholesky decomposition of a Hermitian positive-definite matrix
    1≥n←≢⍵:⍵*0.5
    p←⌈n÷2
    q←⌊n÷2
    X←(p,p)↑⍵ ⊣ Y←(p,-q)↑⍵ ⊣ Z←(-q,q)↑⍵
    L0←∇ X
    L1←∇ Z-T+.×Y ⊣ T←(+⍉Y)+.×⌹X
    ((p,n)↑L0)⍪(T+.×L0),L1
}

The recursive block matrix technique can be used for triangular matrix inversion [4], LU decomposition [5], and QR decomposition [6].

Proof of Correctness

The algorithm can be stated as a block matrix equation:

  ┌───┬───┐          ┌──────────────┬──────────────┐
  │ X │ Y │          │   L0 ← ∇ X   │       0      │
∇ ├───┼───┤  ←→  L ← ├──────────────┼──────────────┤ 
  │+⍉Y│ Z │          │    T+.×L0    │L1 ← ∇ Z-T+.×Y│
  └───┴───┘          └──────────────┴──────────────┘

where T←(+⍉Y)+.×⌹X. To verify that the result is correct, we need to show that A≡L+.×+⍉L and that L is lower triangular. For the first, we need to show:

┌───┬───┐     ┌──────┬───────┐     ┌────────┬────────┐
│ X │ Y │     │  L0  │   0   │     │  +⍉L0  │+⍉T+.×L0│
├───┼───┤  ≡  ├──────┼───────┤ +.× ├────────┼────────┤
│+⍉Y│ Z │     │T+.×L0│   L1  │     │    0   │  +⍉L1  │
└───┴───┘     └──────┴───────┘     └────────┴────────┘

that is:

(a)  X     ≡ L0 +.× +⍉L0
(b)  Y     ≡ L0 +.× +⍉ T+.×L0
(c)  (+⍉Y) ≡ (T+.×L0) +.× +⍉L0
(d)  Z     ≡ ((T+.×L0) +.× (+⍉T+.×L0)) + (L1+.×+⍉L1)

(a) holds because L0 is the Cholesky decomposition of X.

(b) is seen to be true as follows:
L0 +.× +⍉ T+.×L0
L0 +.× +⍉ ((+⍉Y)+.×⌹X)+.×L0 definition of T
L0 +.× (+⍉L0)+.×(+⍉⌹X)+.×Y +⍉A+.×B ←→ (+⍉B)+.×+⍉A and +⍉+⍉Y ←→ Y
(L0+.×+⍉L0)+.×(+⍉⌹X)+.×Y +.× is associative
X+.×(+⍉⌹X)+.×Y (a)
X+.×(⌹X)+.×Y X and hence ⌹X are Hermitian
I+.×Y associativity; matrix inverse
Y identity matrix

(c) follows from (b) by application of +⍉ to both sides of the equation.

(d) turns on that L1 is the Cholesky decomposition of Z-T+.×Y:

((T+.×L0)+.×(+⍉T+.×L0)) + (L1+.×+⍉L1)
((T+.×L0)+.×(+⍉T+.×L0)) + Z-T+.×Y
((T+.×L0)+.×(+⍉L0)+.×+⍉T) + Z-T+.×Y
(T+.×X+.×+⍉T) + Z-T+.×Y
(T+.×X+.×+⍉(+⍉Y)+.×⌹X) + Z-T+.×Y
(T+.×X+.×(+⍉⌹X)+.×Y) + Z-T+.×Y
(T+.×X+.×(⌹X)+.×Y) + Z-T+.×Y
(T+.×I+.×Y) + Z-T+.×Y
(T+.×Y) + Z-T+.×Y
Z

Finally, L is lower triangular if L0 and L1 are lower triangular, and they are by induction.

A Complex Example

   ⎕io←0 ⋄ ⎕rl←7*5

   A←t+.×+⍉t←(¯10+?5 5⍴20)+0j1ׯ10+?5 5⍴20
   A
382        17J131  ¯91J¯124 ¯43J0107  20J0035
 17J¯131  314     ¯107J0005 ¯60J¯154  26J¯137
¯91J0124 ¯107J¯05  379       49J0034  20J0137
¯43J¯107  ¯60J154   49J¯034 272       35J0103
 20J¯035   26J137   20J¯137  35J¯103 324

   L←Cholesky A

   A ≡ L +.× +⍉L
1
   0≠L
1 0 0 0 0
1 1 0 0 0
1 1 1 0 0
1 1 1 1 0
1 1 1 1 1

A Personal Note

This way of computing the Cholesky decomposition was one of the topics of [7] and was the connection (through Professor Shlomo Moran) by which I acquired an Erdős number of 2.

References

  1. Wikipedia, Cholesky decomposition, 2014-11-25.
  2. Thomson, Norman, J-ottings 7, The Education Vector, Volume 12, Number 2, 1995, pp. 21-25.
  3. Muller, Antje, Tineke van Woudenberg, and Alister Young, Two Numerical Algorithms in J, The Education Vector, Volume 12, Number 2, 1995, pp. 26-30.
  4. Hui, Roger, Cholesky Decomposition, J Wiki Essay, 2005-10-14.
  5. Hui, Roger, Triangular Matrix Inverse, J Wiki Essay, 2005-10-27.
  6. Hui, Roger, LU Decomposition, J Wiki Essay, 2005-10-31.
  7. Hui, Roger, QR Decomposition, J Wiki Essay, 2005-10-30.
  8. Ibarra, Oscar, Shlomo Moran, and Roger Hui, A Generalization of the Fast LUP Matrix Decomposition Algorithm and Applications, Journal of Algorithms 3, 1982, pp. 45-56.

Quicksort in APL

Quicksort is a classic sorting algorithm invented by C.A.R. Hoare in 1961 [0, 1]. It has been known for some time that quicksort has a terse rendition in APL [2]. To get right to it, here is the code:

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

The “pivot” ⍵⌷⍨?≢⍵ is randomly chosen. ((∇<S)⍪=S⍪(∇>S)) is a fork, selecting the parts of which are less than the pivot, equal to the pivot, and greater than the pivot. The function is recursively applied to the first and the last of these three parts.

      ⎕io←0 ⋄ ⎕rl←7*5

      ⎕←x←?13⍴20
3 2 19 16 11 4 18 17 9 17 7 3 1
      Q x
1 2 3 3 4 7 9 11 16 17 17 18 19

The variant Q1 obtains by enclosing each of the three parts. Its result exhibits 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 result of Q1 can be different on the same argument, as illustrated below:

      Q1←{1≥≢⍵:⍵ ⋄ S←{⍺⌿⍨⍺ ⍺⍺ ⍵} ⋄ ⍵((⊂∘∇<S)⍪(⊂=S)⍪(⊂∘∇>S))⍵⌷⍨?≢⍵}

      Q1 x
┌──────┬───┬────────────────────────────────────┐
│┌─┬─┬┐│3 3│┌─┬─┬──────────────────────────────┐│
││1│2│││   ││4│7│┌┬─┬─────────────────────────┐││
│└─┴─┴┘│   ││ │ │││9│┌──┬──┬─────────────────┐│││
│      │   ││ │ │││ ││11│16│┌─────────┬──┬──┐││││
│      │   ││ │ │││ ││  │  ││┌┬─────┬┐│18│19│││││
│      │   ││ │ │││ ││  │  ││││17 17│││  │  │││││
│      │   ││ │ │││ ││  │  ││└┴─────┴┘│  │  │││││
│      │   ││ │ │││ ││  │  │└─────────┴──┴──┘││││
│      │   ││ │ │││ │└──┴──┴─────────────────┘│││
│      │   ││ │ │└┴─┴─────────────────────────┘││
│      │   │└─┴─┴──────────────────────────────┘│
└──────┴───┴────────────────────────────────────┘

      Q1 x
┌───────────────────────┬─┬─────────────────────────┐
│┌──────────────────┬─┬┐│9│┌──┬──┬─────────────────┐│
││┌┬─┬─────────────┐│7│││ ││11│16│┌───────────┬──┬┐││
││││1│┌┬─┬────────┐││ │││ ││  │  ││┌┬─────┬──┐│19││││
││││ │││2│┌┬───┬─┐│││ │││ ││  │  ││││17 17│18││  ││││
││││ │││ │││3 3│4││││ │││ ││  │  ││└┴─────┴──┘│  ││││
││││ │││ │└┴───┴─┘│││ │││ ││  │  │└───────────┴──┴┘││
││││ │└┴─┴────────┘││ │││ │└──┴──┴─────────────────┘│
││└┴─┴─────────────┘│ │││ │                         │
│└──────────────────┴─┴┘│ │                         │
└───────────────────────┴─┴─────────────────────────┘

The enlist of the result of Q1 x is the same as Q x, the sort of x:

Q1 x
1 2 3 3 4 7 9 11 16 17 17 18 19
      Q x
1 2 3 3 4 7 9 11 16 17 17 18 19

This note is meant to explore the workings of a classical algorithm. To actually sort data in Dyalog, it is more convenient and more efficient to use {⍵[⍋⍵]}. Earlier versions of this text appeared in [3, 4].

References

  1. Hoare, C.A.R, Algorithm 63: Partition, Communications of the ACM, Volume 4, Number 7, 1961-07.
  2. Hoare, C.A.R, Algorithm 64: Quicksort, Communications of the ACM, Volume 4, Number 7, 1961-07.
  3. Hui, Roger K.W., and Kenneth E. Iverson, J Introduction and Dictionary, 1991-2014; if. entry.
  4. Hui, Roger K.W., Quicksort, J Wiki Essay, 2005-09-28.
  5. Hui, Roger K.W. Sixteen APL Amuse-Bouches, 2014-11-02.

Ken Iverson’s Favourite APL Expression?

What was Ken Iverson’s favourite APL expression? I don’t know that he had one and if he had I don’t know what it was, but if I have to guess …

From Sixteen APL Amuse-Bouches:

The expression (0,x)+(x,0) or its commute, which generates the next set of binomial coefficients, is present in the document that introduced APL\360 in 1967 [20, Fig.1] and the one that introduced J in 1990 [21, Gc&Gd]; in Elementary Functions: An Algorithmic Treatment in 1966 [22, p.69], in APL\360 User’s Manual in 1968 [23, A.5], in Algebra: An Algorithmic Treatment in 1972 [24, p.141], in Introducing APL to Teachers in 1972 [25, p.22], in An Introduction to APL for Scientists and Engineers in 1973 [26, p.19], in Elementary Analysis in 1976 [27, ex.1.68], in Programming Style in APL in 1978 [28, §6], in Notation as a Tool of Thought in 1980 [29, A.3], in A Dictionary of APL in 1987 [30, m∇n], and probably others.

The expression in action:

   ⎕←x←,1
1
   ⎕←x←(0,x)+(x,0)
1 1
   ⎕←x←(0,x)+(x,0)
1 2 1
   ⎕←x←(0,x)+(x,0)
1 3 3 1
   ⎕←x←(0,x)+(x,0)
1 4 6 4 1
   ⎕←x←(0,x)+(x,0)
1 5 10 10 5 1
   ⎕←x←(0,x)+(x,0)
1 6 15 20 15 6 1

It is easily seen from the expression that the n-th vector of binomial coefficients is a palindrome and that its sum is 2*n.