Array languages like APL are famous for using many pretty squiggles for\nrepresenting built-in functions. These are called "primitives" and are\ndifferent from user-defined functions, because the interpreters know about them\nand can do fancy optimizations with them, like recognizing special combinations\n[1].
\nAPL programmers have come to rely on these optimizations, for example the\nvectorization of some combinations like +/ (sum). if you were to feed your\ntypical APL interpreter the function {⍺+⍵}/ - which is functionally\nequivalent to +/ [2] since {⍺+⍵} is\njust a dfn [3] over + - the\ninterpreter would execute it much slower, not just because executing the\nabstraction does two function calls per element (one for the dfn, one for the\nprimitive +) but because it does any at all.
Any respectable APL has a tight loop implemented in native code dedicated to\nthis exact combination, which is why APL can sometimes compete with languages\nthat compile to native code.
\nTable of Contents
\nfold is a higher order function generalizing the concept of iterating over
an ordered collection repeatedly calling a binary function that takes as\ninputs an accumulator and an element of the collection, and returns a new\naccumulator.
scan is the same, but instead of returning the final value of the
accumulator it returns a collection of all intermediate values the\naccumulator took over the course of the iteration.
In both cases an initial value may or may not be provided, and in the second\ncase the first element is treated as an initial value, and is returned\nunchanged, concatenated with the scan of the rest of the collection.
\nAPL represents these functions with the operators / and \\, so that
10 = +/ 1 2 3 4\n 1 3 6 10 = +\\ 1 2 3 4In the rest of this article, I’ll be using a fantasy dialect of APL modified to\nhave better behaved fold and scan semantics [4] to write some functional equations. If you don’t know APL,\ndon’t worry about it, just pretend this is the scan from any traditional\nfunctional language like Lisp, Haskell, Ocaml, etc. The main equations will\nalso be translated to Haskell.
\nWhen I say APLers rely on these, I really mean it! scans, folds, reverse scans,\nwith the operators min, max, plus, times are all familiar, if not even routine\noperations for the average python programmer who is comfortable with numpy,\nbecause they really are useful.
\nFolding (from the right) with subtraction gives an alternating sum, that’s cute.
\nAPL programs go beyond that. What does scanning with the less-than function do?\nscanning with not-equal? folding with… modulo? [5]
\nThe issue here is that APL needs escape hatches to handle data dependencies in\nalgorithms, both because it’s syntactically awkward to write loops and because\nexplicitly iterating over an array in APL is mind-numbingly slow when compared\nto the speed at which the primitive functions run.\nScans and folds are a very general way to accumulate information and "move it\naround" in the array, and they’re APL’s most idiomatic way to do so.
\nSay you have a string of text containing quote characters, and you want to tell\nwhich positions in the string are inside a pair of quotes or not: if you were to\nwrite that in C you’d do something like this:
\nbool in_quotes = false;\nfor (size_t i = 0; i < n; ++i) {\n if (string[i] == '"') {\n in_quotes = !in_quotes;\n }\n}While in APL you would write the following, where the result is an array containing the values of\nin_quotes from the previous example. The idea is to scan over the boolean results from comparing\neach character in the string with the " character.
in_quotes ← ≠\\string='"'To see why this works, let’s make the C version branchless by encoding the if statement\ninto arithmetic; to make the similarity more apparent we can write != where ordinarily\n^ would be more natural, on the values 0 and 1 they do the same thing.
bool in_quotes = false;\nfor (size_t i = 0; i < n; ++i) {\n in_quotes = in_quotes != (string[i] == '"');\n}APL has a lot of primitives. it’s the only language I know of which has all non-trivial binary\nboolean functions built in: here they are, next to their truth tables.
\n| Common name | 00 | 01 | 10 | 11 | |
|---|---|---|---|---|---|
∧ | and | 0 | 0 | 0 | 1 |
> | greater than | 0 | 0 | 1 | 0 |
< | less than | 0 | 1 | 0 | 0 |
≠ | not equal (xor) | 0 | 1 | 1 | 0 |
∨ | or | 0 | 1 | 1 | 1 |
⍱ | nor | 1 | 0 | 0 | 0 |
= | equal | 1 | 0 | 0 | 1 |
≥ | greater than or equal | 1 | 0 | 1 | 1 |
≤ | less than or equal (implication) | 1 | 1 | 0 | 1 |
⍲ | nand | 1 | 1 | 1 | 0 |
This is an issue on an implementer’s side, because many primitives mean many\ncombinations, and writing decent code for each of them costs in labour as well\nas executable size (therefore instruction cache). Moreover, what qualifies as\ndecent code for a boolean scan is a much higher bar than the snippet I showed\nabove: APL uses packed bit vectors to represent boolean arrays, meaning that\nthese operations have the potential to consume tens of elements per CPU cycle,\nwhile by a quick benchmark a naive implementation would handle about 0.6\nelements per cycle [6]
\nTelling which of these are really useful only looking at the truth tables is\nnot trivial, and it really sucks as a user to learn that the operation\nthat precisely solves your problem is shoddily implemented.
\nIs there some way that we can reduce the implementation work without sacrificing\nthe performance?
\nThe structure of the scan ⍺ f\\ ⍵ is as follows:
Suppose we want to do some pre-processing with a function g to alter the behaviour of the scan:\nwhat we’d end up with is a data flow of this kind:
The motivation for looking at this is that preprocessing a boolean array with a\nunary scalar function is just about the cheapest thing we could do, in fact there\nis only one non-trivial boolean unary scalar function: boolean negation.\nLet’s keep this observation aside for now, and ask ourselves: is this operation still\na scan? If so, what function operand should scan take to produce this result?
\nIt’s apparent that the right argument to f is being fed through g at each step in the\niteration, so the answer is {⍺ f (g ⍵)}, which brings us to the simple equation
⍺ f\\ g¨ ⍵ ←→ ⍺ (f∘g)\\ ⍵Equivalent Haskell
\nscanl f x (map g ys) = scanl (\\x y -> f x (g y)) x ysIf g is any function, where ∘ is APL’s composition operator "jot" that satisfies
⍺ f∘g ⍵ ←→ ⍺ f (g ⍵)It’s a little trickier to see what effect post-processing with g would have:
It may look exactly like before, but notice that the input to each subsequent f call\nis not the same as the value that is sent to the output array, so this is not a scan!\nTo fix this, let’s assume that g has a left inverse g-1 and insert it in between each pair of calls:
This is a scan diagram using the function {g ((g⍣¯1) ⍺) f ⍵}, except the initial value x has been preprocessed with g.\nthis gives the equation:
g¨ ⍺ f\\ ⍵ ←→ (g ⍺) (g⍤(g⍣¯1⍛f))\\ ⍵Equivalent Haskell
\nmap g (scanl f x ys) = scanl (\\x y -> g (f (ig x) y)) (g x) ysAgain some composition operators have been used:
\n(g⍣¯1) g ⍵ ←→ ⍵ ⍝ APL can sometimes find inverses by itself!\n⍺ (g⍤f) ⍵ ←→ g (⍺ f ⍵)\n⍺ (g⍛f) ⍵ ←→ (g ⍺) f ⍵And finally we can put these two relations together to find a third.
\ng ⍺ f\\ g ⍵ ←→ (g ⍺) g⍤(g⍣¯1⍛f∘g)\\ ⍵Equivalent Haskell
\nmap g (scanl f x (map g ys)) = scanl (\\x y -> g (f (ig x) (g y))) (g x) ysIn particular, we are interested in seeing which pairs of primitives are\nrelated by one of these three relations, where g is the negation function ~.\nWe start with the truth tables of our primitives and analyze how they change under the\ntransformations described above.
We represent all of the truth tables at once as the array a with shape 10 2 2
a\n┌┌→──┐\n↓↓0 0│\n││0 1│\n││ │\n││0 0│\n││1 0│\n││ │\n││0 1│\n││0 0│\n││ │\n││0 1│\n││1 0│\n││ │\n││0 1│\n││1 1│\n││ │\n││1 0│\n││0 0│\n││ │\n││1 0│\n││0 1│\n││ │\n││1 0│\n││1 1│\n││ │\n││1 1│\n││0 1│\n││ │\n││1 1│\n││1 0│\n└└~──┘Composing on the right argument with ~ is the same as reversing the last axis with\nthe ⌽ function:
pre ← ⌽[2] aNegating the output is just a matter of negating each entry of the table; if we do this, we must also\nnegate the left argument, which we do by reversing the first axis of each truth table (the middle axis here).
\n post ← ~ ⌽[1] aAnd of course we can do both transformations at once
\n both ← ~ ⌽[1] ⌽[2] aNow we can find out who’s related to whom by searching for the truth tables we got in the original array:
\n index ← ↑(⊂a)⍳¨a pre post both┌→──────────────────┐\n↓0 1 2 3 4 5 6 7 8 9│\n│1 0 5 6 7 2 3 4 9 8│\n│7 4 9 3 1 8 6 0 5 2│\n│4 7 8 6 0 9 3 1 2 5│\n└~──────────────────┘Now that we know this, we can index into the list of the symbols the truth tables referred to:
\n '∧><≠∨⍱=≥≤⍲'[index]┌→─────────┐\n↓∧><≠∨⍱=≥≤⍲│\n│>∧⍱=≥<≠∨⍲≤│\n│≥∨⍲≠>≤=∧⍱<│\n│∨≥≤=∧⍲≠><⍱│\n└──────────┘Now, how do we interpret this? along the top row are the\noriginal symbols; the second row are the symbols that\nare related to the top row by the first scan equation we derived\nabove, i.e. ⍺ p\\ ~⍵ ←→ ⍺ q\\ ⍵.\nThe third and fourth row, similarly, are related to the first\nvia ~ (~⍺) p\\ ⍵ ←→ ⍺ q\\ ⍵ and ~ (~⍺) p\\ ~⍵ ←→ ⍺ q\\ ⍵
We can see that the set ∧, <, ≠ covers all columns at\nleast once, therefore if we provide fast implementations of each of\nthese, we can derive decent implementations of the others by just\napplying these symmetries. This cuts down on implementation time\nand executable size, without completely sacrificing performance.
None of this is completely new, by the way. The trick ⍺ ≠\\ ⍵ ←→ ⍺ =\\ ~⍵\nis well known among K programmers, since K doesn’t have a builtin\nnot-equal operator. In the same vein, many of these shortcuts\nare employed in the CBQN interpreter for BQN for the\nless common\ncombinations.
I originally discussed this on the APL farm discord server, where some of the users\nmade interesting comments:
\nIn the K language, the form x f\\ y, is equivalent to the x a\\ y form, where\na is a matrix containing the values of f on pairs of non-negative integers.\nThis form advances the state machine\n[7]\nwith transitions defined in a over the sequence y. In this interpretation, the concrete\ntransformations we applied to the truth tables take the place of the\ncomposition operators in the original equations we defined.
( x a\\ ~y) = x ( |'a)\\ y\n(~ x a\\ y) = (~x) (~| a)\\ y\n(~ x a\\ ~y) = (~x) (~||'a)\\ yEquivalent Haskell
\nrunFsm a x ys = scan (\\x y -> (a !! x) !! y) x ys\n\n runFsm a x $ map (1 -) ys = runFsm x ( map reverse a) ys\nmap (1 -) $ runFsm a x ys = runFsm (1 - x) (map (map (1 -)) $ reverse a) ys\nmap (1 -) $ runFsm a x $ map (1 -) ys = runFsm (1 - x) (map (map (1 -)) $ map reverse $ reverse a) ysThe transformations we gave permute the rows or columns of the tables\nand/or negate the entries. This implies that the number of ones in the tables either remains\nequal to 2 or flip-flops between 1 and 3 [8]. Moreover, since the transformations commute and are self-inverses,\nour relation is an equivalence relation; we deduce that = and ≠ are in a class of their\nown, with one transformation fixing both operations, and the others fill at least two more.
On these, the transformation ~⌽[0]a can be applied so that there are 3 0s,\nthen ⌽[1]a can be applied to ensure that the single 1 entry occurs on the second column.\nThe only two possibilities are distinguished by whether the 1 value is at\na[1;1] (which gives ∧) or at a[0;1] (which gives <)
A xor scan on a machine word can be implemented efficiently in a few instructions by a typical prefix doubling\ntechnique, similar to the one described by Bit Twiddling Hacks
\nuint64_t xorscan_word (uint64_t word) {\n word ^= word << 1;\n word ^= word << 2;\n word ^= word << 4;\n word ^= word << 8;\n word ^= word << 16;\n word ^= word << 32;\n return word;\n}We only need to handle the carry from the previous word. The typical way to do this would be\nto xor it with the first bit in word, but this is actually sub-optimal: notice that carry\nis needed at the very start of each iteration and only updated at the very end, which ties all iterations\nin a very long dependency chain.
void xorscan_bad (bool initial, uint64_t *words, size_t n) {\n uint64_t carry = (uint64_t)initial;\n for (size_t i = 0; i < n; ++i) {\n uint64_t word = xorscan_word(words[i] ^ carry);\n carry = word >> 63;\n words[i] = word;\n }\n}This compiles to the following core loop, which runs at around 3-4 bits per cycle.\n(all the timings I’ll give are for memory that is hot in cache, otherwise memory latency\nwould dominate runtime)
\nAssembly output
\n.loop:\n xor (%rdx),%rdi\n add $0x8,%rdx\n lea (%rdi,%rdi,1),%rax\n xor %rdi,%rax\n lea 0x0(,%rax,4),%rcx\n xor %rcx,%rax\n mov %rax,%rcx\n shl $0x4,%rcx\n xor %rcx,%rax\n mov %rax,%rcx\n shl $0x8,%rcx\n xor %rcx,%rax\n mov %rax,%rcx\n shl $0x10,%rcx\n xor %rcx,%rax\n mov %rax,%rcx\n shl $0x20,%rcx\n xor %rcx,%rax\n mov %rax,%rdi\n mov %rax,-0x8(%rdx)\n shr $0x3f,%rdi\n cmp %rdx,%rsi\n jne .loopOn a CPU with super-scalar execution, it would be nice if we could delay the usage of the previous result until\nthe latest possible time.\nFortunately xor is associative, so we can just broadcast the carry bit to every position in the final result: now\nthe CPU is free to start executing xorscan_word for the next iteration before the result of the previous one is ready.
void xorscan_good (bool initial, uint64_t *words, size_t n) {\n uint64_t carry = - (uint64_t)initial;\n for (size_t i = 0; i < n; ++i) {\n uint64_t word = carry ^ xorscan_word(words[i]);\n carry = - (word >> 63);\n words[i] = word;\n }\n}This is much faster, processing up to 10.5 bits per cycle.
\nAssembly output
\n.loop:\n mov (%rcx),%rdx\n add $0x8,%rcx\n lea (%rdx,%rdx,1),%rax\n xor %rdx,%rax\n lea 0x0(,%rax,4),%rdx\n xor %rdx,%rax\n mov %rax,%rdx\n shl $0x4,%rdx\n xor %rdx,%rax\n mov %rax,%rdx\n shl $0x8,%rdx\n xor %rdx,%rax\n mov %rax,%rdx\n shl $0x10,%rdx\n xor %rdx,%rax\n mov %rax,%rdx\n shl $0x20,%rdx\n xor %rdx,%rax\n xor %rsi,%rax\n mov %rax,%rsi\n mov %rax,-0x8(%rcx)\n sar $0x3f,%rsi\n cmp %rcx,%rdi\n jne .loopIf the CPU supports a carry-less multiply instruction (many modern CPUs do!) then xorscan_word can be implemented as carry-less multiplication with the word 0xffffffffffffffff (all ones).
All of the above (except for the carry-less multiply) can be applied verbatim, by replacing ^ with &; however there is a faster\nmethod available.
Firstly, a single-word and scan can be implemented much more simply than xor\n[9].
\nuint64_t andscan_word (uint64_t word) {\n return word & ~(word + 1);\n}This works because adding one turns the first run of ones into a run of\nzeros, then the negation turns the initial run of zeros back into a run of ones,\nas well as negating everything else. The final bitwise and\nsees two matching runs of ones, a matching pair of zero bits, and mismatched\nbits everywhere else.
\nHere is an example, with the bits ordered from least to most significant.
\nword | 11100011010101011 |
word+1 | 00010011010101011 |
~(word + 1) | 11101100101010100 |
word & ~(word + 1) | 11100000000000000 |
To extend this to work on an array of words, we observe that at most one word\nwill need an actual scan, because everything before it is known to be all ones\nalready, and everything after it can be unconditionally set to zero.
\nThe full implementation might be something like the following.
\nvoid andscan (bool initial, uint64_t *words, size_t n) {\n size_t i = 0;\n if (initial) {\n for (; i < n; ++i) {\n uint64_t word = words[i];\n if (~word) {\n words[i] = andscan_word(word);\n ++i;\n break;\n }\n }\n }\n for (; i < n; ++i) {\n words[i] = 0;\n }\n}This produces two separate loops, the first and most intensive of them can do 62 bits per cycle, while the\nother gets compiled into a call to memset.
\nAssembly output
\n.loop:\n add $0x8,%rcx\n cmp %rax,%rdx\n je .found\n mov (%rcx),%rsi\n inc %rax\n cmp $0xffffffffffffffff,%rsi\n je .loopThis is the most complicated case, and the solution I’ll present comes from the\nsimdjson\nlibrary, which uses it to find escaped characters. In fact, the intuition behind\nwhat <\\ does is that it turns off every 1 bit that is at an odd offset in\na contiguous sequence of 1s, which is exactly what you need to figure out\nwhich backslashes in a string work as escapes and which are escaped.
uint64_t ltscan_word (uint64_t word) {\n uint64_t guard = word << 1 | 0xaaaaaaaaaaaaaaaa;\n uint64_t parity = guard - word;\n return (parity ^ 0xaaaaaaaaaaaaaaaa) & word;\n}Here is an example, again in little endian bit order.
\nword | 1110000111100111 |
word << 1 | 0111000011110011 |
guard | 0111010111110111 |
parity | 1110010000010000 |
parity ^ odd | 1011000101000101 |
result | 1010000101000101 |
Shifting and setting the odd bits controls whether or not a borrow chain is\ninitiated in the subtraction: the start of a run of 1s is set to 0 if it\nis at an even position, otherwise it stays 1. The last element of a run is\nalso moved past the end to stop the borrow chain. Then subtracting the original\nbit mask creates a chain of borrows which leaves every bit set in the first\ncase, and just leaves every bit unset in the second case. If we now flip every\nodd bit, we will have the even bits set in the runs that start at even\npositions, the odd bits set in the runs that start at odd positions, and\nsome garbage outside of these runs which we clean up with bitwise and.
Like the xor case, the trivial way to chain these together is to clear the\nfirst bit whenever the last bit of the previous result is set, though\nMarshall’s\nnotes\non this also mention a smarter way to reduce the length of the critical path.
\nvoid ltscan (bool init, uint64_t *words, size_t n) {\n uint64_t carry = (uint64_t)init;\n for (size_t i = 0; i < n; ++i) {\n uint64_t word = ltscan_word(words[i] &~ carry);\n carry = word >> 63;\n words[i] = word;\n }\n}This particular version can do 8.8 bits per cycle.
\nAssembly output
\n.loop:\n andn (%rsi),%rdi,%rdi\n add $0x8,%rsi\n lea (%rdi,%rdi,1),%rax\n or %rdx,%rax\n sub %rdi,%rax\n xor %rdx,%rax\n and %rdi,%rax\n mov %rax,%rdi\n mov %rax,-0x8(%rsi)\n shr $0x3f,%rdi\n cmp %rsi,%rcx\n jne .loop