leah blogs: December 2021

Yesterday, Advent of Code had an interesting problem: given the heightmap of a lava cave, compute the lowest points and the size of their basins (connected regions).

Let’s do this in BQN again, as this problem teaches some good ways to think in array languages.

First, let’s load the input data into a matrix:

   d ← > '0' -˜ •FLines"day09"

We subtract the ASCII lines from the character 0 to get numerical rows. The merge function (>) then converts this list-of-lists into a 2-dimensional array. For the sample data, we get:

┌─
╵ 2 1 9 9 9 4 3 2 1 0
  3 9 8 7 8 9 4 9 2 1
  9 8 5 6 7 8 9 8 9 2
  8 7 6 7 8 9 6 7 8 9
  9 8 9 9 9 6 5 6 7 8
                      ┘

A low point is a point that is lower than every orthogonally adjacent point. Thanks to array progamming, we can solve this for the whole matrix at once without any loops!

The core idea is to shift the array into each cardinal direction, and then compute the minimum of these arrays. If the original array is smaller then the array of the minimums, it’s a low point.

By default, shifting («, ») in BQN inserts zeroes for numerical arrays. But since we are looking for the minimum, we need to shift in a value that is higher than any. We can simply use ∞.

So, to shift in ∞ from the left, we use:

   ∞»˘d
┌─
╵ ∞ 2 1 9 9 9 4 3 2 1
  ∞ 3 9 8 7 8 9 4 9 2
  ∞ 9 8 5 6 7 8 9 8 9
  ∞ 8 7 6 7 8 9 6 7 8
  ∞ 9 8 9 9 9 6 5 6 7
                      ┘

Unfortunately, shifting from the top is not so easy:

   ∞»d
Error: shift: =𝕨 must be =𝕩 or ¯1+=𝕩 (0≡=𝕨, 2≡=𝕩)
at ∞»d

We need would need to make a list of ∞ long enough to agree with the array width:

   (∞¨⊏d)»d
┌─
╵ ∞ ∞ ∞ ∞ ∞ ∞ ∞ ∞ ∞ ∞
  2 1 9 9 9 4 3 2 1 0
  3 9 8 7 8 9 4 9 2 1
  9 8 5 6 7 8 9 8 9 2
  8 7 6 7 8 9 6 7 8 9
                      ┘

However, since we need to do this on every side, we can also look at the problem differently: we shift in from the left under rotation by 0, 90, 180, 270 degrees.

How do we “rotate” a matrix? We reverse the rows (⌽) and then transpose (⍉) it.

   ⍉⌽d
┌─
╵ 9 8 9 3 2
  8 7 8 9 1
  9 6 5 8 9
  9 7 6 7 9
  9 8 7 8 9
  6 9 8 9 4
  5 6 9 4 3
  6 7 8 9 2
  7 8 9 2 1
  8 9 2 1 0
            ┘

By using the repeat modifier (⍟) we can easily rotate several times.

   (⍉∘⌽⍟(↕4)) d
┌─
· ┌─                      ┌─            ┌─                      ┌─
  ╵ 2 1 9 9 9 4 3 2 1 0   ╵ 9 8 9 3 2   ╵ 8 7 6 5 6 9 9 9 8 9   ╵ 0 1 2 9 8
    3 9 8 7 8 9 4 9 2 1     8 7 8 9 1     9 8 7 6 9 8 7 6 7 8     1 2 9 8 7
    9 8 5 6 7 8 9 8 9 2     9 6 5 8 9     2 9 8 9 8 7 6 5 8 9     2 9 8 7 6
    8 7 6 7 8 9 6 7 8 9     9 7 6 7 9     1 2 9 4 9 8 7 8 9 3     3 4 9 6 5
    9 8 9 9 9 6 5 6 7 8     9 8 7 8 9     0 1 2 3 4 9 9 9 1 2     4 9 8 9 6
                        ┘   6 9 8 9 4                         ┘   9 8 7 8 9
                            5 6 9 4 3                             9 7 6 7 9
                            6 7 8 9 2                             9 8 5 6 9
                            7 8 9 2 1                             1 9 8 7 8
                            8 9 2 1 0                             2 3 9 8 9
                                      ┘                                     ┘
                                                                              ┘

Finally, we perform the shift operation under (⌾) the rotation, that is, BQN rotates the array, does the shift, and knows how to undo the rotation!

┌─
· ┌─                      ┌─                      ┌─                      ┌─
  ╵ ∞ 2 1 9 9 9 4 3 2 1   ╵ 3 9 8 7 8 9 4 9 2 1   ╵ 1 9 9 9 4 3 2 1 0 ∞   ╵ ∞ ∞ ∞ ∞ ∞ ∞ ∞ ∞ ∞ ∞
    ∞ 3 9 8 7 8 9 4 9 2     9 8 5 6 7 8 9 8 9 2     9 8 7 8 9 4 9 2 1 ∞     2 1 9 9 9 4 3 2 1 0
    ∞ 9 8 5 6 7 8 9 8 9     8 7 6 7 8 9 6 7 8 9     8 5 6 7 8 9 8 9 2 ∞     3 9 8 7 8 9 4 9 2 1
    ∞ 8 7 6 7 8 9 6 7 8     9 8 9 9 9 6 5 6 7 8     7 6 7 8 9 6 7 8 9 ∞     9 8 5 6 7 8 9 8 9 2
    ∞ 9 8 9 9 9 6 5 6 7     ∞ ∞ ∞ ∞ ∞ ∞ ∞ ∞ ∞ ∞     8 9 9 9 6 5 6 7 8 ∞     8 7 6 7 8 9 6 7 8 9
                        ┘                       ┘                       ┘                       ┘
                                                                                                  ┘

Now we insert (´) the minimum function (⌊) between these arrays and compute the minimum at each position:

   ⌊´{∞⊸»˘⌾(⍉∘⌽⍟𝕩)d}¨↕4
┌─
╵ 1 2 1 7 4 3 2 1 0 1
  2 1 5 6 7 4 3 2 1 0
  3 5 6 5 6 7 4 7 2 1
  7 6 5 6 7 6 5 6 7 2
  8 7 6 7 6 5 6 5 6 7
                      ┘

The positions where the original array d is still smaller are the low points, and we store them for part 2.

   l ← d < ⌊´{∞⊸»˘⌾(⍉∘⌽⍟𝕩)d}¨↕4
┌─
╵ 0 1 0 0 0 0 0 0 0 1
  0 0 0 0 0 0 0 0 0 0
  0 0 1 0 0 0 0 0 0 0
  0 0 0 0 0 0 0 0 0 0
  0 0 0 0 0 0 1 0 0 0
                      ┘

To finish part 1, we need to compute the risk level for each low point, which is 1 plus the height. So compute that:

   (1+d)×l
┌─
╵ 0 2 0 0 0 0 0 0 0 1
  0 0 0 0 0 0 0 0 0 0
  0 0 6 0 0 0 0 0 0 0
  0 0 0 0 0 0 0 0 0 0
  0 0 0 0 0 0 6 0 0 0
                      ┘

Finally, we compute the end result by deshaping (⥊) the array into a single long list and summing it up:

   +´⥊(1+d)×l
15

This concludes part 1.

Part 2 is less straight forward. We need to compute the basins around every low point, which is the area limited by the points of height 9.

Since we need to compute the size for each basin in the end, we need to know which point belongs to which basin. To get started, we first give each low point a unique number. One way to do this is to assign an index to each position and just use those that are used in the low point array.

In BQN, the shape (≢) of an array is the list of sizes for each axis:

   ≢d
⟨ 5 10 ⟩

We count up to the product of these values and add one (to avoid numbering a basin with 0, which will be used for the basin limits). Then, we multiply this with the array of low points so all ones in it get turned into a unique basin index. We perform the multiplication under deshape, so we keep the shape of the input data:

   s ← (1+↕×´≢d) ×⌾⥊ l
┌─
╵ 0 2  0 0 0 0  0 0 0 10
  0 0  0 0 0 0  0 0 0  0
  0 0 23 0 0 0  0 0 0  0
  0 0  0 0 0 0  0 0 0  0
  0 0  0 0 0 0 47 0 0  0
                         ┘

Here’s the main idea how to solve part 2: we incrementally grow these areas by adding their neighbors until the whole array is filled. For one step, we do this with a function Rise:

   Rise ← (d≠9)⊸×(»⌈«⌈«˘⌈»˘⌈⊣)

Here, shifting in zeroes is good enough, so we can use the monadic shift functions. The train at the end computes the maximum (⌈) of the four shifted versions and the array itself (⊣). We multiply it with a matrix that is 0 where the depth is 9, so the basin limits will be constantly zero.

Let’s run it once and twice to see how it works:

   Rise s
┌─
╵ 2  2  0  0 0  0  0  0 10 10
  0  0 23  0 0  0  0  0  0 10
  0 23 23 23 0  0  0  0  0  0
  0  0 23  0 0  0 47  0  0  0
  0  0  0  0 0 47 47 47  0  0
                              ┘
   Rise⍟2 s
┌─
╵ 2  2  0  0  0  0  0 10 10 10
  2  0 23 23  0  0  0  0 10 10
  0 23 23 23 23  0  0  0  0 10
  0 23 23 23  0  0 47 47  0  0
  0  0  0  0  0 47 47 47 47  0
                               ┘

As you can see, row 1 colum 3 does not get filled by basin #2 since it has height 9.

Now we could just iterate this step a often enough, or actually only until we reach a fixpoint; that is, applying Rise again doesn’t change the value anymore.

A simple way to implement a fixpoint operator is

_fix ← { 𝕩 ≡ 𝔽 𝕩 ? 𝕩 ; 𝕊 𝔽 𝕩 }

Let’s compute the filled map:

   Rise _fix s
┌─
╵  2  2  0  0  0 10 10 10 10 10
   2  0 23 23 23  0 10  0 10 10
   0 23 23 23 23 23  0 47  0 10
  23 23 23 23 23  0 47 47 47  0
   0 23  0  0  0 47 47 47 47 47
                                ┘

Now, we just need to compute the sizes of the basins, which means computing a histogram of the basin numbers. We deshape the array again, and only keep all values bigger than zero:

   m ← ⥊ Rise _fix s
   (m>0)/m
⟨ 2 2 10 10 10 10 10 2 23 23 23 10 10 10 23 23 23 23 23 47 10 23 23 23 23 23 47 47 47 23 47 47 47 47 47 ⟩

With the nice group indices function (⊔), we can count how often each value appears:

   ≠¨⊔(m>0)/m
⟨ 0 0 3 0 0 0 0 0 0 0 9 0 0 0 0 0 0 0 0 0 0 0 0 14 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 9 ⟩

Finally, let’s compute the three largest values by sorting in descending order (∨) and taking the first three entries (↑):

   3↑∨≠¨⊔(m>0)/m
⟨ 14 9 9 ⟩

We then multiply this together and get the result for part 2:

   ×´3↑∨≠¨⊔(m>0)/m
1134

Yesterday, Advent of Code had an interesting problem: modeling the growth rate of lanternfish. Read the link for the full details, but the core algorithm is this:

Each day, a 0 becomes a 6 and adds a new 8 to the end of the list, while each other number decreases by 1 if it was present at the start of the day.

For part 1, you need to compute the total number of lanternfish after 80 days, and for part 2 after 256 days. While the first part is possible with a native list based version (around 300k), the second part yields over 1.6 trillion for my input, so storing every lanternfish on its own is infeasible.

Luckily, with a bit of thinking we quickly realize that the position and order of the lanternfish is irrelevant and we can thus just count how many at each age we have and sum them up at the end for the answer.

I decided to tackle this problem using the array language BQN and will try to explain how it works.

First, let’s convert the sample input data into a histogram: For reasons that will be obvious in a minute, we compare each input element with the list of numbers from 0 to 8 (↕9).

   input ← ⟨3,4,3,1,2⟩
   input =⌜ ↕9
┌─
╵ 0 0 0 1 0 0 0 0 0
  0 0 0 0 1 0 0 0 0
  0 0 0 1 0 0 0 0 0
  0 1 0 0 0 0 0 0 0
  0 0 1 0 0 0 0 0 0
                    ┘

If we sum this table by column now, we have a count which age appears how often (note that ages over 8 never happen):

   d ← +˝ input =⌜ ↕9
⟨ 0 1 1 2 1 0 0 0 0 ⟩

Next, we need to figure out what happens in a step. We can model one part (number decreases by 1, or becomes 8 when it was 0) by rotating (⌽) the array to the left:

   1 ⌽ d
⟨ 1 1 2 1 0 0 0 0 0 ⟩
   2 ⌽ d
⟨ 1 2 1 0 0 0 0 0 1 ⟩

But we also need to add the first element to column 6, e.g. by adding an appropriately scaled vector. First, the vector, which we get by taking the numbers of 0 to 8 and comparing them to 6:

   6=↕9
⟨ 0 0 0 0 0 0 1 0 0 ⟩

Now we can scale this vector by the first element (⊑):

   (6=↕9) × ⊑ ⟨9,8,7,6,5,4,3,2,1⟩
⟨ 0 0 0 0 0 0 9 0 0 ⟩

Finally, we put both things together. We can use a hook if we bind (⊸) the argument of the rotation to the operator:

Step ← 1⊸⌽ + (6=↕9)×⊑

Now, we simply need to repeat (⍟) this step 80 (resp. 256) times and count up the results to solve the puzzle!

   +´ Step⍟80 d
5934
   +´ Step⍟256 d
26984457539

This agrees with the provided example.

Of course, computing 256 steps is very fast, but we can wonder if there is not a more efficient solution. Let’s say we wanted to compute many lanternfish populations quickly.

The key realization is that Step is a linear function, which means:

(n × Step v) ≡ (Step n × v)
(Step v) + (Step w) ≡ (Step v + w)

We can thus compute Step for every basis vector of our 9-dimensional histogram space; let’s take the unit vectors here:

   step256 ← {+´Step⍟256 𝕩}˘ =⌜˜↕9
⟨ 6703087164 6206821033 5617089148 5217223242 4726100874
  4368232009 3989468462 3649885552 3369186778 ⟩

And now we can solve part 2 for any input by using a dot-product (pairwise multiplication followed by addition):

   step256 +´∘× d
26984457539

This yields the same result as above, but only does 9 multiplications and 9 additions.

We can also compute step256 faster by using a bit of matrix theory. Since we know the operation of step on the basis vectors, we can compute the matrix m corresponding to this linear operator:

   m ← Step˘ =⌜˜↕9
┌─
╵ 0 0 0 0 0 0 1 0 1
  1 0 0 0 0 0 0 0 0
  0 1 0 0 0 0 0 0 0
  0 0 1 0 0 0 0 0 0
  0 0 0 1 0 0 0 0 0
  0 0 0 0 1 0 0 0 0
  0 0 0 0 0 1 0 0 0
  0 0 0 0 0 0 1 0 0
  0 0 0 0 0 0 0 1 0
                    ┘

Note how we just applied Step for every column of the identity matrix.

We can solve the problem the same way now by computing the 256th power of the matrix and multiplying it to the histogram vector and summing up again.

First, we need the matrix product:

   MP ← +˝∘×⎉1‿∞

Naively, we can compute the power like this (we already have m to the power of 1!):

   m256 ← m MP⍟255 m

Now multiply d to it, and sum up:

   +´ m256 +˝∘× d
26984457539

Note that step256 is just the sum of rows, so we can precompute that:

   step256 ≡ +˝˘ m256
1

If you paid attention, you may now wonder why we replaced 256 row operations with 255 matrix multiplications and claim it’s faster to compute. We need an additional trick: fast matrix exponentiation. Instead of 255 matrix multiplications we just need 8 matrix squarings with itself to compute the 256th power:

   m256 ≡ MP˜⍟8 m
1

This yields an asymptotically faster algorithm to compute the step vector for higher powers. For non-powers of two you need to implement a square-and-multiply algorithm, or look up the optimal result. This is left as an exercise to the reader.

For day 80 we can use only 7 matrix multiplications:

   (MP˜⍟4 (m MP MP˜⍟2 m)) ≡ (m MP⍟79 m)
1

In theory, you could optimize this even further and compute a closed form from the eigenvectors of m, but you’ll end up with nonic roots and I have not found a way to compute them precise enough to be useful. So let’s leave it at that for now.