A coding impromptu

I call it infamous because it made me feel stupid twice: first when I encountered it in code, and again when I read its docs. To understand its behaviour, you need to be familiar with quite a bit of BQN, especially the functional programming and combinators aspects. As a newbie at the time, I found it quite daunting. It took me about five scattered attempts over several months to get it. Looking back, the difficulty wasn't so much BQN's syntax, but my struggle to express complex recursion, which modifier recursion definitely is.

An unrolling of the first two steps reveals that up to 2⋆n evaluations of 𝔽 can occur at recursion level n. This is derived by noting that, within the BQN combinator, the left function of the rightmost atop dictates that the 𝔽 for the subsequent step is, in accordance to 𝔽_𝕣_𝔾:

_w0_ ← {𝔽⍟𝔾∘𝔽_𝕣_𝔾∘𝔽⍟𝔾𝕩}
_w1_ ← {(𝔽⍟𝔾∘𝔽)⍟𝔾∘(𝔽⍟𝔾∘𝔽)_w0_𝔾∘(𝔽⍟𝔾∘𝔽)⍟𝔾𝕩}
_w2_ ← {((𝔽⍟𝔾∘𝔽)⍟𝔾∘(𝔽⍟𝔾∘𝔽))⍟𝔾∘((𝔽⍟𝔾∘𝔽)⍟𝔾∘(𝔽⍟𝔾∘𝔽))_w0_𝔾∘((𝔽⍟𝔾∘𝔽)⍟𝔾∘(𝔽⍟𝔾∘𝔽))⍟𝔾𝕩}

Another way to clarify the concept is to implement the same logic both as a function and as a 1-modifier, and then compare these implementations with the two 2-modifiers (one exhibiting linear and the other a logarithmic number of stack frames):

Whiles ← {F‿G𝕊𝕩:
  Wfun ← {𝕎⍟G∘𝕎˙⊸𝕊∘𝕎⍟G𝕩}
  _wom ← {𝔽⍟G∘𝔽_𝕣∘𝔽⍟G𝕩}
  _wtmlog_ ← {𝔽⍟𝔾∘𝔽_𝕣_𝔾∘𝔽⍟𝔾𝕩}
  _wtmlin_ ← {𝕊∘𝔽⍟𝔾𝕩}
  ⟨f Wfun 𝕩, f _wom 𝕩, f _wtmlog_ g 𝕩, f _wtmlin_ g ⎊"SO"𝕩⟩
}

Let’s test it with a simple iteration that exceeds CBQN’s recursion limit, triggering a stack overflow:

⟨1⊸+, 5000⊸≥⟩ Whiles 0

⟨ 5001 5001 5001 "SO" ⟩

One thing that has always bothered me in BQN is the absence of a high-performance matrix product implementation. While the current inner product idiom guarantees numerical accuracy, it remains quite slow for larger matrices. Of course, wrapping BLAS directly is always an option:

blasFFI ← (⊣•FFI·(1-˜+`׬)∘=⟜⊏⊸⊔⊢)´ ⟨
  "/lib/libcblas.so"∾˜•BQN 1⊑•Sh "nix-instantiate"‿"--eval-only"‿"--expr"‿"(import <nixpkgs> {}).blas.outPath"
  " & cblas_dgemm u32 u32 u32 i32 i32 i32 f64 *f64 i32 *f64 i32 f64 &f64 i32"
⟩
Dgemm ← {BlasFFI 101‿111‿111‿m‿n‿k‿1‿𝕨‿k‿𝕩‿n‿0‿(m‿n⥊0)∾⊢´m‿k‿·‿n←𝕨∾○≢𝕩}

This approach has roughly the same overhead as NumPy's dot, so anyone applauding NumPy's performance is essentially praising BLAS. However, it is possible to implement a simple blocked matrix product directly in BQN, using a basic back-of-the-envelope calculation for block sizes (sorry, no fancy optimized kernels here). For matrix sizes exceeding my machine's cache, speed-ups of about six times are achievable:

mat‿mbt ← ⟨⋈˜2⥊500, ⋈˜5⥊600⟩ /¨⊸⊔¨ ma‿mb ← •rand.Range⟜0¨1e3×⟨1‿1, 3‿3⟩
>⟨ma‿ma‿mat, mb‿mb‿mbt⟩ {𝕎˜•_timed𝕩}¨¨˜ <⟨Dgemm, +˝∘×⎉1‿∞, ∾(+˝+˝∘×⎉1‿∞¨)⎉1‿∞⟩

┌─                                                            
╵         0.008988871        0.646108393 0.37081367400000004  
  0.16528436400000002 45.110128999000004   7.460860705000001  
                                                             ┘

Moreover, there is only a modest 10-character leap from +˝∘×⎉1‿∞ to ∾(+˝+˝∘×⎉1‿∞¨)⎉1‿∞. As a bonus, here's a convenient function to compute powers of a square matrix 𝕩 (particularly useful in graph theory) using blocks of size 𝕨, assuming it evenly divides the matrix dimensions:

MPow ← {∾(+˝+˝∘×⎉1‿∞¨)⎉1‿∞˜𝕩⊔˜/¨⥊⟜𝕨¨⌈𝕨÷˜≢𝕩}

Here I could have used a fancier but slower under 𝔽˜⌾((/¨⥊⟜𝕨¨⌈𝕨÷˜≢𝕩)⊸⊔). A naïve search for the optimal block size yields:

50‿100‿300‿600‿1000‿1500 {𝕨⊸MPow•_timed𝕩}¨ <3e3‿3e3 •rand.Range 0

⟨ 15.883247281000001 12.451545102 8.224263037 7.7123262310000005 15.417112082000001 34.447925623 ⟩

For more fine-grained control, here is the code for arbitrary block sizes (not necessarily divisors of the length), which pads the matrices with zeros. Let's take a look at the local minimum:

MPow ← {𝕩≢⊸↑∾(+˝+˝∘×⎉1‿∞¨)⎉1‿∞˜𝕩(⥊⟜𝕨¨∘⊢/¨⊸⊔𝕨⊸×↑⊣)⌈𝕨÷˜≢𝕩}
(300+50×↕8) {𝕨⊸MPow•_timed𝕩}¨ <3e3‿3e3 •rand.Range 0

⟨ 8.30279774 10.112563361000001 9.781014477000001 9.670085717000001 7.556631647000001 10.970897867000001 7.570657628 10.231164773000001 ⟩

The blocked algorithm is only about 2 seconds slower than the analogous inner product in Dyalog APL 19 +.×⌿?2 3000 3000⍴1e10, whereas the standard one is 10 times slower.

This is a very efficient procedure that finds prefix strings in linear time. The imperative implementation reads:

ZI ← {𝕊s:
  l‿r‿z ← 0⚇0 0‿0‿s
  z ⊣ {
    v ← r(⊢1⊸+•_while_{(𝕩+𝕨)<≠s ? =´⟨𝕩,𝕩+𝕨⟩⊑¨<s ; 0}<◶({z⊑˜𝕩-l}⌊-+1)‿0)𝕩
    r <◶@‿{𝕊: l↩𝕩-v+1 ⋄ r↩𝕩} 𝕩+v-1
    z v⌾(𝕩⊸⊑)↩
  }¨ ↕≠s
}
ZI "abacabadabacaba"

⟨ 15 0 1 0 3 0 1 0 7 0 1 0 3 0 1 ⟩

Two algorithmic improvements can be made here, namely only iterate over indices where the character found is equal to the first character, and only search to extend the count if it goes up to the end of r:

ZFun ← {𝕊s:
  CountEq ← { 1⊸+•_while_((≠𝕨)⊸≤◶⟨⊑⟜𝕨≡⊑⟜𝕩,0⟩) 0 }
  l←r←0 ⋄ Ulr ← {(r⌈↩𝕨+𝕩)>r ? l↩𝕨 ⋄ 𝕩; 𝕩}
  SearchEq ← ⊣ Ulr ⊢ + + CountEq○(↓⟜s) ⊢
  Set ← {i𝕊𝕩: ((r-i) (i SearchEq 0⌈⊣)⍟≤ (i-l)⊑𝕩)⌾(i⊸⊑) 𝕩 }
  (⌽1↓/⊑⊸=s) Set´˜ ↑˜≠s
}

I came up with two array versions, with quadratic and cubic time complexities respectively:

ZAQ ← ¯1↓↓(+´·∧`⊣=≠⊸↑)¨<
ZAC ← (+´∧`)¨<=↕∘≠{«⍟𝕨𝕩}⌜<
(ZAQ≡ZAC)◶@‿ZAC "abacabadabacaba"

⟨ 15 0 1 0 3 0 1 0 7 0 1 0 3 0 1 ⟩

With further refinements, the earlier solutions can be transformed into:

ZAQ‿ZAC ← {(+´∧`)¨𝕏}¨ ⟨≠↑↓=⌽∘↑, <=«⍟(↕∘≠)⟩

This problem can be solved in \(O(n\log n)\) using dynamic programming. Here is an imperative implementation which is quadratic, but can be optimized:

LISI ← {
  k‿dp ← ¯1‿(∞¨𝕩)
  {i ← ∧´◶(⊑⊐⟜0)‿{𝕊:k+↩1} dp<𝕩 ⋄ dp 𝕩⌾(i⊸⊑)↩}¨ 𝕩
  +´∞>dp
}
LISI¨ ⟨0‿1‿0‿3‿2‿3, 10‿9‿2‿5‿3‿7‿101‿18, 7‿7‿7‿7‿7⟩

⟨ 4 4 1 ⟩

A more elegant array solution is:

LISA ← +´∞≠∞¨{𝕨⌾((⊑𝕩⍋𝕨-1)⊸⊑)𝕩}´⌽
LISA¨ ⟨0‿1‿0‿3‿2‿3, 10‿9‿2‿5‿3‿7‿101‿18, 7‿7‿7‿7‿7⟩

⟨ 4 4 1 ⟩

Let's )explain this optimized version, so we can truly appreciate its beauty:

 +´∞≠∞¨{𝕨⌾((⊑𝕩⍋𝕨-1)⊸⊑)𝕩}´⌽ 
 │ │ │ ││    │ │ │  │ │  │ 
 │ │ │ {┼────┼─┼─┼──┼─┼─´│ 
 │ │ ∞¨ │    │ │ │  │ │ ││ 
 │ │  └─┼────┼─┼─┼──┼─┼─┼⌽ 
 │ ∞≠───┼────┼─┼─┼──┼─┼─┘  
 +´ │   │    │ │ │  │ │    
  └─┘   │    │ │ │  │ │    
        │    │ 𝕨-1  │ │    
        │    𝕩⍋─┘   │ │    
        │   ⊑─┘     │ │    
        │   └──────⊸⊑ │    
        𝕨⌾─────────┘  │    
         ├────────────𝕩    
╶────────┘

The full expression is structured as a two-train: we sum all finite entries from the result of the rightmost three-train. The three-train is a right fold over the reversed input, with an initial array of ∞ and the same length as the input. In each step of the fold, we modify the right argument using under: we perform a binary search with strict comparison to find where the next element should go. The element is either placed at the end of the unfilled region, or it replaces the first element that is greater than 𝕨. Since BQN uses a based array model, we pick the enclosed atom from this operation's result. So it goes³.

This problem is the archetypal example of backtracking. Initially, I tried to solve it using a function to place the queens in the full board, hoping that it would lead to a more array oriented solution:

8 {((∨⌜´0⊸=)∨(0=-⌜´)∨0=+⌜´) 𝕩-¨<↕𝕨} 2‿3

┌─                 
╵ 0 1 0 1 0 1 0 0  
  0 0 1 1 1 0 0 0  
  1 1 1 1 1 1 1 1  
  0 0 1 1 1 0 0 0  
  0 1 0 1 0 1 0 0  
  1 0 0 1 0 0 1 0  
  0 0 0 1 0 0 0 1  
  0 0 0 1 0 0 0 0  
                  ┘

This resulted in a more complicated algorithm, so I decided to go for the classical Wirth implementation:

NQ ← {𝕊n:
  V‿P ← {⊣𝕏(⊢∾-⋈+)´∘⊢}¨ ⟨∨´⊑¨˜, {1⌾(𝕩⊸⊑)𝕨}¨⟩
  {n≠𝕩 ? +´(𝕨V⊢)◶⟨(𝕩+1)𝕊˜𝕨P⊢,0⟩∘(𝕩⋈⊢)¨ ↕n ; 1
  }˜´ (0⋈0×·↕¨⊢∾·⋈˜+˜)n 
}

Which nicely compares with the OEIS sequence:

a000170 ← 1‿0‿0‿2‿10‿4‿40‿92
a000170 ≡ NQ¨ 1+↕8

1

And of course, in the implementation above I could have used a single array instead of three, but I find the resulting validation and position functions very aesthetic the way they are.

The Boyer–Moore algorithm allows for finding the majority element (element that appears more than ⌊𝕩÷2 times in the array) in linear time. If such element exists, then it is equal to the mode of the data, and for this task we have a nice array solution. The original implementation could be expressed as:

BM ← {v←0 ⋄ I←⊢⊣=◶{𝕊:v+↩1}‿{𝕊:v-↩1} ⋄ 0{𝕊:v=0}◶⟨I,I˜⊣⟩´𝕩}
BM 6‿1‿3‿1‿3‿3‿4‿3‿3‿5

3

The previous fold tracks the majority element as state, a more elegant approach maintains the number of votes:

BM ← {e←@ ⋄ 0{𝕩=0 ? e↩𝕨⋄1 ; 𝕩+¯1⋆e≢𝕨}´𝕩 ⋄ e}
BM 6‿1‿3‿1‿3‿3‿4‿3‿3‿5

3

Some time ago, while working on performance optimization of linear algebra operations with Boolean arrays, I encountered an interesting summation property for an array \(a\) of length \(n\):

It turns out that the RHS can be elegantly transformed into a scan, giving rise to a beautiful identity that applies to all natural numbers, not just Booleans as I initially thought:

(+`≡·+´/≤⟜<⊒˜) •rand.Range˜ 1e3

1

This identity holds because ⊒˜ represents the indices i of the list, and since +´(/𝕩)=i ←→ i⊑𝕩, the fold sums all the elements in 𝕩 up to i, for i in the range of the length of the list. Ergo, a scan.

Studying tree algorithms in APL, I learned about the depth vector representation. If the nested object in consideration is a string, the best approach is using boolean masks. However, when dealing with a BQN list, recursion becomes necessary to determine the depth of nested elements. Here’s how it can be implemented:

{=◶⟨⋈0, 1+·∾𝕊¨⟩𝕩} ⟨1, ⟨2, ⟨3⟩, ⟨4, ⟨5, ⟨6, 7⟩⟩⟩⟩, 1⟩

⟨ 1 2 3 3 4 5 5 1 ⟩

This metric is one of the reasons for the deplorable state of modern academia, and the headaches for outsiders trying to get in. Consider that Peter Higgs has an estimated h-index of only 12. By contrast, a random professor nowadays boasts an h-index ten times as high, and exponentially less impact. Enough of ranting, let's concentrate on finding an elegant way to implement this useless thing:

HL ← (+´∘«⊒˜≤+`⌾⌽)·/⁼≠⊸⌊
HS ← +´∨≥1+⊒˜
(HL≡HS)◶@‿HL 14‿14‿11‿9‿5‿5‿1‿1‿1‿1‿0

5

If someone ever published that much, sorting will eventually be slower:

HL‿HS {𝕎•_timed𝕩}¨< 1e8 •rand.Range 1e3

⟨ 0.083824959 0.21801262700000001 ⟩

A testament to the idea that the simplest solution in BQN is often the most efficient: I initially clip my citations array with {≠¨⊔≠∘𝕩¨⌾(≥⟜≠∘𝕩⊸/)𝕩}, which is just /⁼≠⊸⌊.

This is a classical interview problem that can be solved in linear time. Interestingly, it admits a very elegant array solution:

(+´⊢-˜⌈`⌾⌽⌊⌈`) [0,1,0,2,1,0,1,3,2,1,2,1]

6

That is, we take the minimum of max-scans from the left and from the right, and subtract the corresponding height. Reducing the resulting array gives the amount of trapped water.

A closely related problem is container with most water, which unfortunately is not so easy to solve in linear time using an array approach (one can easily implement the imperative two pointers solution in BQN, but it will probably be slow). Here are two solutions, one \(O(n^2)\) and the other \(O(n\log n)\), both tacit:

⟨⌈´∘⥊⌊⌜˜×·-⌜˜⊒˜, ⌈´∨×(⌈`⊸-⌈⊢-⌊`)∘⍒⟩ {10 𝕎•_timed𝕩}¨< •rand.Range˜1e4

⟨ 0.080050875 4.14558e¯5 ⟩

The Levenshtein (or edit) distance is a measure of the similarity between two strings. It is defined by the following recurrence, which is the basis of dynamic programming algorithms like Wagner-Fisher:

There is an elegant implementation of a variation of the Wagner–Fischer algorithm in the BQNcrate. It has been particularly challenging for me to understand it—not due to the clarity of the primitives, but rather because of the clever transformation employed. I believe that this variant can be derived by shifting the distance matrix. Given two strings \(s\) and \(t\) of lengths \(n\) and \(m\), respectively, we define a new distance matrix as follows:

Under this transformation, the recurrence relation becomes:

The above recurrence can be easily identified in the 3-train's middle function, which is folded over the table of the costs (table comparing the characters). Note that we compare insertions and substitutions, and then we can do a min scan over the result to get the deletions, which gives a vectorised implementation.

The only part I can't quite piece together is the construction of the cost table, which is done by reversing \(t\). Given that the final result for \(p_{ij}\) ​ is located in the bottom-right corner and we use foldr, I would have expected \(s\) to be the one reversed instead. However, both approaches work, as demonstrated by the following code:

_l ← {¯1⊑(1⊸+⥊+)○≠(⌊`⊢⌊⊏⊸»∘⊢-0∾1+⊣)˝𝔽}
T ← ⌽⊸(=⌜)_l≡=⌜⟜⌽_l
T○{@+97+𝕩•rand.Range 25}´ 1e4‿1e5

1

I suspect the above can be explained by the following properties of the Levenshtein distance:

• \(L(s,t) = L(t,s)\)
• \(L(s,t) = L(\text{rev}(s),\text{rev}(t))\)
• \(L(\text{rev}(s),t) = L(s,\text{rev}(t))\)

If you know why both formulations work, please let me know!

This function computes the real roots of an arbitrary cubic equation. Initially, the equation is transformed into its depressed form via an appropriate substitution. Depending on the sign of the discriminant, the roots are then determined using Cardano's method when the discriminant is positive, or Viète’s trigonometric method when it is negative. In the case where the discriminant is zero, the proportionality to the square of the Vandermonde polynomial implies that a repeated root is present, the roots are resolved through direct analytical methods. We have chosen those methods to avoid using complex numbers, which are not yet supported in BQN.

Cub ← {a‿b‿c‿d:
  (b÷3×a)-˜•math{
    𝕩>0 ? +´𝕩(𝕗.Cbrt+⋈-)⟜√˜-q÷2;
    𝕩=0 ? 0⊸=◶⟨¯1‿2‿2÷˜·𝕗.Cbrt×⟜4,3⊸⥊⟩q;
    (2×√-p÷3)×𝕗.Cos(2×π×↕⊸÷3)-˜3÷˜𝕗.Acos(√-3÷p)×1.5×q÷p
  }(27÷˜p⋆3)+4÷˜×˜q←(d÷a)-(27÷˜3⋆˜b÷a)+3÷˜b×a÷˜p←(c÷a)-3÷˜×˜b÷a
}

The above implementation only works for the case where a≢0, it will yield NaN otherwise. Here are some tests for the four possible branches:

Cub¨ ⟨1‿0‿¯7‿6, 1‿¯1‿¯8‿12, 1‿¯6‿12‿¯8, 1‿3‿0‿¯1⟩ 

⟨ ⟨ 2.0000000000000004 1 ¯3.0000000000000004 ⟩ ⟨ ¯2.9999999999999996 1.9999999999999998 1.9999999999999998 ⟩ ⟨ 2 2 2 ⟩ ⟨ 0.532088886237956 ¯0.6527036446661387 ¯2.879385241571817 ⟩ ⟩

I put some effort golfing this QR decomposition implementation, and I got a very satisfying 98 chars one-liner. Ungolfed a bit, it looks like this:

QR ← +˝∘×⎉1‿∞{
  1=⊢´≢𝕩 ? 𝕩⊸÷⟜⊑⊸⋈√+˝×˜𝕩;
  ∾˘{(q𝔽𝕨)⋈(r𝔽t)∾0𝔽⍟k𝕩}´𝕊𝔽{𝕘-𝕩𝔽t↩𝕩⍉⊸𝔽𝕘}(k↓˘𝕩)⊑q‿r←𝕊𝕩↑˘˜k←⌈2÷˜⊢´≢𝕩⊣t←@
}

The function works like this: it recursively computes the QR decomposition of a matrix by first handling the base case (normalizing a single column) then splitting the matrix into two halves. The first half is decomposed into \(Q_0\) and \(R_0\), and the second half is orthogonalized against \(Q_0\) by subtracting its projection, yielding a residual matrix that is itself decomposed into \(Q_1\) and \(R_1\). Finally, the overall orthogonal matrix \(Q\) is formed by horizontally concatenating \(Q_0\) and \(Q_1\), and the upper triangular \(R\) is assembled as a block matrix combining \(R_0\), the projection coefficients, and \(R_1\):

We can test it with random matrices:

(⊢∾⟜<m-+˝∘×⎉1‿∞´) QR m ← 3‿3•rand.Range 0

┌─                                                                                                                                        
· ┌─                                                             ┌─                                                            ┌─         
  ╵  0.8157427013276365 ¯0.577946856084976 0.02326535562123689   ╵ 0.9106163258394209 0.7411115590785274  0.7652096291273813   ╵ 0 0 0    
    0.32843727859545113 0.4297133155667652 ¯0.8411155809122974                      0  0.709988720748101 0.15322713799622295     0 0 0    
      0.476122672490509 0.6937751061879561  0.5403547934222346                      0                  0 0.36577814222564664     0 0 0    
                                                               ┘                                                             ┘         ┘  
                                                                                                                                         ┘

Below is an implementation of the radix-2 Cooley–Tukey FFT algorithm. The function leverages BQN's headers to define the inverse transform in a succinct way using the property:

We also define a namespace for dealing with complex numbers, in particular the Cis function:

z ← {
  _p ⇐ {(-´𝔽¨)⋈(+´𝔽¨)⟜⌽}
  C‿E ⇐ ⟨⋈⟜-´˘, •math{𝕗.Cos≍˘𝕗.Sin}⟩
}
FFT ← {𝕊⁼: z.C{≠÷˜·𝔽𝔾∘𝔽}𝕊𝕩; (1=≠)◶⟨(+∾-)⟜(⊢×z._p˘·z.E∘-π×↕⊸÷∘≠)´(𝕊¨⊢⊔˜2|⊒˜), ⊢⟩𝕩}

Let's confirm that the inverse returns back the original list:

(+´∘⥊⊢-FFT⁼∘FFT) 0•rand.Range˜2⋈˜2⋆10

1.914614300435602e¯14

We could also compare with the discrete Fourier transform, which despite being \(O(N^2)\) has should have a nice array formulation:

DFT ← ≍˘´<˘{𝔽∘⍉+˝∘×⎉1‿∞ z._p˜·𝔽1‿0⍉⁼·z.E ¯2×π×≠÷˜·×⌜˜⊒˜}
(+´∘⥊FFT-DFT) 0•rand.Range˜2⋈˜2⋆10

¯2.8412011632283907e¯10

In the DFT code above, I got into a big mess with the complex numbers, because the z namespace was too tightly coupled with the FFT implementation. I had to do a bunch of enclosing and coupling to get the same shape. With proper complex numbers support it would be something like:

DFT ← ⊢+˝∘×⎉1‿∞˜·⋆¯2×π×≠÷˜·×⌜˜⊒˜

The n-mode product is a key ingredient for computing the Tucker decomposition of a tensor. For this we can use the HOSVD algorithm: a method that has been rediscovered several times. For example, in the nuclear quantum dynamics community it is known as POTFIT and was published before the often cited De Lathauwer paper, see arXiv:1309.5060 for a discussion. For a tensor \(\mathcal{X}\) and a matrix \(U\) we define:

In BQN's parlance, we can express it as:

{+˝∘×⎉1‿∞⟜𝕩⌾(⍉⍟𝕗)𝕨}

A beautiful example of notation as a tool of thought, in my opinion: this deferred 1-modifier (itself a compact melange of six modifiers) computes the 𝕗-mode product of a tensor 𝕨 and a matrix 𝕩. It works by moving the 𝕗-axis to the front, then multiplying 𝕨 and 𝕩 without the need for explicit unfolding, courtesy of the rank operator, and moving the last axis of the result back to 𝕗, all gracefully managed by under.