Power in Cornwall

st-ives-imageAn old folk riddle:

As I was going to St Ives
I met a man with seven wives
Every wife had seven sacks
Every sack had seven cats
Every cat had seven kits
Kits, cats, sacks, wives –
How many were going to St Ives?

One – no calculation required – everyone else was coming from St Ives! From which we see the riddle is older than the railways.

So, what if the narrator met the polygamist on a train? Let’s suppose the polygamist is to be counted too. We start with him. And his wives.

In a scalar language we can loop 5 times, multiplying by 7 each time an initial seed of 1, and adding it to an accumulator.

count = 1; // me
for(i=0;i<5;i++) {
    count+= 7**i;
}
printf(count);

An array solution starts by listing the travellers.

      1 7            ⍝ a man with 7 wives
1 7
      +/1 7          ⍝ count them
8
      +/1 7 49       ⍝ 7 wives had 7 sacks
57

Hang on. There’s a pattern here.

      7*0 1 2        ⍝ 7 raised to consecutive powers
1 7 49
      +/7*0 1 2
57

The APL * means raise to a power. (For multiplication, we use the multiply sign × – how crazy is that?)

      +/7*0 1 2 3 4  ⍝ wives, sacks, cats, kits
2801
      +/7*⍳5         ⍝ ie first 5 integers
2801
      1++/7*⍳5       ⍝ and me
2802

But for vectors u and v and functions f and g, f/u g v is a degenerate case of u f.g v, which leads the way to expressions that work for arrays with more dimensions, so we could equally write:

      1+7+.*⍳5
2802

(Crowded train, I’d say.)

Works in q as well, where the syntax requires parens around the reduction:

q)1+(+/)7 xexp til 5
2802f
q)

You can experiment with the APL expressions in this post at TryAPL.org.

St Ives, by the way, is a charming seaside town in Cornwall. More power to it.

Leaping not looping

Leaping not looping

Today two things. We’ll solve a simple problem in all-at-once array style. And we’ll use an operator to do it. The problem is determining whether a year is a leap year. Simples. But first the operator.

Operator has a special meaning in APL. An operator modifies one function to produce another. (APL here follows the usage of Heaviside, an early user of vectors in mathematics.) For example, the reduction operator / modifies the addition function + so that +/ is the sum function:

      +/1 2 3 4
10

(You can experiment with the expressions in this post at TryAPL.org)

Similarly, ×/ is product and ⌊/ and ⌈/ are respectively minimum and maximum.

    ×/1 2 3 4
24
    ⌊/1 2 3 4
1
    ⌈/1 2 3 4
4

The J language neatly catches the relationship between functions and operators: it calls them verbs and adverbs.

Here is the outer product operator .∘ at work, modifying + and × into functions that return addition and multiplication tables:

      10 20 ∘.+ 1 2 3 4
11 12 13 14
21 22 23 24
      10 20 ∘.× 1 2 3 4
10 20 30 40
20 40 60 80

Our simple problem is to flag which in a list of years are leap years. You probably know the rule: a year is a leap year if it is divisible by 4 — unless it is also divisible by 100, in which case it is not — unless it is also divisible by 400, in which case it is.

Neither 1901 nor 2001 are leap years, because they are not divisible by 4. 1904 and 2004 are both leap years, because they are divisible by 4 and not by 100. 1900 was not a leap year, because it was divisible by 100 but not 400. And 2000 was a leap year because it was divisible by 400.

We can express this quite conventionally in a control structure:

     ∇ Z←isLeap1 years;i
[1]    Z←(≢years)⍴0
[2]    :For i :In ⍳≢years
[3]        :If Z[i]←0=4|years[i]
[4]        :AndIf 0=100|years[i]
[5]            Z[i]←0=400|years[i]
[6]        :EndIf
[7]    :EndFor
     ∇

We’re using the modulus function | to test the remainder after dividing the year by, successively, 400, 100 and 4.

      isLeap1 1900 1901 1904 2000 2001 2004
0 0 1 1 0 1

That succeeds, but uses no array thinking. Array thinking goes more like this. Divide all the years by all the divisors to see the remainders:

      400 100 4∘.| 1900 1901 1904 2000 2001 2004
300 301 304 0 1 4
  0   1   4 0 1 4
  0   1   0 0 1 0

We’re interested only in the sign of the remainders:

      ×400 100 4∘.| 1900 1901 1904 2000 2001 2004
1 1 1 0 1 1
0 1 1 0 1 1
0 1 0 0 1 0

A column for each year, each column three flags. For leap years all the flags must be down (divisible by 400) or only the bottom one (divisible by 4 only). But we can treat each column of flags as a binary number and decode it as a decimal number:

      2⊥×400 100 4∘.| 1900 1901 1904 2000 2001 2004
4 7 6 0 7 6

In this scheme the leap years show up as either 0 or 6. The membership function will flag them:

      (2⊥×400 100 4∘.| 1900 1901 1904 2000 2001 2004) ∊ 0 6
0 0 1 1 0 1

To make a function from this expression, embrace it and replace the data with (‘wha’ever’) placeholder:

      {(2⊥×400 100 4∘.| ⍵) ∊ 0 6}

Or, for a finishing touch, we can use the commute operator to switch the arguments of the membership function and so lose the annoying parentheses:


      isLeap2←{0 6 ∊⍨2⊥×400 100 4∘.| ⍵}
      isLeap2 1900 1901 1904 2000 2001 2004
0 0 1 1 0 1

Now here’s an array-ish thing. You probably spotted the simple pattern in the list of years: the two centuries and the first and fourth years after them. More neatly expressed as an outer product of +:

      1900 2000 ∘.+ 0 1 4
1900 1901 1904
2000 2001 2004

Pleasingly, isLeap2 returns its result in the same shape as its argument:

      isLeap2 1900 2000 ∘.+ 0 1 4
0 0 1
1 0 1

We’ve touched here on three operators: reduce, commute, and outer product. Operators give us great expressive power in describing ways to apply functions across or through arrays. Moreover, you can write your own. For a more elaborate example, see Industrial FP: A case study.

The array thinking characteristically gave us some over-computing: isLeap2 divides every year in the array by three numbers. isLeap1 makes fewer comparisons. On the other hand, its greater code volume has to be parsed by the interpreter, and under the outer product operator, the interpreter has more scope for optimisation or even parallelisation.

Thanks to Tommy Luff for help testing the sample problem.

Let it snow

Let it snow

You might have seen John Scholes’ quietly famous implementation (in Dyalog APL) of Conway’s Game of Life. We’re going to explore just one of the array-language features in it, indexing one array with another.

You can experiment with the expressions at tryapl.org. To animate them you need a Dyalog interpreter, with its editor and delay function.

The principle seems simple enough.

      'abcdef'[1 3 4 6]
acdf

The hidden surprise is that the result of such an expression has the same shape as the array in the index brackets.

      2 2⍴1 3 4 6
1 3
4 6
      'abcdef'[2 2⍴1 3 4 6]
ac
df

snow00

Snow Crash, the title of Neal Stephenson’s acclaimed debut novel, is a reference to the random visual ‘white noise’ on the screen of a completely crashed computer. Let’s make some snow.

      '⎕⌹'[?2 2⍴2]
⌹⌹
⌹⎕
      '⎕⌹'[?2 2⍴2]
⎕⌹
⌹⌹

Now let’s follow John’s example and animate it.

      pic←''
      ⎕ed 'pic'
      {} {pic∘←'⎕⌹'[?30 100⍴2]⊣⎕DL÷⍵ ⋄ ⍵} ⍣40 ⊢8

snow01

Five seconds of ‘snow crash’. Perhaps it would look more snowy with snowflakes – for which we’ll use the five-pointed APL star.

      {} {pic∘←' *'[?30 100⍴2]⊣⎕DL÷⍵ ⋄ ⍵} ⍣40 ⊢8

snow02

Nice. Except real snow doesn’t flicker like that. It falls. It drifts.

Park that thought. Raise the temperature.

Oh, westron wind, when wilt thou blow
That the small rain down may rain?
Christ, that my love were in my arms
And I in my bed again!

wrote Henry VIII. Here we go.

      {} {pic∘←' |'[?30 100⍴2]⊣⎕DL÷⍵ ⋄ ⍵} ⍣40 ⊢8

snow03

Uh. Wasn’t that wind supposed to be blowing from the west?

{} {pic∘←' \'[?30 100⍴2]⊣⎕DL÷⍵ ⋄ ⍵} ⍣40 ⊢8

snow04

Enough rain. Back to snowflakes. And fewer of them.

      pic←' *'[1+1=?30 100⍴20]

snow05

Flakes must fall.We need a skyfall algorithm. Drop the bottom row, prefix new flakes at the top:

      pic←{' *'[1+1=?100/20]⍪¯1 0↓⍵} pic

snow06

And some animation.

      next←{' *'[1+1=?100/20]⍪¯1 0↓⍵} 
      {} {pic∘←next ⍵⊣⎕DL÷8} ⍣40 ⊢pic

This snow falls straight down, which isn’t quite right. We expect a little wind. From the west, perhaps.

next←{
   (L T)←{' *'[1+⍵]}¨1=?((⍴⍵)-0 1)⍴¨20
   L,T⍪¯1 ¯1↓⍵
} 
{} {pic∘←next ⍵⊣⎕DL÷8} ⍣200 ⊢pic

The snowfield is 3-D. Nearer flakes should be larger – and fewer – than distant ones.

next←{
   (L T)←{'⍟**∘∘∘...... '[⍵]}¨13⌊?(0 ¯1+⍴⍵)⍴¨100
   L,T⍪¯1 ¯1↓⍵
}
{} {pic∘←next ⍵⊣⎕DL÷8} ⍣200 ⊢pic

snow08

There we go. 25 seconds of simulated snowfall.

We’re missing some features. Turbulence produces a certain amount of random jiggling for the flakes. Distant flakes should fall more slowly than near flakes. We’ve reached the limit of our simple representation of a snowfield as 2-D character array. We might (in a later post) represent the flakes differently, and project their representations onto a visual plane.

Watch the whole thing on YouTube!

For now, a Merry Christmas and – let it snow!

Arrays, actually

Arrays, actually

What is it about the array languages? What does it mean to ‘think in arrays’? This blog aims to exhibit what is distinctive about the array languages descended from Iverson Notation, originally developed at Harvard in the 1950s.

We’re not trying here to teach you these languages. We suppose you already know a few languages and are curious about these. We aim to satisfy your curiosity and help you decide whether you want to learn and use array languages. We will explain what we think are the key concepts in the code we show.

The languages shown here are Dyalog APL, J and kdb. You can experiment with APL online at TryAPL.org.

Dan Baronet & Stephen Taylor