Solving the 2014 APL Problem Solving Competition – Cryptography Problem 2

This post is the continuation of the series where we examine some of the problems selected for the 2014 APL Problem Solving Competition. In this post we’ll continue looking at the cryptography problems from Phase II that we started looking at in a previous blog post.

Cryptography Problem 2 – Book Cipher Variation

Task 1 – Let’s Get Normal

The first task is to normalise some text by weeding out non-alphabetic characters, collapsing consecutive spaces and converting to upper case. It’s possible to do this all by hand with APL, but it’s much easier to use the ⎕R operator to search for and replace regular expressions.

Taking the transformations in turn, here’s how to convert non-alphabetic characters in message to spaces:

('[^[:alpha:]]'⎕R' ')⍵

Here’s how to convert multiple consecutive spaces to a single space:

(' +'⎕R' ')⍵

And here’s how to convert every alphabetic character to upper case:

('.'⎕R'\u&')⍵

We can combine the first two of these, by converting any sequence of one or more non-alphabetic characters to a single space, giving the following implementation:

Normalise←{
    text←⎕SE.UnicodeFile.ReadText ⍵
    ('[^[:alpha:]]+' '.'⎕R' ' '\u&'⍠'Mode' 'D')text
}

The option 'Mode' 'D' tells ⎕R to operate in Document mode, which processes the whole file at once instead of line by line, as we are not interested in preserving the original line breaks. Here it is in action:

      70↑bor←Normalise'/home/jay/Desktop/BillOfRights.txt'
THE PREAMBLE TO THE BILL OF RIGHTS CONGRESS OF THE UNITED STATES BEGUN

Task 2 – Encryption

In this cipher there are lots of different ways of encoding each character of the message, and we are free to pick any of them. In order to try to “minimise the number of duplicated pairs in the result”, we simply pick randomly whenever we have a free choice. The function pickone helps with this. Given a boolean vector, it first uses {⍵/⍳⍴⍵} to get a vector of the indices of all the 1 bits, and then uses {⍵[?≢⍵]} to choose one of these indices at random.

In this coding of BookEncrypt, the anonymous inner dfn encodes a single character of the message into a (word offset) pair. These pairs are joined together with ⊃,/, a common pattern for catenating strings. The Disclose is required because, in Dyalog, reduction always reduces the rank of its argument, so ,/ on a vector of strings returns a scalar: the enclose of the catenated strings.

BookEncrypt←{
    pickone←{⍵[?≢⍵]}∘{⍵/⍳⍴⍵}
    b←' ',⍺                     ⍝ b has a space wherever a word starts in the key.
    s←{⍵/⍳⍴⍵}b=' '              ⍝ Get the indices of all the word starts.
    ⊃,/{
         p←pickone b=⍵          ⍝ Choose a random occurrence of letter ⍵
         p-←s                   ⍝ and get its offset within each word.
         w←pickone{(0<⍵)∧⍵≤20}p ⍝ Choose a word with a reasonable offset
         w(p[w])                ⍝ and return the (word offset) pair.
    }¨⍵                         ⍝ ... for each letter in the message.
}

Here it is in action:

      ⊢cipher←bor BookEncrypt 'MYSECRETMESSAGE'
480 11 523 11 440 6 115 5 78 16 579 18 696 20 330 16 544 4 658 17 400 9 661 11
      246 18 186 4 482 13

Task 3 – Decryption

Decryption is simpler then encryption, because there is no need to make random choices. All we have to do is:

  • Find the index of the start of each word in the key, as before.
  • Split the input into pairs of numbers.
  • For each pair, find the character in the key at the specified offset from the start of the specified word.

There are various ways to split the input into pairs of numbers. Here, we do it with the Rank operator (). Encrypting an N-character message gives a vector of 2×N numbers. To split it into pairs we first reshape it into a matrix with N rows and 2 columns; and then use f⍤1 to apply f to the rank-1 subarrays of this matrix, which are its row vectors.

Here’s the code:

BookDecrypt←{
    b←' ',¯1↓⍺                  ⍝ b has a space wherever a word starts in the key.
    s←{⍵/⍳⍴⍵}b=' '              ⍝ Get the indices of all the word starts.
    {
        (w o)←⍵                 ⍝ Get word number and offset
        b[s[w]+o]               ⍝ and find the character at that position
    }⍤1⊢(0.5×≢⍵)2⍴⍵             ⍝ ... for each pair of numbers in the input.
}

And here it is in action:

      bor BookDecrypt cipher
MYSECRETMESSAGE

To be continued…

Solving the 2014 APL Problem Solving Competition – Cryptography Problem 1

This post is a continuation of the series where we examine some of the problems selected for the 2014 APL Problem Solving Competition. I’ll start by looking at the cryptography problems from Phase II.

Cryptography Problem 1 – Vigenère Cipher

The cipher is described using a large table of letters, but you don’t need to create this table in APL. Instead, you can convert each letter of the key and the corresponding letter of the message to numbers; add them together; and then convert the result back to a letter. Some subtleties:

  • We want the result of the addition to “wrap around” at the end of the alphabet, so 25 + 1 = 26 (Z) but 25 + 2 = 27 which should wrap round to 1 (A). This is called modular arithmetic.
  • We can implement modular arithmetic easily using APL’s Residue primitive, but this works most naturally when the alphabet is numbered from 0 to 25, rather than from 1 to 26. To use 0-based numbering, I’ll set ⎕IO←0 in the code below.
  • An idiomatic way of converting letters to numbers is to look them up in the ⎕A, the upper case alphabet: ⎕A⍳⍵. To convert the other way we can simply index into ⎕A: ⎕A[⍵]. (For more comprehensive conversions between characters and numbers see the ⎕UCS system function.)

The code is then straightforward:

VigEncrypt←{
    ⎕IO←0                       ⍝ use 0-based numbering
    key←(⍴⍵)⍴⎕A⍳⍺               ⍝ convert key to numbers and expand to length of message
    ⎕A[26|(⎕A⍳⍵)+key]           ⍝ add key to message and convert back to characters
}

To encrypt the message we added the key (modulo 26). To decrypt, we need to subtract the key (modulo 26) instead, so we can go from VigEncrypt to VigDecrypt just by changing one character in the code. Note that the result of the subtraction might be negative, but Residue (unlike the remainder or modulo operations in some other programming languages you may have used) will still give us the correct result in the range 0 to 25. For example, 26|¯3 is 23.

VigDecrypt←{
    ⎕IO←0                       ⍝ use 0-based numbering
    key←(⍴⍵)⍴⎕A⍳⍺               ⍝ convert key to numbers and expand to length of message
    ⎕A[26|(⎕A⍳⍵)-key]           ⍝ subtract key from message and convert back to characters
}

Here are the functions in action:

      key←'KEIVERSON'
      key VigEncrypt'APLISFUNTOUSE'
KTTDWWMBGYYAZ
      key VigDecrypt key VigEncrypt'APLISFUNTOUSE'
APLISFUNTOUSE

To be continued…

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.

Three-and-a-bit

The most obvious expression for computing π in APL is ○1. But what if you can’t remember how works, or your O key is broken, or you feel like taking the road less travelled? With thanks to Wikipedia’s excellent list of Approximations of π, here are some short sweet APL expressions for three-and-a-bit:

      3                                 ⍝ very short
3
      4                                 ⍝ not so sweet
4
      s←*∘0.5                           ⍝ let's allow ourselves some square roots
      +/s 2 3
3.14626436994197234232913506571557
      31*÷3
3.141380652391393004493075896462748
      +/1.8*1 0.5                       ⍝ Ramanujan
3.141640786499873817845504201238766
      s 7+s 6+s 5
3.141632544503617840472137945142766
      ÷/7 4*7 9
3.141567230224609375
      9⍟995
3.141573605337628094187009177086444
      355÷113
3.141592920353982300884955752212389
      s s 2143÷22                       ⍝ Ramanujan again
3.141592652582646125206037179644022
      +∘÷/3 7 15 1 292                  ⍝ continued fraction
3.14159265301190260407226149477373
      ÷/63 25×17 7+15×s 5
3.14159265380568820189839000630151
      (1E100÷11222.11122)*÷193
3.141592653643822210363178893440074
      (⍟744+640320*3)÷s 163             ⍝ Ramanujan yet again
3.141592653589793238462643383279727

This last one is accurate to more places than I ever learned in my youth!

Technical note: to get plenty of precision, these examples were evaluated with 128-bit decimal floating-point, by setting ⎕FR←1287 and ⎕PP←34.

For more on continued fractions, see cfract in the dfns workspace.