Permutations

I started composing a set of APL exercises in response to a query from a new APL enthusiast who attended Morten’s presentation at Functional Conf, Bangalore, October 2014. The first set of exercise are at a low level of difficulty, followed by another set at an intermediate level. One of the intermediate exercises is:

Permutations

a. Compute all the permutations of ⍳⍵ in lexicographic order. For example:

   perm 3
0 1 2
0 2 1
1 0 2
1 2 0
2 0 1
2 1 0

b. Write a function that checks whether is a solution to perm ⍺, without computing perm ⍺. You can use the function assert. For example:

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

pcheck←{
  assert 2=⍴⍴⍵:
  assert (⍴⍵)≡(!⍺),⍺:
  …
  1
}

   6 pcheck perm 6
1

   4 pcheck 2 4⍴0 1 2 3, 0 1 3 2
assertion failure
pcheck[2] assert(⍴⍵)≡(!⍺),⍺:
         ∧

c. What is the index of permutation in perm ⍺? Do this without computing all the permutations. For example:

   7 ip 1 6 5 2 0 3 4    ⍝ index from permutation
1422

(The left argument in this case is redundant, being the same as ≢⍵.)

d. What is the -th permutation of ⍳⍺? Do this without computing all the permutations. For example:

   7 pi 1442             ⍝ permutation from index
1 6 5 2 0 3 4

   (perm 7) ≡ 7 pi⍤0 ⍳!7
1

The Anagram Kata

Coincidentally, Gianfranco Alongi was attempting in APL the anagrams kata from Cyber Dojo:

Write a program to generate all potential anagrams of an input string.

For example, the potential anagrams of “biro” are
biro bior brio broi boir bori
ibro ibor irbo irob iobr iorb
rbio rboi ribo riob roib robi
obir obri oibr oirb orbi orib

This is essentially the same program/exercise/kata, because the potential anagrams are 'biro'[perm 4]. You can compare solutions in other languages to what’s here (google “anagrams kata”).

Spoiler Alert

Deriving a Solution

I am now going to present solutions to part a of the exercise, generating all permutations of ⍳⍵.

Commonly, in TDD (test-driven development) you start with a very simple case and try to extend it successively to 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 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 entry 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 advancements 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 the shortest known solution, 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

    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 functions be written to work in either index-origin and therefore have ⎕io sprinkled hither, thither, and yon.

1987, Some Uses of { and }

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

Written in Dictionary APL, wherein: ⍪⌿⍵ ←→ ⊃⍪⌿⊂⍤¯1⊢⍵ and differs from its definition in Dyalog APL; is equivalent to in 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 and in a Vector article. The lessons translate directly into Dyalog APL.

2008, http://dfns.dyalog.com/n_pmat.htm

   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 a or 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] :-)​​

In Praise of Magic Functions: Part II

Part I of this post was concerned with the development speed and execution speed of magic functions and should be read before this post.

Benefitting from Future and Past Improvements

Looking at the magic function for key in Part I, we see that its performance depends on the following APL primitives, listed with information on when they were last improved (or when they will be improved):

expression factor version
⍋i 2-18 12.1
⎕DR 1-3.8 14.1
⍳⍨x 2-4 14.1
⍳⍨i 3-6 14.1
↑x 5-12 15.0
⊂[a]x 1-1.75 15.0
b⊂[⎕IO]x 1.3-42 14.1

Key was introduced in version 14.0 and x{⍺,f⌿⍵}⌸y is currently computed by the general case. Comparing the speed in version 14.0 to 15.0 as of June 2015:

   x←?1e6⍴5e4
   y←?1e6⍴0

   1 1 cmpx 'x{⍺,+/⍵}⌸y'    ⍝ 14.0
0.1735

   1 1 cmpx 'x{⍺,+/⍵}⌸y'    ⍝ 15.0
0.12855

   0.1735 ÷ 0.12855
1.34967

a factor of 1.3 improvement without touching the key code. Special code had been planned for {⍺,f/⍵}⌸, and we can get an idea of the amount of speed up:

   cmpx 'x{⍺,+/⍵}⌸y' '(∪x),⍤0⊢x{+/⍵}⌸y'    ⍝ 15.0
 x{⍺,+/⍵}⌸y       → 1.25E¯1 |   0% ⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕
 (∪x),⍤0⊢x{+/⍵}⌸y → 1.16E¯2 | -91% ⎕⎕⎕

Perhaps the special code for {⍺,f/⍵}⌸ can be a magic function?

Magic functions also benefit from past improvements. The primitive ⍕b (in the Suggestivity subsection below) has a magic function implementation: MAGIC("(¯1+2×≢⍵)⍴(2 2⍴'0 1 ')[⍵;]⊣⎕io←0"). Then I thought, why not get an extra level of performance by doing it in native C? Benchmarks done after that work showed that the native C implementation was faster by a factor of 25 but the magic function was faster by a factor of 40. What happened? On reflection, the speed of the code depends crucially on the speed of x[b;…;] and it turned out that x[b;…;] was previously sped up (version 13.2). Apparently the x[b;…;] code is faster than a casually-written C program.

Finally, I expect magic functions are prime candidates for being fed to the APL compiler when it becomes possible to do so.

APL as a Tool of Implementation

Being APL code, magic functions have the characteristic terseness of APL and all that that entails. K.E. Iverson’s 1979 Turing Lecture Notation as a Tool of Thought lists five important characteristics of notation. We examine how they play out in this context.

Ease of Expression
Notation as a Tool of Thought entitles this characteristic as “ease of expressing constructs arising in problems”, explicated as follows (§1.1):

If it is to be effective as a tool of thought, a notation must allow convenient expression not only of notions arising directly from a problem, but also of those arising in subsequent analysis, generalization, and specialization.

For magic functions, the analogous concept would be facility in program development, debugging, tracing, profiling, benchmarking, and other activities required in implementation. Dyalog APL has these in spades.

Suggestivity
It was first identified that ⍕b can be sped up by doing (¯1+2×≢⍵)⍴(2 2⍴'0 1 ')[⍵;]⊣⎕io←0, the idea being that if f is an expensive function and x is from a small domain D, then f x can be computed by (f D)[⍵;]. (A similar idea is described in §4.3 of Notation as a Tool of Thought.) For the idea to be effective, f does not have to be very expensive, just more expensive than indexing, and x is sufficiently large. A few speed-ups in the pattern have been identified:

⍕b
x⍕b
float ⍳ b
float ⍳ int8
scalar×b and b×scalar
scalar*b

All of these can be implemented by magic functions. As well, the utility of the pattern (f D)[i;] motivates extra attention on x[i;…;]. Fortuitously, x[b;…;] has previously been made fast.

Subordination of Detail
Writing in APL means not having to deal with many details involved in writing in C in the interpreter, including:

  • workspace compaction
  • arrays changing their internal datatype
  • call by name
  • allocating and releasing memory
  • index errors AKA memory read and write exceptions
  • different numeric types
  • different character types
  • loops and nested loops
  • declarations
  • etc.

Economy
Economy refers to the size of the vocabulary and the complexity of the grammar. Here, the comparison is not just between APL and C (in which comparison APL wins), but between APL and the collection of programs, utilities, macros, data structures, typedefs, calling conventions, …, accumulated by the C source code over the years. These have not been as rigorously defined and executed as APL.

Amenability to Formal Manipulation
The opening paragraph of the present text had a different magic function for ∧.= in a previous version. In repeated readings of the MSS and in meditation on the ideas I realized that the line of code can be made simpler:

MAGIC("((≢⍺)↑i)∘.=(≢⍺)↓i←⍺⍳⍺⍪⍉⍵");    ⍝ old version
MAGIC("(≢⍺)(↑∘.=↓)⍺⍳⍺⍪⍉⍵");           ⍝ current version

I daresay that such simplification would be harder to conceive with a more verbose statement of the algorithm.

In Praise of Magic Functions: Part I

A magic function is an APL-coded (dfn) computation in the C source of the interpreter. For example, some cases of ∧.= are coded as MAGIC("(≢⍺)(↑∘.=↓)⍺⍳⍺⍪⍉⍵"). The rubric “magic function” includes magic functions and magic operators.

Acknowledgments. Magic functions in Dyalog are made possible due to work done by John Scholes and Richard Smith.

Development Speed

A magic function implementation takes an order of magnitude less time to write than a C-coded implementation. A case in point is ∘.≡ (and ∘.≢) on enclosed character vectors, as recounted in the blog post A Speed-Up Story. The chronology is recovered from the G-mail, Mantis and SVN systems.

2014-10-27 04:24 Nicolas initial e-mail describing the problem
2014-10-27 07:33 Roger proposes i∘.=i←a⍳a as a solution
2014-10-27 07:36 Jay objects that solution should check ⎕CT
2014-10-27 07:40 Roger responds to objection
2014-10-27 10:29 Roger creates Mantis issue
2014-10-27 13:57 Roger SVN commit
2014-10-27 18:12 Roger reports factor 6-64 speed-up and submits blog post
2014-10-28 16:33 Roger SVN commit to fix a bug

After checking that ⎕CT is not required, the main processing in the C source is as follows:

  • 1 iff a “selfie”, i.e. the left and right arguments are equal as C pointers
  • eq 1 if ∘.≡ ; 0 if ∘.≢
  • 1 iff the right argument has rank 1
#define CASE(x,y,z)  ((x<<2)+(y<<1)+(z))

switch(CASE(c,eq,r))
{
  case CASE(0,0,0): MAGIC("((,⍺)⍳⍺)∘.≠((,⍺)⍳⍵)"); break;
  case CASE(0,0,1): MAGIC("(  ⍺ ⍳⍺)∘.≠(  ⍺ ⍳⍵)"); break;
  case CASE(0,1,0): MAGIC("((,⍺)⍳⍺)∘.=((,⍺)⍳⍵)"); break;
  case CASE(0,1,1): MAGIC("(  ⍺ ⍳⍺)∘.=(  ⍺ ⍳⍵)"); break;
  case CASE(1,0,0): MAGIC("∘.≠⍨(,⍵)⍳⍵");          break;
  case CASE(1,0,1): MAGIC("∘.≠⍨  ⍵ ⍳⍵");          break;
  case CASE(1,1,0): MAGIC("∘.=⍨(,⍵)⍳⍵");          break;
  case CASE(1,1,1): MAGIC("∘.=⍨  ⍵ ⍳⍵");          break;
}

Execution Speed

Development speed for a magic function need not be at the expense of execution speed. As indicated above, ∘.≡ is 6 to 64 times faster than the previous (C coded) implementation. Factors for a few other magic function implementations:

expression factor version
x∘.≡y and x∘.≢y 6-64 14.1
x∧.=y and x∨.≠y 1.5-4 14.1
,/y 1-∞ 15.0
x⊥y 1-26 15.0
x⊤y 1-3.3 15.0

One can always make a magic function faster by hand-translating it from APL into C, and in so doing save on the tokenizing (scanning) and parsing costs. However, the effort increases the development time and the code size (see Special Cases below) for (I argue) not much reduction in execution time. I also expect that such translation can be done automatically by the APL compiler in the future.

Accurate estimates for the amount of speed up obtain readily from APL benchmarks. For example, x⍕b where x is a scalar or 2-element vector and b is a Boolean array, is a candidate for a magic function implementation, and:

   f←{((¯1↓s),(⊃⌽⍴t)×⊃⌽s←⍴⍵)⍴(t←⍺⍕⍪0 1)[⍵;]}

   b←1=?10⍴2     ⍝ small argument
   cmpx '2 f b' '2⍕b'
 2 f b → 5.83E¯6 |    0% ⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕
 2⍕b   → 1.23E¯5 | +110% ⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕

   b←1=?1e6⍴2    ⍝ large argument
   cmpx '8 2 f b' '8 2⍕b'
 8 2 f b → 9.75E¯3 |     0% ⎕
 8 2⍕b   → 3.35E¯1 | +3339% ⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕

We can say with confidence that x⍕b can be made 2 to 34 times faster before actually proceeding with the implementation.

Magic functions can be used when execution speed is not paramount. For example a problem was reported with !∘y⍣¯1:

   !∘25⍣¯1 ⊢ 1e4
DOMAIN ERROR
   !∘25⍣¯1⊢10000
  ∧

The problem can be solved readily with an APL function using Newton iteration, with a relatively complicated computation for the initial estimate:

   f←{-∘⍵∘⍺⍺{⍵-(⊃t)ׯ0.01÷-/t←⍺⍺ ⍵×1 1.01}⍣≡⊃⍋|⍵-(⍳!⊢)⌊0.5+⍺⍺ 1}
   ⎕←x←!∘25 f 1e4
3.85142
   x!25
10000

General Case vs. Special Cases

Magic functions are well-suited for implementing operators. An operator has potentially tens or hundreds of cases, and it would be onerous, not cost-effective and unnecessary to write code for each case. Instead, the general case can be implemented quickly and succinctly by a magic function, and the saved effort can be devoted to cases deemed worthy of attention. The rank and key operators are done this way. The magic function for the general case of key:

MAGIC(
  "0=⎕nc'⍺':⍵ ∇ ⍳≢⍵    ⋄"  // monadic case: f⌸ x ←→ x f⌸ ⍳≢x
  "⍺ ⍺⍺{⍺ ⍺⍺ ⍵}{        "
  "  ⎕ml←1             ⋄"  // needed by ↑⍵ and ⍺⊂⍵ 
  "  j←⍋i⊣⎕dr i←⍳⍨⍺    ⋄"
  "  ↑ (⊂[1↓⍳⍴⍴⍺]((⍳≢i)=⍳⍨i)⌿⍺) ⍺⍺¨ (2≠/¯1,i[j])⊂[⎕io]⍵⌷⍨⊂j"
  "}⍵                   "
);

Before execution reaches the general case, special cases are detected and are implemented by special code, usually in C. These special cases are:

operand version comment
{f⌿⍵} and ⊢∘(f⌿) 14.0 for f one of + ⌈ ⌊ or (for Boolean right arguments) one of ∧ ∨ = ≠; also / instead of for vector right arguments
{⍺(f⌿⍵)} 15.0 for f as above
{⍺,f⌿⍵} 15.0 for f as above and numeric left arguments
{≢⍵} and ⊢∘≢ 14.0
{⍺(≢⍵)} 14.1
{⍺,≢⍵} 14.1 for numeric left arguments
{≢∪⍵} 14.1
{⍺(≢∪⍵)} 14.1
{⍺,≢∪⍵} 14.1 for numeric left arguments
{⊂⍵} and ⊢∘⊂ 14.0
{⍺⍵} 14.1
{⍺} and 14.1 ⊣⌸x ←→ ∪x if were extended like

Additional special cases can be implemented as required.

Special Cases

Since magic functions are so terse, sometimes it is worthwhile to make special cases which would not be made special cases if the implementation were more verbose and/or require more effort. The implementation of ∘.≡ (in the Development Speed section above) illustrates the point. In general, x∘.≡y ←→ ((,⍺)⍳⍺)∘.=((,⍺)⍳⍵). However, if x is a vector, the two ravels can be elided; moreover, if x and y are equal as C pointers, the two uses of can be reduced to one use (not only that, but to a “selfie” ⍵⍳⍵ if a vector). So instead of the one case:

MAGIC("((,⍺)⍳⍺)∘.=((,⍺)⍳⍵)");

we have

switch(CASE(c,eq,r))
{
    …
  case CASE(0,1,0): MAGIC("((,⍺)⍳⍺)∘.=((,⍺)⍳⍵)"); break;
  case CASE(0,1,1): MAGIC("(  ⍺ ⍳⍺)∘.=(  ⍺ ⍳⍵)"); break;
    …
  case CASE(1,1,0): MAGIC("∘.=⍨(,⍵)⍳⍵");          break;
  case CASE(1,1,1): MAGIC("∘.=⍨  ⍵ ⍳⍵");          break;
}

The extra cases are not too burdensome, and their detection is done in scalar code at which C is better than APL. The following benchmarks illustrate the difference that the special cases make:

   t ←' zero one two three four five six seven eight nine'
   t,←' zéro un deux trois quatre cinq six sept huit neuf'
   t,←' zero eins zwei drei vier fünf sechs sieben acht neun'
   u←1↓¨(' '=t)⊂t

   x←u[?300⍴≢u]    ⍝ vector selfie vs. non-selfie
   y←(⍴x)⍴x
   cmpx 'x∘.≡x' 'x∘.≡y'
 x∘.≡x → 1.48E¯4 |   0% ⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕
 x∘.≡y → 1.93E¯4 | +30% ⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕

   x1←(1,⍴x)⍴x     ⍝ matrix selfie vs. non-selfie
   y1←(⍴x1)⍴x1
   cmpx 'x1∘.≡x1' 'x1∘.≡y1'
 x1∘.≡x1 → 1.50E¯4 |   0% ⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕
 x1∘.≡y1 → 1.97E¯4 | +31% ⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕

To be continued…Part II will describe the benefits from future improvements and APL as a tool of implementation.

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]

Solving the 2014 APL Problem Solving Competition – Cryptography Problem 3

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 conclude looking at the cryptography problems from Phase II that we started looking at in a previous blog post and continued in a further blog post.

Cryptography Problem 3 – Playfair Cipher

Task 1 – Squaring Off

The first task is to convert a string into a 5×5 Playfair table. The solution makes straightforward use of APL primitives:

  • Unique () to remove duplicate characters from the string.
  • Without (~) to find the rest of the alphabetic characters that are not mentioned in the string, and again to remove the character J.
PlayfairTable←{
    k←∪⍵
    5 5⍴(k,⎕A~k)~'J'
}

Here it is in action:

      ⊢table←PlayfairTable'KENNETHEIVERSON'
KENTH
IVRSO
ABCDF
GLMPQ
UWXYZ

Task 2 – Encryption

To encrypt a message we need to take two characters at a time, find their coordinates in the 5×5 Playfair table, swap their column coordinates and look up the letters at the new coordinates. There are lots of tricks in the code below of which we’ll describe just a few:

  • To process the message two characters at a time, appending to the result as we go, we use a tail-recursive dfn whose left argument accumulates the result. For an introduction to this technique see this implementation of the Fibonacci function.
  • To find letters in the Playfair table we first look them up in the ravel (i.e. the linearised form) of the table, and then use Encode () to convert this linear index into a pair of coordinates. Decode () does the opposite, converting a pair of coordinates back into a linear index.
  • Mix () and Split () are used to convert between two different representations of the coordinates of the two characters we’re encoding: either as a flat 2×2 matrix, or as a nested 2-vector of 2-vectors. The choice of representation is largely a matter of taste, and it might be fun to play with this part of the code. You could tweak it to work entirely with the flat representation, or entirely with the nested representation, rather than converting back and forth between them.
PlayfairEncrypt←{
    ⎕IO←0                       ⍝ To aid arithmetic modulo 5.
    t←⍺                         ⍝ The Playfair table.
    m←⍵,(2|≢⍵)/'Z'              ⍝ Add a Z if the message length is odd.
    ((m='J')/m)←'I'             ⍝ Convert J to I in the message.
    ''{                         ⍝ Start with an empty accumulator.
        0=≢⍵:⍺                  ⍝ Finished: return the accumulated encrypt.
        1=≢⍵:⍺ ∇ ⍵,'X'          ⍝ Only one character left: add an X.
        p←2↑⍵
        =/p:⍺ ∇(⊃⍵),'X',1↓⍵     ⍝ Duplicate character: insert an X.
        c←(⍴t)⊤(,t)⍳p           ⍝ Coords of each letter in the table.
        c←{
            =/1↑⍵:↑(⍴t)|0 1+↓⍵  ⍝ Same row: move to the right.
            =/1↓⍵:↑(⍴t)|1 0+↓⍵  ⍝ Same column: move down.
            0 1⌽⍵               ⍝ Else swap column coords.
        }c
        p←(,t)[(⍴t)⊥c]          ⍝ Look up new coords in table.
        (⍺,p)∇ 2↓⍵              ⍝ Tail recursion.
    }m
}

Here it is in action:

      ⊢cipher←table PlayfairEncrypt'HELLOWORLD'
KNMWQVZVVMCY

(Note that this result is slightly different from that given in the original problem description, because of some confusion about the rules for handling duplicate letters and odd message lengths, which were clarified later in the student competition forums.)

Task 3 – Decryption

Decryption is very similar to encryption. The differences are:

  • Letters found in the same row of the table need to move left instead of right, and letters in the same column need to move up instead of down.
  • We don’t need to worry about the input having an odd length, or containing the letter J.

Hence the code for PlayfairDecrypt is shorter than but very similar to PlayfairEncrypt:

PlayfairDecrypt←{
    ⎕IO←0                       ⍝ To aid arithmetic modulo 5.
    t←⍺                         ⍝ The Playfair table.
    ''{                         ⍝ Start with an empty accumulator.
        0=≢⍵:⍺                  ⍝ Finished: return the accumulated encrypt.
        p←2↑⍵
        c←(⍴t)⊤(,t)⍳p           ⍝ Coords of each letter in the table.
        c←{
            =/1↑⍵:↑(⍴t)|0 ¯1+↓⍵ ⍝ Same row: move to the left.
            =/1↓⍵:↑(⍴t)|¯1 0+↓⍵ ⍝ Same column: move up.
            0 1⌽⍵               ⍝ Else swap column coords.
        }c
        p←(,t)[(⍴t)⊥c]          ⍝ Look up new coords in table.
        (⍺,p)∇ 2↓⍵              ⍝ Tail recursion.
    }⍵
}

Here it is in action:

      table PlayfairDecrypt cipher
HELXLOWORLDX

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…