Solving the 2014 APL Problem Solving Competition – it’s as easy as 1 1 2 3…

Competition LogoThe winners of the 2014 APL Problem Solving Competition were recognized at the Dyalog ’14 user meeting and had a chance to share their competition experiences. Emil Bremer Orloff won the student competition and received $2500 USD and an expenses-paid trip to Dyalog ’14, while Iryna Pashenkovska took first place among the non-student entries and received a complimentary registration to Dyalog ’14. Our congratulations to them for their fine efforts! This post is the first of a series where we will examine some of the problems selected for this year’s competition. The problem we’ll discuss first is Problem 3 from Phase I dealing with the Fibonacci sequence.

Problem 3 – Tell a Fib Write a dfn that takes an integer right argument and returns that number of terms in the Fibonacci sequence. Test cases:

      {your_solution} 10
 1 1 2 3 5 8 13 21 34 55
      {your_solution} 1
 1
      {your_solution} 0   ⍝ should return an empty vector

Essentially, you start with the sequence 1 1 and concatenate the sum of the last 2 numbers to the end and repeat until the sequence has the correct number of terms. In Python, a solution might look something like this:

def fib(n): # return the first n elements of the Fibonacci sequence
    result = []
    a, b = 0, 1
    while 0 < n:
        result.append(b)
        a, b = b, a+b
        n -= 1
    return result

You can write nearly the same code in APL:

r←fibLooping n;a;b
  r←⍬
  (a b)←0 1
:While 0<n
  r,←b
  (a b)←b,a+b
  n-←1
:EndWhile

While it’s possible to write APL code that looks like Python (or many other languages), one of the judging criteria for the competition is the effective use of APL syntax – in this case, by avoiding the explicit loop. Two ways to do this are 1) using recursion and 2) using the power operator ().

Recursive Solution

      fibRecursive←{⍵<3:⍵⍴1 ⋄ {⍵,+/¯2↑⍵}∇⍵-1}

The neat thing about recursion is that the function keeps calling itself with a “simpler” argument until it reaches a base case and then passes the results back up the stack. Here the base case occurs when the argument () is less than 3 and the function returns ⍵⍴1 – it’s either an empty vector, 1, or 1 1 for values of of 0, 1 and 2 respectively. It’s easy to illustrate the recursive process by adding some code (⎕←) to display information at each level of recursion.

      fibRecursive←{⍵<3:⎕←⍵⍴1 ⋄ {⎕←⍵,+/¯2↑⍵}∇⎕←⍵-1}
      fibRecursive 10
9
8
7
6
5
4
3
2
1 1
1 1 2
1 1 2 3
1 1 2 3 5
1 1 2 3 5 8
1 1 2 3 5 8 13
1 1 2 3 5 8 13 21
1 1 2 3 5 8 13 21 34
1 1 2 3 5 8 13 21 34 55

The recursive solution above is termed “stack recursive” in that it creates a new level on the function call stack for each recursive call. We can modify the code slightly to implement a “tail recursive” solution which Dyalog can detect and optimize by avoiding creating those additional function call stack levels. You can see the effect of this as each level computes its result immediately:

      fibTail←{⍺←0 1 ⋄ ⍵=0:⎕←1↓¯1↓⍺ ⋄ (⎕←⍺,+/¯2↑⍺)∇ ⎕←⍵-1}
      fibTail 10
        fibTail 10
9
0 1 1
8
0 1 1 2
7
0 1 1 2 3
6
0 1 1 2 3 5
5
0 1 1 2 3 5 8
4
0 1 1 2 3 5 8 13
3
0 1 1 2 3 5 8 13 21
2
0 1 1 2 3 5 8 13 21 34
1
0 1 1 2 3 5 8 13 21 34 55
0
0 1 1 2 3 5 8 13 21 34 55 89
1 1 2 3 5 8 13 21 34 55

Power Operator Solution

      fibPower←{⍵⍴({⍵,+/¯2↑⍵}⍣(0⌈⍵-1))1}

The power operator is defined as {R}←{X}(f⍣g)Y. When the right operand g is a numeric integer scalar – in this case (0⌈⍵-1), the power operator applies its left operand function f{⍵,+/¯2↑⍵} cumulatively g times to its argument Y, which in this case is 1. Similar to the recursive solution, we can add ⎕← to see what happens at each application:

      fibPower←{⍵⍴({⎕←⍵,+/¯2↑⍵}⍣(0⌈⍵-1))1}
      fibPower 10
1 1
1 1 2
1 1 2 3
1 1 2 3 5
1 1 2 3 5 8
1 1 2 3 5 8 13
1 1 2 3 5 8 13 21
1 1 2 3 5 8 13 21 34
1 1 2 3 5 8 13 21 34 55
1 1 2 3 5 8 13 21 34 55

In discussing this blog post with my fellow Dyalogers, Roger Hui showed me a rather neat construct:

      fibRoger←{1∧+∘÷\⍵⍴1}     ⍝ Roger's original version
      fibBrian←{1∧÷(+∘÷)\⍵⍴1}  ⍝ tweaked to make result match contest examples
      fibRoger 10
1 2 3 5 8 13 21 34 55 89
      fibBrian 10
1 1 2 3 5 8 13 21 34 55

I’ve added the parentheses around +∘÷ to make it a bit clearer what the left operand to the scan operator \ is. Let’s break the code down…

      ⍵⍴1 ⍝ produces a vector of ⍵ 1's 
1 1 1 1 1 ⍝ for ⍵=5

Then we apply scan with the function +∘÷, which in this case has the effect of adding 1 to the reciprocal of the previous result:

      (+∘÷)\1 1 1 1 1
1 2 1.5 1.666666667 1.6

Note that the denominator of each term is the corresponding element of the Fibonacci series…

1 ÷ 1 = 1
2 ÷ 1 = 2
3 ÷ 2 = 1.5
5 ÷ 3 = 1.6666666667
8 ÷ 5 = 1.6

To extract them, take the reciprocal and then the LCM (lowest common multiple) with 1.

      1∧÷1 2 1.5 1.666666667 1.6
1 1 2 3 5

What happens if you want to accurately compute larger Fibonacci numbers? There’s a limit to the precision based on the internal number representations. Because fibRoger and fibBrian delve into floating point representation, they’re the most limited. (Roger points out that the J language has native support for extended precision and does not suffer from this limitation.)

Dyalog has a system variable, ⎕FR, that sets the how floating-point operations are performed using either IEEE 754 64-bit floating-point operations or IEEE 754-2008 128-bit decimal floating-point operation. For most applications, 64-bit operations are perfectly acceptable, but in applications that demand a very high degree of precision, 128-bit operations can be used, albeit at some performance degradation. Using 64-bit operations, fibBrian loses accuracy after the 35th term while using 128-bits extends the accuracy to 68 terms.

Even setting ⎕FR to use 128-bit operations and the print precision, ⎕PP, to its maximum, we can only compute up to the 164th element 34-digit 8404037832974134882743767626780173 before being forced into scientific notation.

Can we go further easily? Yes we can! Using Dyalog for Microsoft Windows’ .NET integration, we can seamlessly make use of the BigInteger class in the System.Numerics .NET namespace.
First we need to specify the namespace and where to look for it…

      ⎕USING←'System.Numerics,system.numerics.dll'

Then we make a trivial change to our code to use BigInteger instead of native data types…

      fibRecursive←{⍵<3:⍵⍴⎕NEW BigInteger 1 ⋄ {⍵,+/¯2↑⍵}∇ ⍵-1}

And we’re done!

So the next time someone walks up to you and asks “What can you tell me about the 1,000th element of the Fibonacci series?”, you can confidently reply that it has 209 digits and a value of
43466557686937456435688527675040625802564660517371780402481729
08953655541794905189040387984007925516929592259308032263477520
96896232398733224711616429964409065331879382989696499285160037
04476137795166849228875.

For other APL Fibonacci implementations, check out this page.

A Speed-Up Story

The first e-mail of the work week came from Nicolas Delcros. He wondered whether anything clever can be done with ∘.≡ on enclosed character strings. I immediately thought of using “magic functions“, an implementation technique whereby interpreter facilities are coded as dfns. I thought of magic functions because the APL expressions involved in this case are particularly terse:

   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'
   b←⌽a←1↓¨(' '=t)⊂t

   cmpx 'a∘.≡a' '∘.=⍨⍳⍨a'
  a∘.≡a   → 9.69E¯5 |   0% ⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕
  ∘.=⍨⍳⍨a → 8.21E¯6 | -92% ⎕⎕
   cmpx 'a∘.≡b' '(a⍳a)∘.=(a⍳b)'
  a∘.≡b         → 9.70E¯5 |   0% ⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕
  (a⍳a)∘.=(a⍳b) → 1.43E¯5 | -86% ⎕⎕⎕⎕

   y←⌽x←300⍴a

   cmpx 'x∘.≡x' '∘.=⍨⍳⍨x'
  x∘.≡x   → 9.53E¯3 |   0% ⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕
  ∘.=⍨⍳⍨x → 1.52E¯4 | -99% ⎕
   cmpx 'x∘.≡y' '(x⍳x)∘.=(x⍳y)'
  x∘.≡y         → 9.55E¯3 |   0% ⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕
  (x⍳x)∘.=(x⍳y) → 1.95E¯4 | -98% ⎕

The advantage will be even greater if/when we speed up . (And, obviously, the idea is also applicable to ∘.≢ : just replace with and = with .)

Jay Foad objected that the comparisons above aren’t quite fair as the more verbose expressions should check that either ⎕ct←0 or that the arguments do not contain any elements subject to tolerant comparison. I countered that the checking in C would not greatly affect the benchmarks as the time to do the checking is O(m+n) but the benefits are O(m×n) .

Since the factors are so promising and the coding relatively easy, I went ahead and did the work, with the following results:

14.1 14.0 ratio cost of checking
a∘.≡a 9.16e¯6 9.69e¯5 10.58 1.09
a∘.≡b 1.53e¯5 9.70e¯5 6.34 1.08
x∘.≡x 1.48e¯4 9.53e¯3 64.39 1.00
x∘.≡y 1.94e¯4 9.55e¯3 49.23 1.01

The last column (the cost of checking that tolerant comparison is not involved) is under 10% and decreases as the argument sizes increase.

This work is another illustration of the ubiquity and practical usefulness of the “selfie” concept — in the new (14.1) implementation, x∘.≡y is faster when the left and right arguments are the same than when they are not. In particular, the selfie x⍳x or ⍳⍨x occurs twice, and bolsters the point made in item 3 of Sixteen APL Amuse-Bouches:

x⍳x are like ID numbers; questions of identity on x can often be answered more efficiently on x⍳x than on x itself.

Finally, after all that, I recommend that one should consider using x⍳y before using x∘.≡y . The latter requires O(m×n) space and O(m×n) time, and is inherently inefficient.

Selective Expression Tracing

Operator trace from dfns.dws is a simple and low-tech tool for displaying intermediate results of functions, during the evaluation of an expression. For example, what’s going on here?

      (+⌿ ÷ ≢) 1 2 3 4
2.5

We can guess that the above fork is computing the average of the right argument. To see how this works:

      )copy dfns trace
...

⍝ Rename to reduce clutter

      t ← trace

⍝ Bind each "tine" of the fork with t to see what's happening

      (+⌿t ÷t ≢t) 1 2 3 4
≢  1 2 3 4  =>  4
+⌿  1 2 3 4  =>  10
10  ÷  4  =>  2.5
2.5

⍝ But how is that +⌿ reduction evaluated?

      (+t⌿ ÷ ≢) 1 2 3 4
3  +  4  =>  7
2  +  7  =>  9
1  +  9  =>  10
2.5

As a second example, what does the atop (~⍷) do in the following function to remove superfluous '.' characters from its argument?

      {('··'(~⍷)⍵)/⍵} 'squeeze·····me'
squeeze·me

      {('··'(~⍷)t⍵)/⍵} 'squeeze·····me'
··  ~⍷  squeeze·····me  =>  1 1 1 1 1 1 1 0 0 0 0 1 1 1
squeeze·me

…and so forth. Notice how t can be injected to the right of any of the functions in an expression.

For more examples, see the notes for trace.

Data-driven Conditionals

ACKNOWLEDGEMENT: My thanks to Andy Shiers for simplifying the examples

Visitors to APL from other languages are sometimes a bit sniffy about our use of numbers (1 0) for Boolean values True and False. They’re further bemused by our fondness for moving conditional processing out of code and into data. For example, instead of:

   0=2|⍵: ⍵÷2 ⋄ ⍵       ⍝ even number: half

we might say:

   ⍵ ÷ 1+0=2|⍵          ⍝ half if even

At first blush, this seems a bit loony; such transformations appear to produce code which is both slower and more obscure! Here’s another example:

   2|⍵: 1+3×⍵ ⋄ ⍵       ⍝ odd: one more than its triple

vs:

   p ← 2|⍵              ⍝ parity
   p + ⍵ × 1+2×p        ⍝ one more than triple, if odd

But think arrays! Data-driven processing of conditionals is inherently parallel and so performs significantly better for large arrays.

The Collatz conjecture asserts that the sequence “half if even, else one plus triple if odd” always converges to 1. For example, for the number 7 this sequence is:

   7 22 11 34 17 52 26 13 40 20 10 5 16 8 4 2 1

To test this up to a million we can borrow time, a coarse-grained timing operator, from dfns.dws:

      scalar←{                                  
          {                 ⍝ scalar function
              2|⍵:1+3×⍵     ⍝ odd              
              ⍵÷2           ⍝ even             
          }⍣{⍺=1}¨⍵         ⍝ applied under each
      }


      vector←{         
          {                 ⍝ vector function                                              
              p←2|⍵         ⍝ parity                   
              u←⍵÷2-p       ⍝ half of evens            
              v←p+u×1+2×p   ⍝ 1+3× odds                
              ∪v~1          ⍝ without 1s and duplicates
          }⍣{⍺≡⍬}⍵          ⍝ applied to whole vector                           
      }  


      scalar time ⍳1e6      ⍝ scalar processing 8¾ minutes
08:46.03
 
      vector time ⍳1e6      ⍝ array processing 8¾ seconds
08.74

In the above, time and are high order “operators”.

Defined operator time takes a function left operand, which it applies to its argument, returning the execution time in minutes and seconds.

Primitive operator (power or fixpoint) applies its left operand function repeatedly until its right operand function, which is supplied with the result () and argument of each application, returns True (1).

For further discussion of the Collatz function, see http://dfns.dyalog.com/n_osc.htm.

Why Dyalog ’14 was a Shear Delight

rot55

Recent versions of Microsoft Windows support touch screens, which of course means that applications can respond to events originating from touches. Microsoft calls these events “gestures”.

Dyalog decided to add support for gestures to version 14.1 and so projects were planned, designs were designed, code was coded and, at Dyalog ’14 (#Dyalog14), a demonstration was demonstrated. The video of that demonstration is here

The three most interesting gestures are “pan”, “zoom” and “rotate”; I needed a good demo for the rotate gesture. I had an idea but was unsure how to implement it, so I decided to ask the team. Here’s the email I sent (any typos were included in the original email!):

“Does any have code (or an algorithm, but code ‘d be better) to rotate a matrix by an arbitrary angle. The resulting matrix may then be larger than the original but will be filled with an appropriate value.
Does the question makes sense? For the conference I want to demonstrate 14.1 support for Windows touch gestures, one of which is a “rotate” gesture, and I thought that a good demo would be to rotate a CD cover image when the gesture is detected. So I need to be able to rotate the RGB values that make up the image.”

There were a number of helpful responses, mainly about using rotation matrices, but what intrigued me was Jay Foad’s almost throw away comment: “It’s also possible to do a rotation with repeated shears”

That looked interesting, so I had a bit of a Google and found an article on rotation by shearing written by Tobin Fricke. The gist of this article is that shearing a matrix three times (first horizontally, then vertically and then again horizontally) effects a rotation. An interesting aspect of this (as compared to the rotation method) is that no information is lost. Every single pixel is moved and there is no loss of quality of the image.

I’m no mathematician, but I could read my way though enough of that to realise that it would be trivial in APL. The only “real” code needed is to figure how much to shear each line in the matrix.

Converting Fricke’s α=-tan(Θ/2) and ϐ=sin(Θ), we get (using for Θ):

a←-3○⍺÷2
b←1○⍺

and so, for a given matrix, the amounts to shift are given by:

⌊⍺×(⍳≢⍵)-0.5×≢⍵

where is either a or b, depending on which shift we are doing.

Wrapping this in an operator (where ⍺⍺ is either or ):

shear←{
(⌊⍺×(⍳≢⍵)-0.5×≢⍵)⍺⍺ ⍵
}

Our three shears of matrix , of size s, can therefore be written as:

a⌽shear b⊖shear a⌽shear(3×s)↑(-2×s)↑⍵

Note the resize and pad of the matrix before we start so that we don’t replicate the edges of the bitmap. The padding elements are zeros and display as black pixels. They have been removed from the images below for clarity.

The final rotate function is:

rotate←{
shear←{(⌊⍺×(⍳≢⍵)-0.5×≢⍵)⍺⍺ ⍵}
a←-3○⍺÷2 ⋄ b←1○⍺ ⋄ s←⍴⍵
a⌽shear b⊖shear a⌽shear(3×s)↑(-2×s)↑⍵
}

Let’s look at a 45° rotate of a bitmap called “jd”.

jd.CBits←(○.25)rotate jd.CBits

Prior to performing any shears our image is:

rot0

After the first shear it becomes:

rot1

After the second shear it is:

rot2

And after the third shear it is:

rot3

All nicely rotated.

Thursday 25 September at Dyalog ’14

This was originally posted to Catalyst PR’s Facebook page and is reproduced here to make it accessible to people who don’t have Facebook accounts.

For information on the presentations at Dyalog ’14, see http://www.dyalog.com/user-meetings/dyalog14.htm.

On Thursday 25 September we were in for an impressive presentation called School Laser System with NFC Registration and Alerts by Chris “Ziggi” Paul, Laser Learning Ltd (U.K.).

Ziggi has previously developed other applications in Dyalog APL – most notably an ePortfolio-based nursery training solution for the Childcare company (see Dyalog’s case study at http://www.dyalog.com/case-studies/childcare.htm) which won the Nursery Supplier/Innovator 2008 Award at the Nursery Management Today Awards.

Since then the company has grown substantially and have branched out into Carer’s Careers – a company that is focused on training people who work with caring for others – and now the School Laser System with Near Field Communications (NFC).

The background for the School Laser System with NFC is that the principals are now also working with a new Academy free school. The School has 728 pupils, aged 4-11, and they need to be accounted for at all times during the day. The school is located over two main sites and has an external catering hall and a swimming pool. It also provides before-school and after-school clubs and mini-bus services.

The school “LASER” system is designed as a complete communication and information system between parents and school staff; recently it has been expanded to NFC in order to provide teachers and staff with up-to-date, relevant information. Each pupil in the school is issued with an NFC wristband and tags for their bags and key items of the uniform. Staff use NFC-enabled mobile phones or NFC readers connected to PCs to log a pupil at key points, for example, on and off a bus or at the dining hall. The NFC tags let the lunch staff know what food order the parents have requested each day and can give a full itinerary of clubs and school activities. Alerts of allergies or other critical medical information is displayed when a tag is scanned.

Ziggi showed us the parents’ interface to the system, as well as that used by the staff. The security around the system is obviously very high, using encryption and secure log-in. There is no data stored on the wristbands, bag tags or uniform items. If a child has forgotten to bring their wristband, a substitute – which only works for that day – is issued.

I believe that the system has great potential. There is a market not only for other UK infant and junior schools to implement this applications, but catering firms supplying school dinners could also benefit greatly from using this application.

NOTE for readers not familiar with the requirements for UK school staff concerning reporting, tracking, tracing, accounting for, etc. all children at all times – and who may have concerns regarding “big brother” tracking of children’s movements during the time they’re at school or attending after school clubs, etc.: This system is developed for tracking children aged 4-11. The application saves the school staff an enormous amount of admin. Parents love it, and as for issues around data protection, for example, medical information, no data is stored on the items the children carry around (wristband, bag tag or uniform items). It is already common practice that children with allergies and other medical conditions have their own little poster with their picture, age, condition and actions required in case of an incident posted in the classroom and the teacher room.

This is my final post from Dyalog ’14 – see you at Dyalog ’15, which seems to be going in the direction of Italy (TBC).