Why is matrix multiplication defined so very differently from matrix addition? If we didn’t know these procedures, could we derive them from first principles? What might those principles be?

This post gives a simple semantic model for matrices and then uses it to systematically derive the implementations that we call matrix addition and multiplication. The development illustrates what I call “denotational design”, particularly with type class morphisms. On the way, I give a somewhat unusual formulation of matrices and accompanying definition of matrix “multiplication”.

For more details, see the linear-map-gadt source code.

Edits:

• 2012–12–17: Replaced lost $B$ entries in description of matrix addition. Thanks to Travis Cardwell.
• 2012–12018: Added note about math/browser compatibility.

Note: I’m using MathML for the math below, which appears to work well on Firefox but on neither Safari nor Chrome. I use Pandoc to generate the HTML+MathML from markdown+lhs+LaTeX. There’s probably a workaround using different Pandoc settings and requiring some tweaks to my WordPress installation. If anyone knows how (especially the WordPress end), I’d appreciate some pointers.

Matrices

For now, I’ll write matrices in the usual form: $$\left(\begin{array}{ccc}{A}_{11}& \cdots & A_{1m}\\ ⋮& \ddots & ⋮\\ A_{n1}& \cdots & A_{nm}\end{array}\right)$$

Addition

To add two matrices, we add their corresponding components. If $$A=\left(\begin{array}{ccc}{A}_{11}& \cdots & A_{1m}\\ ⋮& \ddots & ⋮\\ A_{n1}& \cdots & A_{nm}\end{array}\right)\mathrm{\text{and}}B=\left(\begin{array}{ccc}{B}_{11}& \cdots & B_{1m}\\ ⋮& \ddots & ⋮\\ B_{n1}& \cdots & B_{nm}\end{array}\right),$$ then $$A+B=\left(\begin{array}{ccc}{A}_{11}+B_{11}& \cdots & A_{1m}+B_{1m}\\ ⋮& \ddots & ⋮\\ A_{n1}+B_{n1}& \cdots & A_{nm}+B_{nm}\end{array}\right).$$ More succinctly, $$(A+B{)}_{ij}=A_{ij}+B_{ij}.$$

Multiplication

Multiplication, on the other hand, works quite differently. If $$A=\left(\begin{array}{ccc}{A}_{11}& \cdots & A_{1m}\\ ⋮& \ddots & ⋮\\ A_{n1}& \cdots & A_{nm}\end{array}\right)\mathrm{\text{and}}B=\left(\begin{array}{ccc}{B}_{11}& \cdots & B_{1p}\\ ⋮& \ddots & ⋮\\ B_{m1}& \cdots & B_{mp}\end{array}\right),$$ then $$(A\bullet B{)}_{ij}=\sum _{k=1}^m{A}_{ik}\cdot B_{kj}.$$ This time, we form the dot product of each $A$ row and $B$ column.

Why are these two matrix operations defined so differently? Perhaps these two operations are implementations of more fundamental specifications. If so, then making those specifications explicit could lead us to clear and compelling explanations of matrix addition and multiplication.

Transforming vectors

Simplifying from matrix multiplication, we have transformation of a vector by a matrix. If $$A=\left(\begin{array}{ccc}{A}_{11}& \cdots & A_{1m}\\ ⋮& \ddots & ⋮\\ A_{n1}& \cdots & A_{nm}\end{array}\right)\mathrm{\text{and}}x=\left(\begin{array}{c}{x}_1\\ ⋮\\ x_m\end{array}\right),$$ then $$A\cdot x=\left(\begin{array}{ccccc}{A}_{11}\cdot x_1& +& \cdots & +& A_{1m}\cdot x_m\\ ⋮& & \ddots & & ⋮\\ A_{n1}\cdot x_1& +& \cdots & +& A_{nm}\cdot x_m\end{array}\right)$$ More succinctly, $$(A\cdot x{)}_i=\sum _{k=1}^m{A}_{ik}\cdot x_k.$$

What’s it all about?

We can interpret matrices as transformations. Matrix addition then adds transformations:

$$(A+B)x=Ax+Bx$$

Matrix “multiplication” composes transformations:

$$(A\bullet B)x=A(Bx)$$

What kinds of transformations?

Linear transformations

Matrices represent linear transformations. To say that a transformation (or “function” or “map”) $f$ is “linear” means that $f$ preserves the structure of addition and scalar multiplication. In other words, $$\begin{array}{ccc}\hfill f(x+y)& \hfill =\hfill & fx+fy\hfill \\ \hfill f(c\cdot x)& \hfill =\hfill & c\cdot fx\hfill \end{array}$$ Equivalently, $f$ preserves all linear combinations: $$f(c_1\cdot x_1+\cdots +c_m\cdot x_m)=c_1\cdot f{x}_1+\cdots +c_m\cdot f{x}_m$$

What does it mean to say that “matrices represent linear transformations”? As we saw in the previous section, we can use a matrix to transform a vector. Our semantic function will exactly be this use, i.e., the meaning of matrix is as a function (map) from vectors to vectors. Moreover, these functions will satisfy the linearity properties above.

Representation

For simplicity, I’m going structure matrices in a unconventional way. Instead of a rectangular arrangement of numbers, use the following generalized algebraic data type (GADT):

data a ⊸ b where
  Dot   ∷ InnerSpace b ⇒
          b → (b ⊸ Scalar b)
  (:&&) ∷ VS3 a c d ⇒  -- vector spaces with same scalar field
          (a ⊸ c) → (a ⊸ d) → (a ⊸ c × d)

I’m using the notation “c × d” in place of the usual “(c,d)”. Precedences are such that “×” binds more tightly than “⊸”, which binds more tightly than “→”.

This definition builds on the VectorSpace class, with its associated Scalar type and InnerSpace subclass. Using VectorSpace is overkill for linear maps. It suffices to use modules over semirings, which means that we don’t assume multiplicative or additive inverses. The more general setting enables many more useful applications than vector spaces do, some of which I will describe in future posts.

The idea here is that a linear map results in either (a) a scalar, in which case it’s equivalent to dot v (partially applied dot product) for some v, or (b) a product, in which case it can be decomposed into two linear maps with simpler range types. Each row in a conventional matrix corresponds to Dot v for some vector v, and the stacking of rows corresponds to nested applications of (:&&).

Semantics

The semantic function, apply, interprets a representation of a linear map as a function (satisfying linearity):

apply ∷ (a ⊸ b) → (a → b)
apply (Dot b)   = dot b
apply (f :&& g) = apply f &&& apply g

where, (&&&) is from Control.Arrow.

(&&&) ∷ Arrow (↝) ⇒ (a ↝ b) → (a ↝ c) → (a ↝ (b,c))

For functions,

(f &&& g) a = (f a, g a)

Functions, linearity, and multilinearity

Functions form a vector space, with scaling and addition defined “pointwise”. Instances from the vector-space package:

instance AdditiveGroup v ⇒ AdditiveGroup (a → v) where
  zeroV   = pure   zeroV
  (^+^)   = liftA2 (^+^)
  negateV = fmap   negateV

instance VectorSpace v ⇒ VectorSpace (a → v) where
  type Scalar (a → v) = a → Scalar v
  (*^) s = fmap (s *^)

I wrote the definitions in this form to fit a template for applicative functors in general. Inlining the definitions of pure, liftA2, and fmap on functions, we get the following equivalent instances:

instance AdditiveGroup v ⇒ AdditiveGroup (a → v) where
  zeroV     = λ _ → zeroV
  f ^+^ g   = λ a → f a ^+^ g a
  negateV f = λ a → negateV (f a)

instance VectorSpace v ⇒ VectorSpace (a → v) where
  type Scalar (a → v) = a → Scalar v
  s *^ f = λ a → s *^ f a

In math, we usually say that dot product is “bilinear”, or “linear in each argument”, i.e.,

dot (s *^ u,v) ≡ s *^ dot (u,v)
dot (u ^+^ w, v) ≡ dot (u,v) ^+^ dot (w,v)

Similarly for the second argument:

dot (u,s *^ v) ≡ s *^ dot (u,v)
dot (u, v ^+^ w) ≡ dot (u,v) ^+^ dot (u,w)

Now recast the first of these properties in a curried form:

dot (s *^ u) v ≡ s *^ dot u v

i.e.,

dot (s *^ u)
 ≡ {- η-expand -}
λ v → dot (s *^ u) v
 ≡ {- "bilinearity" -}
λ v → s *^ dot u v
 ≡ {- (*^) on functions -}
λ v → (s *^ dot u) v
 ≡ {- η-contract -}
s *^ dot u

Likewise,

dot (u ^+^ v)
 ≡ {- η-expand -}
λ w → dot (u ^+^ v) w
 ≡ {- "bilinearity" -}
λ w → dot u w ^+^ dot v w
 ≡ {- (^+^) on functions -}
dot u ^+^ dot v

Thus, when “bilinearity” is recast in terms of curried functions, it becomes just linearity. (The same reasoning applies more generally to multilinearity.)

Note that we could also define function addition as follows:

f ^+^ g = add ∘ (f &&& g)

where

add = uncurry (^+^)

This uncurried form will come in handy in derivations below.

Deriving matrix operations

Addition

We’ll add two linear maps using the (^+^) operation from Data.AdditiveGroup.

(^+^) ∷ (a ⊸ b) → (a ⊸ b) → (a ⊸ b)

Following the principle of semantic type class morphisms, the specification simply says that the meaning of the sum is the sum of the meanings:

apply (f ^+^ g) ≡ apply f ^+^ apply g

which is half of the definition of “linearity” for apply.

The game plan (as always) is to use the semantic specification to derive (or “calculate”) a correct implementation of each operation. For addition, this goal means we want to come up with a definition like

f ^+^ g = <rhs>

where <rhs> is some expression in terms of f and g whose meaning is the same as the meaning as f ^+^ g, i.e., where

apply (f ^+^ g) ≡ apply <rhs>

Since Haskell has convenient pattern matching, we’ll use it for our definition of (^+^) above. Addition has two arguments, and our data type has two constructors, there are at most four different cases to consider.

First, add Dot and Dot. The specification

apply (f ^+^ g) ≡ apply f ^+^ apply g

specializes to

apply (Dot b ^+^ Dot c) ≡ apply (Dot b) ^+^ apply (Dot c)

Now simplify the right-hand side (RHS):

apply (Dot b) ^+^ apply (Dot c)
 ≡ {- apply definition -}
dot b ^+^ dot c
 ≡ {- (bi)linearity of dot, as described above -}
dot (b ^+^ c)
 ≡ {- apply definition -}
apply (Dot (b ^+^ c))

So our specialized specification becomes

apply (Dot b ^+^ Dot c) ≡ apply (Dot (b ^+^ c))

which is implied by

Dot b ^+^ Dot c ≡ Dot (b ^+^ c)

and easily satisfied by the following partial definition (replacing “≡” by “=”):

Dot b ^+^ Dot c = Dot (b ^+^ c)

Now consider the case of addition with two (:&&) constructors:

The specification specializes to

apply ((f :&& g) ^+^ (h :&& k)) ≡ apply (f :&& g) ^+^ apply (h :&& k)

As with Dot, simplify the RHS:

apply (f :&& g) ^+^ apply (h :&& k)
 ≡ {- apply definition -}
(apply f &&& apply g) ^+^ (apply h &&& apply k)
 ≡ {- See below -}
(apply f ^+^ apply h) &&& (apply g ^+^ apply k)
 ≡ {- induction -}
apply (f ^+^ h) &&& apply (g ^+^ k)
 ≡ {- apply definition -}
apply ((f ^+^ h) :&& (g ^+^ k))

I used the following property (on functions):

(f &&& g) ^+^ (h &&& k) ≡ (f ^+^ h) &&& (g ^+^ k)

Proof:

(f &&& g) ^+^ (h &&& k)
 ≡ {- η-expand -}
λ x → ((f &&& g) ^+^ (h &&& k)) x
 ≡ {- (&&&) definition for functions -}
λ x → (f x, g x) ^+^ (h x, k x)
 ≡ {- (^+^) definition for pairs -}
λ x → (f x ^+^ h x, g x ^+^ k x)
 ≡ {- (^+^) definition for functions -}
λ x → ((f ^+^ h) x, (g ^+^ k) x)
 ≡ {- (&&&) definition for functions -}
(f ^+^ h) &&& (g ^+^ k)

The specification becomes

apply ((f :&& g) ^+^ (h :&& k)) ≡ apply ((f ^+^ h) :&& (g ^+^ k))

which is easily satisfied by the following partial definition

(f :&& g) ^+^ (h :&& k) = (f ^+^ h) :&& (g ^+^ k)

The other two cases are (a) Dot and (:&&), and (b) (:&&) and Dot, but they don’t type-check (assuming that pairs are not scalars).

Composing linear maps

I’ll write linear map composition as “g ∘ f”, with type

(∘) ∷ (b ⊸ c) → (a ⊸ b) → (a ⊸ c)

This notation is thanks to a Category instance, which depends on a generalized Category class that uses the recent ConstraintKinds language extension. (See the source code.)

Following the semantic type class morphism principle again, the specification says that the meaning of the composition is the composition of the meanings:

apply (g ∘ f) ≡ apply g ∘ apply f

In the following, note that the ∘ operator binds more tightly than &&&, so f ∘ h &&& g ∘ h means (f ∘ h) &&& (g ∘ h).

Derivation

Again, since there are two constructors, we have four possible cases cases. We can handle two of these cases together, namely (:&&) and anything. The specification:

apply ((f :&& g) ∘ h) ≡ apply (f :&& g) ∘ apply h

Reasoning proceeds as above, simplifying the RHS of the constructor-specialized specification.

Simplify the RHS:

apply (f :&& g) ∘ apply h
 ≡ {- apply definition -}
(apply f &&& apply g) ∘ apply h
 ≡ {- see below -}
apply f ∘ apply h &&& apply g ∘ apply h
 ≡ {- induction -}
apply (f ∘ h) &&& apply (g ∘ h)
 ≡ {- apply definition -}
apply (f ∘ h :&& g ∘ h)

This simplification uses the following property of functions:

(p &&& q) ∘ r ≡ p ∘ r &&& q ∘ r

Sufficient definition:

(f :&& g) ∘ h = f ∘ h :&& g ∘ h

We have two more cases, specified as follows:

apply (Dot c ∘ Dot b) ≡ apply (Dot c) ∘ apply (Dot b)

apply (Dot c ∘ (f :&& g)) ≡ apply (Dot c) ∘ apply (f :&& g)

Based on types, c must be a scalar in the first case and a pair in the second. (Dot b produces a scalar, while f :&& g produces a pair.) Thus, we can write these two cases more specifically:

apply (Dot s ∘ Dot b) ≡ apply (Dot s) ∘ apply (Dot b)

apply (Dot (a,b) ∘ (f :&& g)) ≡ apply (Dot (a,b)) ∘ apply (f :&& g)

In the derivation, I won’t spell out as many details as before. Simplify the RHSs:

  apply (Dot s) ∘ apply (Dot b)
≡ dot s ∘ dot b
≡ dot (s *^ b)
≡ apply (Dot (s *^ b))

  apply (Dot (a,b)) ∘ apply (f :&& g)
≡ dot (a,b) ∘ (apply f &&& apply g)
≡ add ∘ (dot a ∘ apply f &&& dot b ∘ apply g)
≡ dot a ∘ apply f ^+^ dot b ∘ apply g
≡ apply (Dot a ∘ f ^+^ Dot b ∘ g)

I’ve used the following properties of functions:

dot (a,b)             ≡ add ∘ (dot a *** dot b)

(r *** s) ∘ (p &&& q) ≡ r ∘ p &&& s ∘ q

add ∘ (p &&& q)       ≡ p ^+^ q

apply (f ^+^ g)       ≡ apply f ^+^ apply g

Implementation:

 Dot s     ∘ Dot b     = Dot (s *^ b)
 Dot (a,b) ∘ (f :&& g) = Dot a ∘ f ^+^ Dot b ∘ g

Cross products

Another Arrow operation handy for linear maps is the parallel composition (product):

(***) ∷ (a ⊸ c) → (b ⊸ d) → (a × b ⊸ c × d)

The specification says that apply distributes over (***). In other words, the meaning of the product is the product of the meanings.

apply (f *** g) = apply f *** apply g

Where, on functions,

p *** q = λ (a,b) → (p a, q b)
        ≡ p ∘ fst &&& q ∘ snd

Simplify the specifications RHS:

  apply f *** apply g
≡ apply f ∘ fst &&& apply g ∘ snd

If we knew how to represent fst and snd via our linear map constructors, we’d be nearly done. Instead, let’s suppose we have the following functions.

compFst ∷ VS3 a b c ⇒ a ⊸ c → a × b ⊸ c
compSnd ∷ VS3 a b c ⇒ b ⊸ c → a × b ⊸ c

specified as follows:

apply (compFst f) ≡ apply f ∘ fst
apply (compSnd g) ≡ apply g ∘ snd

With these two functions (to be defined) in hand, let’s try again.

  apply f *** apply g
≡ apply f ∘ fst &&& apply g ∘ snd
≡ apply (compFst f) &&& apply (compSnd g)
≡ apply (compFst f :&& compSnd g)

Composing with fst and snd

I’ll elide even more of the derivation this time, focusing reasoning on the meanings. Relating to the representation is left as an exercise. The key steps in the derivation:

dot a     ∘ fst ≡ dot (a,0)

(f &&& g) ∘ fst ≡ f ∘ fst &&& g ∘ fst

dot b     ∘ snd ≡ dot (0,b)

(f &&& g) ∘ snd ≡ f ∘ snd &&& g ∘ snd

Implementation:

compFst (Dot a)   = Dot (a,zeroV)
compFst (f :&& g) = compFst f &&& compFst g

compSnd (Dot b)   = Dot (zeroV,b)
compSnd (f :&& g) = compSnd f &&& compSnd g

where zeroV is the zero vector.

Given compFst and compSnd, we can implement fst and snd as linear maps simply as compFst id and compSnd id, where id is the (polymorphic) identity linear map.

Reflections

This post reflects an approach to programming that I apply wherever I’m able. As a summary:

• Look for an elegant what behind a familiar how.
• Define a semantic function for each data type.
• Derive a correct implementation from the semantics.

You can find more examples of this methodology elsewhere in this blog and in the paper Denotational design with type class morphisms.

I’ve been thinking much more about parallel computation for the last couple of years, especially since starting to work at Tabula a year ago. Until getting into parallelism explicitly, I’d naïvely thought that my pure functional programming style was mostly free of sequential bias. After all, functional programming lacks the implicit accidental dependencies imposed by the imperative model. Now, however, I’m coming to see that designing parallel-friendly algorithms takes attention to minimizing the depth of the remaining, explicit data dependencies.

As an example, consider binary addition, carried out from least to most significant bit (as usual). We can immediately compute the first (least significant) bit of the result, but in order to compute the second bit, we’ll have to know whether or not a carry resulted from the first addition. More generally, the $(n+1)$th sum & carry require knowing the $n$th carry, so this algorithm does not allow parallel execution. Even if we have one processor per bit position, only one processor will be able to work at a time, due to the linear chain of dependencies.

One general technique for improving parallelism is speculation—doing more work than might be needed so that we don’t have to wait to find out exactly what will be needed. In this post, we’ll see a progression of definitions for bitwise addition. We’ll start with a linear-depth chain of carry dependencies and end with logarithmic depth. Moreover, by making careful use of abstraction, these versions will be simply different type specializations of a single polymorphic definition with an extremely terse definition.

A full adder

Let’s start with an adder for two one-bit numbers. Because of the possibility of overflow, the result will be two bits, which I’ll call “sum” and “carry”. So that we can chain these one-bit adders, we’ll also add a carry input.

addB ∷ (Bool,Bool) → Bool → (Bool,Bool)

In the result, the first Bool will be the sum, and the second will be the carry. I’ve curried the carry input to make it stand out from the (other) addends.

There are a few ways to define addB in terms of logic operations. I like the following definition, as it shares a little work between sum & carry:

addB (a,b) cin = (axb ≠ cin, anb ∨ (cin ∧ axb))
 where
   axb = a ≠ b
   anb = a ∧ b

I’m using (≠) on Bool for exclusive or.

A ripple carry adder

Now suppose we have not just two bits, but two sequences of bits, interpreted as binary numbers arranged from least to most significant bit. For simplicity, I’d like to assume that these sequences to have the same length, so rather than taking a pair of bit lists, let’s take a list of bit pairs:

add ∷ [(Bool,Bool)] → Bool → ([Bool],Bool)

To implement add, traverse the list of bit pairs, threading the carries:

add [] c     = ([]  , c)
add (p:ps) c = (s:ss, c'')
 where
   (s ,c' ) = addB p c
   (ss,c'') = add ps c'

State

This add definition contains a familiar pattern. The carry values act as a sort of state that gets updated in a linear (non-branching) way. The State monad captures this pattern of computation:

newtype State s a = State (s → (a,s))

By using State and its Monad instance, we can shorten our add definition. First we’ll need a new full adder definition, tweaked for State:

addB ∷ (Bool,Bool) → State Bool Bool
addB (a,b) = do cin ← get
                put (anb ∨ cin ∧ axb)
                return (axb ≠ cin)
 where
   anb = a ∧ b
   axb = a ≠ b

And then the multi-bit adder:

add ∷ [(Bool,Bool)] → State Bool [Bool]
add []     = return []
add (p:ps) = do s  ← addB p
                ss ← add ps
                return (s:ss)

We don’t really need the Monad interface to define add. The simpler and more general Applicative interface suffices:

add []     = pure []
add (p:ps) = liftA2 (:) (addB p) (add ps)

This pattern also looks familiar. Oh — the Traversable instance for lists makes for a very compact definition:

add = traverse addB

Wow. The definition is now so simple that it doesn’t depend on the specific choice of lists. To find out the most general type add can have (with this definition), remove the type signature, turn off the monomorphism restriction, and see what GHCi has to say:

add ∷ Traversable t ⇒ t (Bool,Bool) → State Bool (t Bool)

This constraint is very lenient. Traversable can be derived automatically for all algebraic data types, including nested/non-regular ones.

For instance,

data Tree a = Leaf a | Branch (Tree a) (Tree a)
  deriving (Functor,Foldable,Traversable)

We can now specialize this general add back to lists:

addLS ∷ [(Bool,Bool)] → State Bool [Bool]
addLS = add

We can also specialize for trees:

addTS ∷ Tree (Bool,Bool) → State Bool (Tree Bool)
addTS = add

Or for depth-typed perfect trees (e.g., as described in From tries to trees):

addTnS ∷ IsNat n ⇒
         T n (Bool,Bool) → State Bool (T n Bool)
addTnS = add

Binary trees are often better than lists for parallelism, because they allow quick recursive splitting and joining. In the case of ripple adders, we don’t really get parallelism, however, because of the single-threaded (linear) nature of State. Can we get around this unfortunate linearization?

Speculation

The linearity of carry propagation interferes with parallel execution even when using a tree representation. The problem is that each addB (full adder) invocation must access the carry out from the previous (immediately less significant) bit position and so must wait for that carry to be computed. Since each bit addition must wait for the previous one to finish, we get linear running time, even with unlimited parallel processing available. If we didn’t have to wait for carries, we could instead get logarithmic running time using the tree representation, since subtrees could be added in parallel.

A way out of this dilemma is to speculatively compute the bit sums for both possibilities, i.e., for carry and no carry. We’ll do more work, but much less waiting.

State memoization

Recall the State definition:

newtype State s a = State (s → (a,s))

Rather than using a function of s, let’s use a table indexed by s. Since s is Bool in our use, a table is simply a uniform pair, so we could replace State Bool a with the following:

newtype BoolStateTable a = BST ((a,Bool), (a,Bool))

Exercise: define Functor, Applicative, and Monad instances for BoolStateTable.

Rather than defining such a specialized type, let’s stand back and consider what’s going on. We’re replacing a function by an isomorphic data type. This replacement is exactly what memoization is about. So let’s define a general memoizing state monad:

newtype StateTrie s a = StateTrie (s ⇰ (a,s))

Note that the definition of memoizing state is nearly identical to State. I’ve simply replaced “→” by “⇰”, i.e., memo tries. For the (simple) source code of StateTrie, see the github project. (Poking around on Hackage, I just found monad-memo, which looks related.)

The full-adder function addB is restricted to State, but unnecessarily so. The most general type is inferred as

addB ∷ MonadState Bool m ⇒ (Bool,Bool) → m Bool

where the MonadState class comes from the mtl package.

With the type-generalized addB, we get a more general type for add as well:

add ∷ (Traversable t, Applicative m, MonadState Bool m) ⇒
      t (Bool,Bool) → m (t Bool)
add = traverse addB

Now we can specialize add to work with memoized state:

addLM ∷ [(Bool,Bool)] → StateTrie Bool [Bool]
addLM = add

addTM ∷ Tree (Bool,Bool) → StateTrie Bool (Tree Bool)
addTM = add

What have we done?

The essential tricks in this post are to (a) boost parallelism by speculative evaluation (an old idea) and (b) express speculation as memoization (new, to me at least). The technique wins for binary addition thanks to the small number of possible states, which then makes memoization (full speculation) affordable.

I’m not suggesting that the code above has impressive parallel execution when compiled under GHC. Perhaps it could with some par and pseq annotations. I haven’t tried. This exploration helps me understand a little of the space of hardware-oriented algorithms.

The conditional sum adder looks quite similar to the development above. It has the twist, however, of speculating carries on blocks of a few bits rather than single bits. It’s astonishingly easy to adapt the development above for such a hybrid scheme, forming traversable structures of sequences of bits:

addH ∷ Tree [(Bool,Bool)] → StateTrie Bool (Tree [Bool])
addH = traverse (fromState ∘ add)

I’m using the adapter fromState so that the inner list additions will use State while the outer tree additions will use StateTrie, thanks to type inference. This adapter memoizes and rewraps the state transition function:

fromState ∷ HasTrie s ⇒ State s a → StateTrie s a
fromState = StateTrie ∘ trie ∘ runState

A few recent posts have played with trees from two perspectives. The more commonly used I call "top-down", because the top-level structure is most immediately apparent. A top-down binary tree is either a leaf or a pair of such trees, and that pair can be accessed without wading through intervening structure. Much less commonly used are "bottom-up" trees. A bottom-up binary tree is either a leaf or a single such tree of pairs. In the non-leaf case, the pair structure of the tree elements is accessible by operations like mapping, folding, or scanning. The difference is between a pair of trees and a tree of pairs.

As an alternative to the top-down and bottom-up views on trees, I now want to examine a third view, which is a hybrid of the two. Instead of pairs of trees or trees of pairs, this hybrid view is of trees of trees, and more specifically of bottom-up trees of top-down trees. As we’ll see, these hybrid trees emerge naturally from the top-down and bottom-up views. A later post will show how this third view lends itself to an in-place (destructive) scan algorithm, suitable for execution on modern GPUs.

Edits:

• 2011-06-04: "Suppose we have a bottom-up tree of top-down trees, i.e., t ∷ TB (TT a). Was backwards. (Thanks to Noah Easterly.)
• 2011-06-04: Notation: "f ➶ n" and "f ➴ n".

The post Parallel tree scanning by composition defines "top-down" and a "bottom-up" binary trees as follows (modulo type and constructor names):

data TT a = LT a | BT { unBT ∷ Pair (TT a) } deriving Functor

data TB a = LB a | BB { unBB ∷ TB (Pair a) } deriving Functor

So, while a non-leaf TT (top-down tree) has a pair at the top (outside), a non-leaf TB (bottom-up tree) has pairs at the bottom (inside).

Combining these two observations leads to an interesting possibility. Suppose we have a bottom-up tree of top-down trees, i.e., t ∷ TB (TT a). If t is not a leaf, then t ≡ BB tt where tt is a bottom-up tree whose leaves are pairs of top-down trees, i.e., tt ∷ TB (Pair (TT a)). Each of those leaves of type Pair (TT a) can be converted to type TT a (single tree), simply by applying the BT constructor. Moreover, this transformation is invertible. For convenience, define a type alias for hybrid trees:

type TH a = TB (TT a)

Then the two conversions:

upT   ∷ TH a → TH a
upT   = fmap BT ∘ unBB

downT ∷ TH a → TH a
downT = BB ∘ fmap unBT

  upT ∘ downT
≡ fmap BT ∘ unBB ∘ BB ∘ fmap unBT
≡ fmap BT ∘ fmap unBT
≡ fmap (BT ∘ unBT)
≡ fmap id
≡ id

  downT ∘ upT
≡ BB ∘ fmap unBT ∘ fmap BT ∘ unBB
≡ BB ∘ fmap (unBT ∘ BT) ∘ unBB
≡ BB ∘ fmap id ∘ unBB
≡ BB ∘ id ∘ unBB
≡ BB ∘ unBB
≡ id

Consider a perfect binary leaf tree of depth $n$, i.e., an $n$-deep binary tree with each level full and data only at the leaves (where a leaf is depth $0$ tree.) We can view such a tree as top-down, or bottom-up, or as a hybrid.

Each of these three views is really $n+1$ views:

• Top-down: a depth $n$ tree, or a pair of depth $n-1$ trees, or a pair of pairs of depth $n-2$ trees, etc.
• Bottom-up: a depth $n$ tree, or a depth $n-1$ tree of pairs, or a depth $n-2$ tree of pairs of pairs, etc.
• Hybrid: a depth $n$ tree of depth $0$ trees, or a depth $n-1$ tree of depth $1$ trees, or, …, or a depth $0$ tree of depth $n$ trees.

In the hybrid case, counting from $0$ to $n$, the $k^{th}$ such view is a depth $n-k$ bottom-up tree whose elements (leaf values) are depth $k$ top-down trees. When $k=n$, we have a bottom-up tree whose leaves are all single-leaf trees, and when $k=0$, we have a single-leaf bottom-up tree containing a top-down tree. Imagine a horizontal line at depth $k$, dividing the bottom-up outer structure from the top-down inner structure. The downT function moves the dividing line downward, and the upT function moves the line upward. Both functions are partial.

Generalizing

The role of Pair in the tree types above is simple and regular. We can abstract out this particular type constructor, generalizing to an arbitrary functor. I’ll call this generalization "functor trees". Again, there are top-down and bottom-up versions:

data FT f a = FLT a | FBT { unFBT ∷ f (FT f a) } deriving Functor

data FB f a = FLB a | FBB { unFBB ∷ FB f (f a) } deriving Functor

And a hybrid version, with generalized versions of upT and downT:

type FH f a = FB f (FT f a)

upH   ∷ Functor f ⇒ FH f a → FH f a
upH   = fmap FBT ∘ unFBB

downH ∷ Functor f ⇒ FH f a → FH f a
downH = FBB ∘ fmap unFBT

These definitions specialize to the ones (for binary trees) by substituting Pair for the parameter f.

Depth-typing

The upward and downward view-changing functions above are partial, as they can fail at extreme tree views (at depth $0$ or $n$). We could make this partiality explicit by changing the result type to Maybe (TH a) for binary hybrid trees and to Maybe (FH f a) for the functor generalization. Alternatively, make the tree sizes explicit in the types, as in a few recent posts, including A trie for length-typed vectors. (In those posts, I used the terms "right-folded" and "left-folded" in place of "top-down" and "bottom-up", reflecting the right- or left-folding of functor composition. The "folded" terms led to some confusion, especially in the context of data type folds and scans.) In the depth-typed versions, "leaves" are zero-ary compositions, and "branches" are $(m+1)$-ary compositions for some $m$.

Top-down:

data (➴) ∷ (* → *) → * → (* → *) where
  ZeroT ∷ a → (f ➴ Z) a
  SuccT ∷ IsNat n ⇒ f ((f ➴ n) a) → (f ➴ S n) a

unZeroT ∷ (f ➴ Z) a → a
unZeroT (ZeroT a) = a

unSuccT ∷ (f ➴ S n) a → f ((f ➴ n) a)
unSuccT (SuccT fsa) = fsa

instance Functor f ⇒ Functor (f ➴ n) where
  fmap h (ZeroT a)  = ZeroT (h a)
  fmap h (SuccT fs) = SuccT ((fmap∘fmap) h fs)

Bottom-up:

data (➶) ∷ (* → *) → * → (* → *) where
  ZeroB ∷ a → (f ➶ Z) a
  SuccB ∷ IsNat n ⇒ (f ➶ n) (f a) → (f ➶ S n) a

unZeroB ∷ (f ➶ Z) a → a
unZeroB (ZeroB a) = a

unSuccB ∷ (f ➶ S n) a → (f ➶ n) (f a)
unSuccB (SuccB fsa) = fsa

instance Functor f ⇒ Functor (f ➶ n) where
  fmap h (ZeroB a)  = ZeroB (h a)
  fmap h (SuccB fs) = SuccB ((fmap∘fmap) h fs)

Hybrid:

type H p q f a = (f ➶ p) ((f ➴ q) a)

Upward and downward shift become total functions, and their types explicitly describe how the line shifts between $(p+1)/q$ and $p/(q+1)$:

up   ∷ (Functor f, IsNat q) ⇒ H (S p) q f a → H p (S q) f a
up   = fmap SuccT ∘ unSuccB

down ∷ (Functor f, IsNat p) ⇒ H p (S q) f a → H (S p) q f a
down = SuccB ∘ fmap unSuccT

So what?

Why care about the multitude of views on trees?

• It’s pretty.
• A future post will show how these hybrid trees enable an elegant formulation of parallel scanning that lends itself to an in-place, GPU-friendly implementation.

My last few blog posts have been on the theme of scans, and particularly on parallel scans. In Composable parallel scanning, I tackled parallel scanning in a very general setting. There are five simple building blocks out of which a vast assortment of data structures can be built, namely constant (no value), identity (one value), sum, product, and composition. The post defined parallel prefix and suffix scan for each of these five "functor combinators", in terms of the same scan operation on each of the component functors. Every functor built out of this basic set thus has a parallel scan. Functors defined more conventionally can be given scan implementations simply by converting to a composition of the basic set, scanning, and then back to the original functor. Moreover, I expect this implementation could be generated automatically, similarly to GHC’s DerivingFunctor extension.

Now I’d like to show two examples of parallel scan composition in terms of binary trees, namely the top-down and bottom-up variants of perfect binary leaf trees used in previous posts. (In previous posts, I used the terms "right-folded" and "left-folded" instead of "top-down" and "bottom-up".) The resulting two algorithms are expressed nearly identically, but have differ significantly in the work performed. The top-down version does $\Theta (n\log n)$ work, while the bottom-up version does only $\Theta (n)$, and thus the latter algorithm is work-efficient, while the former is not. Moreover, with a very simple optimization, the bottom-up tree algorithm corresponds closely to Guy Blelloch’s parallel prefix scan for arrays, given in Programming parallel algorithms. I’m delighted with this result, as I had been wondering how to think about Guy’s algorithm.

Edit:

• 2011-05-31: Added Scan and Applicative instances for T2 and T4.

Scanning via functor combinators

In Composable parallel scanning, we saw the Scan class:

class Scan f where
  prefixScan, suffixScan ∷ Monoid m ⇒ f m → (m, f m)

Given a structure of values, the prefix and suffix scan methods generate the overall fold (of type m), plus a structure of the same type as the input. (In contrast, the usual Haskell scanl and scanr functions on lists yield a single list with one more element than the source list. I changed the interface for generality and composability.) The post gave instances for the basic set of five functor combinators.

Most functors are not defined via the basic combinators, but as mentioned above, we can scan by conversion to and from the basic set. For convenience, encapsulate this conversion in a type class:

class EncodeF f where
  type Enc f ∷ * → *
  encode ∷ f a → Enc f a
  decode ∷ Enc f a → f a

and define scan functions via EncodeF:

prefixScanEnc, suffixScanEnc ∷
  (EncodeF f, Scan (Enc f), Monoid m) ⇒ f m → (m, f m)
prefixScanEnc = second decode ∘ prefixScan ∘ encode
suffixScanEnc = second decode ∘ suffixScan ∘ encode

Lists

As a first example, consider

instance EncodeF [] where
  type Enc [] = Const () + Id × []
  encode [] = InL (Const ())
  encode (a : as) = InR (Id a × as)
  decode (InL (Const ())) = []
  decode (InR (Id a × as)) = a : as

And declare a boilerplate Scan instance via EncodeF:

instance Scan [] where
  prefixScan = prefixScanEnc
  suffixScan = suffixScanEnc

I haven’t checked the details, but I think with this instance, suffix scanning has okay performance, while prefix scan does quadratic work. The reason is the in the Scan instance for products, the two components are scanned independently (in parallel), and then the whole second component is adjusted for prefixScan, while the whole first component is adjusted for suffixScan. In the case of lists, the first component is the list head, and second component is the list tail.

For your reading convenience, here’s that Scan instance again:

instance (Scan f, Scan g, Functor f, Functor g) ⇒ Scan (f × g) where
  prefixScan (fa × ga) = (af ⊕ ag, fa' × ((af ⊕) <$> ga'))
   where (af,fa') = prefixScan fa
         (ag,ga') = prefixScan ga

  suffixScan (fa × ga) = (af ⊕ ag, ((⊕ ag) <$> fa') × ga')
   where (af,fa') = suffixScan fa
         (ag,ga') = suffixScan ga

The lop-sidedness of the list type thus interferes with parallelization, and makes the parallel scans perform much worse than cumulative sequential scans.

Let’s next look at a more balanced type.

Binary Trees

We’ll get better parallel performance by organizing our data so that we can cheaply partition it into roughly equal pieces. Tree types allows such partitioning.

Top-down trees

We’ll try a few variations, starting with a simple binary tree.

data T1 a = L1 a | B1 (T1 a) (T1 a) deriving Functor

Encoding and decoding is straightforward:

instance EncodeF T1 where
  type Enc T1 = Id + T1 × T1
  encode (L1 a)   = InL (Id a)
  encode (B1 s t) = InR (s × t)
  decode (InL (Id a))  = L1 a
  decode (InR (s × t)) = B1 s t

instance Scan T1 where
  prefixScan = prefixScanEnc
  suffixScan = suffixScanEnc

Note that these definitions could be generated automatically from the data type definition.

For balanced trees, prefix and suffix scan divide the problem in half at each step, solve each half, and do linear work to patch up one of the two halves. Letting $n$ be the number of elements, and $W(n)$ the work, we have the recurrence $W(n)=2W(n/2)+cn$ for some constant factor $c$. By the Master theorem, therefore, the work done is $\Theta (n\log n)$. (Use case 2, with $a=b=2$, $f(n)=cn$, and $k=0$.)

Again assuming a balanced tree, the computation dependencies have logarithmic depth, so the ideal parallel running time (assuming sufficient processors) is $\Theta (\log n)$. Thus we have an algorithm that is depth-efficient (modulo constant factors) but work-inefficient.

Composition

A binary tree as defined above is either a leaf or a pair of binary trees. We can make this pair-ness more explicit with a reformulation:

data T2 a = L2 a | B2 (Pair (T2 a)) deriving Functor

where Pair, as in Composable parallel scanning, is defined as

data Pair a = a :# a deriving Functor

or even

type Pair = Id × Id

For encoding and decoding, we could use the same representation as with T1, but let’s instead use a more natural one for the definition of T2:

instance EncodeF T2 where
  type Enc T2 = Id + Pair ∘ T2
  encode (L2 a)  = InL (Id a)
  encode (B2 st) = InR (O st)
  decode (InL (Id a)) = L2 a
  decode (InR (O st)) = B2 st

Boilerplate scanning:

instance Scan T2 where
  prefixScan = prefixScanEnc
  suffixScan = suffixScanEnc

for which we’ll need an applicative instance:

instance Applicative T2 where
  pure = L2
  L2 f <*> L2 x = L2 (f x)
  B2 (fs :# gs) <*> B2 (xs :# ys) = B2 ((fs <*> xs) :# (gs <*> ys))
  _ <*> _ = error "T2 (<*>): structure mismatch"

The O constructor is for functor composition.

With a small change to the tree type, we can make the composition of Pair and T more explicit:

data T3 a = L3 a | B3 ((Pair ∘ T3) a) deriving Functor

Then the conversion becomes even simpler, since there’s no need to add or remove O wrappers:

instance EncodeF T3 where
  type Enc T3 = Id + Pair ∘ T3
  encode (L3 a)  = InL (Id a)
  encode (B3 st) = InR st
  decode (InL (Id a)) = L3 a
  decode (InR st)     = B3 st

Bottom-up trees

In the formulations above, a non-leaf tree consists of a pair of trees. I’ll call these trees "top-down", since visible pair structure begins at the top.

With a very small change, we can instead use a tree of pairs:

data T4 a = L4 a | B4 (T4 (Pair a)) deriving Functor

Again an applicative instance allows a standard Scan instance:

instance Scan T4 where
  prefixScan = prefixScanEnc
  suffixScan = suffixScanEnc

instance Applicative T4 where
  pure = L4
  L4 f   <*> L4 x   = L4 (f x)
  B4 fgs <*> B4 xys = B4 (liftA2 h fgs xys)
   where h (f :# g) (x :# y) = f x :# g y
  _ <*> _ = error "T4 (<*>): structure mismatch"

or a more explicitly composed form:

data T5 a = L5 a | B5 ((T5 ∘ Pair) a) deriving Functor

I’ll call these new variations "bottom-up" trees, since visible pair structure begins at the bottom. After stripping off the branch constructor, B4, we can get at the pair-valued leaves by means of fmap, fold, or traverse (or variations). For B5, we’d also have to strip off the O wrapper (functor composition).

Encoding is nearly the same as with top-down trees. For instance,

instance EncodeF T4 where
  type Enc T4 = Id + T4 ∘ Pair
  encode (L4 a) = InL (Id a)
  encode (B4 t) = InR (O t)
  decode (InL (Id a)) = L4 a
  decode (InR (O t))  = B4 t

Scanning pairs

We’ll need to scan on the Pair functor. If we use the definition of Pair above in terms of Id and (×), then we’ll get scanning for free. For using Pair, I find the explicit data type definition above more convenient. We can then derive a Scan instance by conversion. Start with a standard specification:

data Pair a = a :# a deriving Functor

And encode & decode explicitly:

instance EncodeF Pair where
  type Enc Pair = Id × Id
  encode (a :# b) = Id a × Id b
  decode (Id a × Id b) = a :# b

Then use our boilerplate Scan instance for EncodeF instances:

instance Scan Pair where
  prefixScan = prefixScanEnc
  suffixScan = suffixScanEnc

We’ve seen the Scan instance for (×) above. The instance for Id is very simple:

newtype Id a = Id a

instance Scan Id where
  prefixScan (Id m) = (m, Id ∅)
  suffixScan        = prefixScan

Given these definitions, we can calculate a more streamlined Scan instance for Pair:

  prefixScan (a :# b)
≡  {- specification -}
  prefixScanEnc (a :# b)
≡  {- prefixScanEnc definition -}
  (second decode ∘ prefixScan ∘ encode) (a :# b)
≡  {- (∘) -}
  second decode (prefixScan (encode (a :# b)))
≡  {- encode definition for Pair -}
  second decode (prefixScan (Id a × Id b))
≡  {- prefixScan definition for f × g -}
  second decode
    (af ⊕ ag, fa' × ((af ⊕) <$> ga'))
     where (af,fa') = prefixScan (Id a)
           (ag,ga') = prefixScan (Id b)
≡  {- Definition of second on functions -}
  (af ⊕ ag, decode (fa' × ((af ⊕) <$> ga')))
   where (af,fa') = prefixScan (Id a)
         (ag,ga') = prefixScan (Id b)
≡  {- prefixScan definition for Id -}
  (af ⊕ ag, decode (fa' × ((af ⊕) <$> ga')))
   where (af,fa') = (a, Id ∅)
         (ag,ga') = (b, Id ∅)
≡  {- substitution -}
  (a ⊕ b, decode (Id ∅ × ((a ⊕) <$> Id ∅)))
≡  {- fmap/(<$>) for Id -}
  (a ⊕ b, decode (Id ∅ × Id (a ⊕ ∅)))
≡  {- Monoid law -}
  (a ⊕ b, decode (Id ∅ × Id a))
≡  {- decode definition on Pair -}
  (a ⊕ b, (∅ :# a))

Whew! And similarly for suffixScan.

Now let’s recall the Scan instance for Pair given in Composable parallel scanning:

instance Scan Pair where
  prefixScan (a :# b) = (a ⊕ b, (∅ :# a))
  suffixScan (a :# b) = (a ⊕ b, (b :# ∅))

Hurray! The derivation led us to the same definition. A "sufficiently smart" compiler could do this derivation automatically.

With this warm-up derivation, let’s now turn to trees.

Scanning trees

Given the tree encodings above, how does scan work? We’ll have to consult Scan instances for some of the functor combinators. The product instance is repeated above. We’ll also want the instances for sum and composition. Omitting the suffixScan definitions for brevity:

data (f + g) a = InL (f a) | InR (g a)

instance (Scan f, Scan g) ⇒ Scan (f + g) where
  prefixScan (InL fa) = second InL (prefixScan fa)
  prefixScan (InR ga) = second InR (prefixScan ga)

newtype (g ∘ f) a = O (g (f a))

instance (Scan g, Scan f, Functor f, Applicative g) ⇒ Scan (g ∘ f) where
  prefixScan = second (O ∘ fmap adjustL ∘ zip)
             ∘ assocR
             ∘ first prefixScan
             ∘ unzip
             ∘ fmap prefixScan
             ∘ unO

This last definition uses a few utility functions:

zip ∷ Applicative g ⇒ (g a, g b) → g (a,b)
zip = uncurry (liftA2 (,))

unzip ∷ Functor g ⇒ g (a,b) → (g a, g b)
unzip = fmap fst &&& fmap snd

assocR ∷ ((a,b),c) → (a,(b,c))
assocR   ((a,b),c) =  (a,(b,c))

adjustL ∷ (Functor f, Monoid m) ⇒ (m, f m) → f m
adjustL (m, ms) = (m ⊕) <$> ms

Let’s consider how the Scan (g ∘ f) instance plays out for top-down vs bottom-up trees, given the functor-composition encodings above. The critical definitions:

type Enc T2 = Id + Pair ∘ T2

type Enc T4 = Id + T4 ∘ Pair

Focusing on the branch case, we have Pair ∘ T2 vs T4 ∘ Pair, so we’ll use the Scan (g ∘ f) instance either way. Let’s consider the work implied by that instance. There are two calls to prefixScan, plus a linear amount of other work. The meanings of those two calls differ, however:

• For top-down trees (T2), the recursive tree scans are in fmap prefixScan, mapping over the pair of trees. The first prefixScan is a pair scan and so does constant work. Since there are two recursive calls, each working on a tree of half size (assuming balance), plus linear other work, the total work $\Theta (n\log n)$, as explained above.
• For bottom-up trees (T4), there is only one recursive recursive tree scan, which appears in first prefixScan. The prefixScan in fmap prefixScan is pair scan and so does constant work but is mapped over the half-sized tree (of pairs), and so does linear work altogether. Since there only one recursive tree scan, at half size, plus linear other work, the total work is then proportional to $n+n/2+n/4+\dots \approx 2n=\Theta (n)$. So we have a work-efficient algorithm!

Looking deeper

In addition to the simple analysis above of scanning over top-down and over bottom-up, let’s look in detail at what transpires and how each case can be optimized. This section may well have more detail than you’re interested in. If so, feel free to skip ahead.

Top-down

Beginning as with Pair,

  prefixScan t
≡  {- specification -}
  prefixScanEnc t
≡  {- prefixScanEnc definition -}
  (second decode ∘ prefixScan ∘ encode) t
≡  {- (∘) -}
  second decode (prefixScan (encode t))

Take T2, with T3 being quite similar. Now split into two cases for the two constructors of T2. First leaf:

  prefixScan (L2 m)
≡  {- as above -}
  second decode (prefixScan (encode (L2 m)))
≡  {- encode for L2 -}
  second decode (prefixScan (InL (Id m)))
≡  {- prefixScan for functor sum -}
  second decode (second InL (prefixScan (Id m)))
≡  {- prefixScan for Id -}
  second decode (second InL (m, Id ∅))
≡  {- second for functions -}
  second decode (m, InL (Id ∅))
≡  {- second for functions -}
  (m, decode (InL (Id ∅)))
≡  {- decode for L2 -}
  (m, L2 ∅)

Then branch:

  prefixScan (B2 (s :# t))
≡  {- as above -}
  second decode (prefixScan (encode (B2 (s :# t))))
≡  {- encode for L2 -}
  second decode (prefixScan (InR (O (s :# t))))
≡  {- prefixScan for (+) -}
  second decode (second InR (prefixScan (O (s :# t))))
≡  {- property of second -}
  second (decode ∘ InR) (prefixScan (O (s :# t)))

Focus on the prefixScan application:

  prefixScan (O (s :# t)) =
≡  {- prefixScan for (∘) -}
 ( second (O ∘ fmap adjustL ∘ zip) ∘ assocR ∘ first prefixScan
 ∘ unzip ∘ fmap prefixScan ∘ unO ) (O (s :# t))
≡  {- unO/O -}
  ( second (O ∘ fmap adjustL ∘ zip) ∘ assocR ∘ first prefixScan
  ∘ unzip ∘ fmap prefixScan ) (s :# t)
≡  {- fmap on Pair -}
  (second (O ∘ fmap adjustL ∘ zip) ∘ assocR ∘ first prefixScan ∘ unzip)
    (prefixScan s :# prefixScan t)
≡  {- expand prefixScan -}
  (second (O ∘ fmap adjustL ∘ zip) ∘ assocR ∘ first prefixScan ∘ unzip)
    ((ms,s') :# (mt,t'))
      where (ms,s') = prefixScan s
            (mt,t') = prefixScan t
≡  {- unzip -}
  (second (O ∘ fmap adjustL ∘ zip) ∘ assocR ∘ first prefixScan)
    ((ms :# mt), (s' :# t')) where ⋯
≡  {- first -}
  (second (O ∘ fmap adjustL ∘ zip) ∘ assocR)
    (prefixScan (ms :# mt), (s' :# t')) where ⋯
≡  {- prefixScan for Pair -}
  (second (O ∘ fmap adjustL ∘ zip) ∘ assocR)
    ((ms ⊕ mt, (∅ :# ms)), (s' :# t')) where ⋯
≡  {- assocR -}
  (second (O ∘ fmap adjustL ∘ zip))
    (ms ⊕ mt, ((∅ :# ms), (s' :# t'))) where ⋯
≡  {- second -}
  ( ms ⊕ mt
  , (O ∘ fmap adjustL ∘ zip) ((∅ :# ms), (s' :# t')) ) where ⋯
≡  {- zip -}
  ( ms ⊕ mt
  , (O ∘ fmap adjustL) ((∅,s') :# (ms,t')) )  where ⋯
≡  {- fmap for Pair -}
  ( ms ⊕ mt
  , O (adjustL (∅,s') :# adjustL (ms,t')) )  where ⋯
≡  {- adjustL -}
  ( ms ⊕ mt
  , O (((∅ ⊕) <$> s') :# ((ms ⊕) <$> t')) )  where ⋯
≡  {- Monoid law (left identity) -}
  ( ms ⊕ mt
  , O ((id <$> s') :# ((ms ⊕) <$> t')) )  where ⋯
≡  {- Functor law (fmap id) -}
  ( ms ⊕ mt
  , O (s' :# ((ms ⊕) <$> t')) )
      where (ms,s') = prefixScan s
            (mt,t') = prefixScan t

Continuing from above,

  prefixScan (B2 (s :# t))
≡  {- see above -}
  second (decode ∘ InR) (prefixScan (O (s :# t)))
≡  {- prefixScan focus from above -}
  second (decode ∘ InR)
    ( ms ⊕ mt
    , O (s' :# ((ms ⊕) <$> t')) )
        where (ms,s') = prefixScan s
              (mt,t') = prefixScan t
≡  {- definition of second on functions -}
    (ms ⊕ mt, (decode ∘ InR) (O (s' :# ((ms ⊕) <$> t')))) where ⋯
≡  {- (∘) -}
    (ms ⊕ mt, decode (InR (O (s' :# ((ms ⊕) <$> t'))))) where ⋯
≡  {- decode for B2 -}
    (ms ⊕ mt, B2 (s' :# ((ms ⊕) <$> t'))) where ⋯

This final form is as in Deriving parallel tree scans, changed for the new scan interface. The derivation saved some work in wrapping & unwrapping and method invocation, plus one of the two adjustment passes over the sub-trees. As explained above, this algorithm performs $\Theta (n\log n)$ work.

I’ll leave suffixScan for you to do yourself.

Bottom-up

What happens if we switch from top-down to bottom-up binary trees? I’ll use T4 (though T5 would work as well):

data T4 a = L4 a | B4 (T4 (Pair a))

The leaf case is just as with T2 above, so let’s get right to branches.

  prefixScan (B4 t)
≡  {- as above -}
  second decode (prefixScan (encode (B4 t)))
≡  {- encode for L2 -}
  second decode (prefixScan (InR (O t)))
≡  {- prefixScan for (+) -}
  second decode (second InR (prefixScan (O t)))
≡  {- property of second -}
  second (decode ∘ InR) (prefixScan (O t))

As before, now focus on the prefixScan call.

  prefixScan (O t) =
≡  {- prefixScan for (∘) -}
 ( second (O ∘ fmap adjustL ∘ zip) ∘ assocR ∘ first prefixScan
 ∘ unzip ∘ fmap prefixScan ∘ unO ) (O t)
≡  {- unO/O -}
  ( second (O ∘ fmap adjustL ∘ zip) ∘ assocR ∘ first prefixScan
  ∘ unzip ∘ fmap prefixScan ) t
≡  {- prefixScan on Pair (derived above) -}
  (second (O ∘ fmap adjustL ∘ zip) ∘ assocR ∘ first prefixScan ∘ unzip)
    fmap (λ (a :# b) → (a ⊕ b, (∅ :# a))) t
≡  {- unzip/fmap -}
  (second (O ∘ fmap adjustL ∘ zip) ∘ assocR ∘ first prefixScan)
    ( fmap (λ (a :# b) → (a ⊕ b)) t
    , fmap (λ (a :# b) → (∅ :# a))   t )
≡  {- first on functions -}
  (second (O ∘ fmap adjustL ∘ zip) ∘ assocR)
    ( prefixScan (fmap (λ (a :# b) → (a ⊕ b)) t)
    , fmap (λ (a :# b) → (∅ :# a))   t )
≡  {- expand prefixScan -}
  (second (O ∘ fmap adjustL ∘ zip) ∘ assocR)
    ((mp,p'), fmap (λ (a :# b) → (∅ :# a)) t)
   where (mp,p') = prefixScan (fmap (λ (a :# b) → (a ⊕ b)) t)
≡  {- assocR -}
  (second (O ∘ fmap adjustL ∘ zip))
    (mp, (p', fmap (λ (a :# b) → (∅ :# a)) t))
   where ⋯
≡  {- second on functions -}
  (mp, (O ∘ fmap adjustL ∘ zip) (p', fmap (λ (a :# b) → (∅ :# a)) t))
    where ⋯
≡  {- fmap/zip/fmap -}
  (mp, O (liftA2 tweak p' t))
    where tweak s (a :# _) = adjustL (s, (∅ :# a))
          (mp,p') = prefixScan (fmap (λ (a :# b) → (a ⊕ b)) t)
≡  {- adjustL, then simplify -}
  (mp, O (liftA2 tweak p' t))
    where tweak s (a :# _) = (s :# s ⊕ a)
          (mp,p') = prefixScan (fmap (λ (a :# b) → (a ⊕ b)) t)

Now re-introduce the context of prefixScan (O t):

  prefixScan (B4 t)
≡  {- see above -}
  second (decode ∘ InR) (prefixScan (O t))
≡  {- see above -}
  second (decode ∘ InR)
    (mp, O (liftA2 tweak p' t))
      where ⋯
≡  {- decode for T4 -}
  (mp, B4 (liftA2 tweak p' t))
    where p = fmap (λ (e :# o) → (e ⊕ o)) t
          (mp,p') = prefixScan p
          tweak s (e :# _) = (s :# s ⊕ e)

Notice how much this bottom-up tree scan algorithm differs from the top-down algorithm derived above. In particular, there’s only one recursive tree scan (on a half-sized tree) instead of two, plus linear additional work, for a total of $\Theta (n)$ work.

Guy Blelloch’s parallel scan algorithm

In Programming parallel algorithms, Guy Blelloch gives the following algorithm for parallel prefix scan, expressed in the parallel functional language NESL:

function scan(a) =
if #a ≡ 1 then [0]
else
  let es = even_elts(a);
      os = odd_elts(a);
      ss = scan({e+o: e in es; o in os})
  in interleave(ss,{s+e: s in ss; e in es})

This algorithm is nearly identical to the T4 scan algorithm above. I was very glad to find this route to Guy’s algorithm, which had been fairly mysterious to me. I mean, I could believe that the algorithm worked, but I had no idea how I might have discovered it myself. With the functor composition approach to scanning, I now see how Guy’s algorithm emerges as well as how it generalizes to other data structures.

Nested data types and parallelism

Most of the recursive algebraic data types that appear in Haskell programs are regular, meaning that the recursive instances are instantiated with the same type parameter as the containing type. For instance, a top-down tree of elements of type a is either a leaf or has two trees whose elements have that same type a. In contrast, in a bottom-up tree, the (single) recursively contained tree is over elements of type (a,a). Such non-regular data types are called "nested". The two tree scan algorithms above suggest to me that nested data types are particularly useful for efficient parallel algorithms.

The post Deriving list scans gave a simple specification of the list-scanning functions scanl and scanr, and then transformed those specifications into the standard optimized implementations. Next, the post Deriving parallel tree scans adapted the specifications and derivations to a type of binary trees. The resulting implementations are parallel-friendly, but not work-efficient, in that they perform $n\log n$ work vs linear work as in the best-known sequential algorithm.

Besides the work-inefficiency, I don’t know how to extend the critical initTs and tailTs functions (analogs of inits and tails on lists) to depth-typed, perfectly balanced trees, of the sort I played with in A trie for length-typed vectors and From tries to trees. The difficulty I encounter is that the functions initTs and tailTs make unbalanced trees out of balanced ones, so I don’t know how to adapt the specifications when types prevent the existence of unbalanced trees.

This new post explores an approach to generalized scanning via type classes. After defining the classes and giving a simple example, I’ll give a simple & general framework based on composing functor combinators.

Edits:

• 2011-03-02: Fixed typo. "constant functor is easiest" (instead of "identity functor"). Thanks, frguybob.
• 2011-03-05: Removed final unfinished sentence.
• 2011-07-28: Replace "assocL" with "assocR" in prefixScan derivation for g ∘ f.

Generalizing list scans

The left and right scan functions on lists have an awkward feature. The output list has one more element than the input list, corresponding to the fact that the number of prefixes (inits) of a list is one more than the number of elements, and similarly for suffixes (tails).

While it’s easy to extend a list by adding one more element, it’s not easy with other functors. In Deriving parallel tree scans, I simply removed the ∅ element from the scan. In this post, I’ll instead change the interface to produce an output of exactly the same shape, plus one extra element. The extra element will equal a fold over the complete input. If you recall, we had to search for that complete fold in an input subtree in order to adjust the other subtree. (See headT and lastT and their generalizations in Deriving parallel tree scans.) Separating out this value eliminates the search.

Define a class with methods for prefix and suffix scan:

class Scan f where
  prefixScan, suffixScan ∷ Monoid m ⇒ f m → (m, f m)

Prefix scans (prefixScan) accumulate moving left-to-right, while suffix scans (suffixScan) accumulate moving right-to-left.

A simple example: pairs

To get a first sense of generalized scans, let’s use see how to scan over a pair functor.

data Pair a = a :# a deriving (Eq,Ord,Show)

With GHC’s DeriveFunctor option, we could also derive a Functor instance, but for clarity, define it explicitly:

instance Functor Pair where
  fmap f (a :# b) = (f a :# f b)

The scans:

instance Scan Pair where
  prefixScan (a :# b) = (a ⊕ b, (∅ :# a))
  suffixScan (a :# b) = (a ⊕ b, (b :# ∅))

As you can see, if we eliminated the ∅ elements, we could shift to the left or right and forgo the extra result.

Naturally, there is also a Fold instance, and the scans produce the fold results as well sub-folds:

instance Foldable Pair where
  fold (a :# b) = a ⊕ b

The Pair functor also has unsurprising instances for Applicative and Traversable.

instance Applicative Pair where
  pure a = a :# a
  (f :# g) <*> (x :# y) = (f x :# g y)

instance Traversable Pair where
  sequenceA (fa :# fb) = (:#) <$> fa <*> fb

We don’t really have to figure out how to define scans for every functor separately. We can instead look at how functors are are composed out of their essential building blocks.

Scans for functor combinators

To see how to scan over a broad range of functors, let’s look at each of the functor combinators, e.g., as in Elegant memoization with higher-order types.

Constant

The constant functor is easiest.

newtype Const x a = Const x

There are no values to accumulate, so the final result (fold) is ∅.

instance Scan (Const x) where
  prefixScan (Const x) = (∅, Const x)
  suffixScan           = prefixScan

Identity

The identity functor is nearly as easy.

newtype Id a = Id a

instance Scan Id where
  prefixScan (Id m) = (m, Id ∅)
  suffixScan        = prefixScan

Sum

Scanning in a sum is just scanning in a summand:

data (f + g) a = InL (f a) | InR (g a)

instance (Scan f, Scan g) ⇒ Scan (f + g) where
  prefixScan (InL fa) = second InL (prefixScan fa)
  prefixScan (InR ga) = second InR (prefixScan ga)

  suffixScan (InL fa) = second InL (suffixScan fa)
  suffixScan (InR ga) = second InR (suffixScan ga)

These definitions correspond to simple "commutative diagram" properties, e.g.,

prefixScan ∘ InL ≡ second InL ∘ prefixScan

Product

Product scannning is a little trickier.

data (f × g) a = f a × g a

Scan each of the two parts separately, and then combine the final (fold) part of one result with each of the non-final elements of the other.

instance (Scan f, Scan g, Functor f, Functor g) ⇒ Scan (f × g) where
  prefixScan (fa × ga) = (af ⊕ ag, fa' × ((af ⊕) <$> ga'))
   where (af,fa') = prefixScan fa
         (ag,ga') = prefixScan ga

  suffixScan (fa × ga) = (af ⊕ ag, ((⊕ ag) <$> fa') × ga')
   where (af,fa') = suffixScan fa
         (ag,ga') = suffixScan ga

Composition

Finally, composition is the trickiest.

newtype (g ∘ f) a = O (g (f a))

The target signatures:

  prefixScan, suffixScan ∷ Monoid m ⇒ (g ∘ f) m → (m, (g ∘ f) m)

To find the prefix and suffix scan definitions, fiddle with types beginning at the domain type for prefixScan or suffixScan and arriving at the range type.

Some helpers:

zip ∷ Applicative g ⇒ (g a, g b) → g (a,b)
zip = uncurry (liftA2 (,))

unzip ∷ Functor g ⇒ g (a,b) → (g a, g b)
unzip = fmap fst &&& fmap snd

assocR ∷ ((a,b),c) → (a,(b,c))
assocR   ((a,b),c) =  (a,(b,c))

adjustL ∷ (Functor f, Monoid m) ⇒ (m, f m) → f m
adjustL (m, ms) = (m ⊕) <$> ms

adjustR ∷ (Functor f, Monoid m) ⇒ (m, f m) → f m
adjustR (m, ms) = (⊕ m) <$> ms

First prefixScan:

gofm                     ∷ (g ∘ f) m
unO                   '' ∷ g (f m)
fmap prefixScan       '' ∷ g (m, f m)
unzip                 '' ∷ (g m, g (f m))
first prefixScan      '' ∷ ((m, g m), g (f m))
assocR                '' ∷ (m, (g m, g (f m)))
second zip            '' ∷ (m, g (m, f m))
second (fmap adjustL) '' ∷ (m, g (f m))
second O              '' ∷ (m, (g ∘ f) m)

Then suffixScan:

gofm                     ∷ (g ∘ f) m
unO                   '' ∷ g (f m)
fmap suffixScan       '' ∷ g (m, f m)
unzip                 '' ∷ (g m, g (f m))
first suffixScan      '' ∷ ((m, g m), g (f m))
assocR                '' ∷ (m, (g m, g (f m)))
second zip            '' ∷ (m, (g (m, f m)))
second (fmap adjustR) '' ∷ (m, (g (f m)))
second O              '' ∷ (m, ((g ∘ f) m))

Putting together the pieces and simplifying just a bit leads to the method definitions:

instance (Scan g, Scan f, Functor f, Applicative g) ⇒ Scan (g ∘ f) where
  prefixScan = second (O ∘ fmap adjustL ∘ zip)
             ∘ assocR
             ∘ first prefixScan
             ∘ unzip
             ∘ fmap prefixScan
             ∘ unO

  suffixScan = second (O ∘ fmap adjustR ∘ zip)
             ∘ assocR
             ∘ first suffixScan
             ∘ unzip
             ∘ fmap suffixScan
             ∘ unO

What’s coming up?

• What might not easy to spot at this point is that the prefixScan and suffixScan methods given in this post do essentially the same job as in Deriving parallel tree scans, when the binary tree type is deconstructed into functor combinators. A future post will show this connection.
• Switch from standard (right-folded) trees to left-folded trees (in the sense of A trie for length-typed vectors and From tries to trees), which reduces the running time from $\Theta (n\log n)$ to $\Theta n$.
• Scanning in place, i.e., destructively replacing the values in the input structure rather than allocating a new structure.

The post Deriving list scans explored folds and scans on lists and showed how the usual, efficient scan implementations can be derived from simpler specifications.

Let’s see now how to apply the same techniques to scans over trees.

This new post is one of a series leading toward algorithms optimized for execution on massively parallel, consumer hardware, using CUDA or OpenCL.

Edits:

• 2011-03-01: Added clarification about "∅" and "(⊕)".
• 2011-03-23: corrected "linear-time" to "linear-work" in two places.

Trees

Our trees will be non-empty and binary:

data T a = Leaf a | Branch (T a) (T a)

instance Show a ⇒ Show (T a) where
  show (Leaf a)     = show a
  show (Branch s t) = "("++show s++","++show t++")"

Nothing surprising in the instances:

instance Functor T where
  fmap f (Leaf a)     = Leaf (f a)
  fmap f (Branch s t) = Branch (fmap f s) (fmap f t)

instance Foldable T where
  fold (Leaf a)     = a
  fold (Branch s t) = fold s ⊕ fold t

instance Traversable T where
  sequenceA (Leaf a)     = fmap Leaf a
  sequenceA (Branch s t) =
    liftA2 Branch (sequenceA s) (sequenceA t)

BTW, my type-setting software uses "∅" and "(⊕)" for Haskell’s "mempty" and "mappend".

Also handy will be extracting the first and last (i.e., leftmost and rightmost) leaves in a tree:

headT ∷ T a → a
headT (Leaf a)       = a
headT (s `Branch` _) = headT s

lastT ∷ T a → a
lastT (Leaf a)       = a
lastT (_ `Branch` t) = lastT t

headT ∘ fmap f ≡ f ∘ headT
lastT ∘ fmap f ≡ f ∘ lastT

  headT (fmap f (Leaf a))
≡  {- fmap on T -}
  headT (Leaf (f a))
≡  {- headT def -}
  f a
≡  {- headT def -}
  f (headT (Leaf a))

  headT (fmap f (Branch s t))
≡  {- fmap on T -}
  headT (Branch (fmap f s) (fmap f t))
≡  {- headT def -}
  headT (fmap f s)
≡  {- induction -}
  f (headT s)
≡  {- headT def -}
  f (headT (Branch s t))

From lists to trees and back

We can flatten trees into lists:

flatten ∷ T a → [a]
flatten = fold ∘ fmap (:[])

Equivalently, using foldMap:

flatten = foldMap (:[])

Alternatively, we could define fold via flatten:

instance Foldable T where fold = fold ∘ flatten

flatten ∷ T a → [a]
flatten (Leaf a)     = [a]
flatten (Branch s t) = flatten s ++ flatten t

We can also "unflatten" lists into balanced trees:

unflatten ∷ [a] → T a
unflatten []  = error "unflatten: Oops! Empty list"
unflatten [a] = Leaf a
unflatten xs  = Branch (unflatten prefix) (unflatten suffix)
 where
   (prefix,suffix) = splitAt (length xs `div` 2) xs

Both flatten and unflatten can be implemented more efficiently.

For instance,

t1,t2 ∷ T Int
t1 = unflatten [1‥3]
t2 = unflatten [1‥16]

*T> t1
(1,(2,3))
*T> t2
((((1,2),(3,4)),((5,6),(7,8))),(((9,10),(11,12)),((13,14),(15,16))))

Specifying tree scans

Prefixes and suffixes

The post Deriving list scans gave specifications for list scanning in terms of inits and tails. One consequence of this specification is that the output of scanning has one more element than the input. Alternatively, we could use non-empty variants of inits and tails, so that the input & output are in one-to-one correspondence.

inits' ∷ [a] → [[a]]
inits' []     = []
inits' (x:xs) = map (x:) ([] : inits' xs)

The cons case can also be written as

inits' (x:xs) = [x] : map (x:) (inits' xs)

tails' ∷ [a] → [[a]]
tails' []         = []
tails' xs@(_:xs') = xs : tails' xs'

For instance,

*T> inits' "abcd"
["a","ab","abc","abcd"]
*T> tails' "abcd"
["abcd","bcd","cd","d"]

Our tree functor has a symmetric definition, so we get more symmetry in the counterparts to inits' and tails':

initTs ∷ T a → T (T a)
initTs (Leaf a)       = Leaf (Leaf a)
initTs (s `Branch` t) =
  Branch (initTs s) (fmap (s `Branch`) (initTs t))

tailTs ∷ T a → T (T a)
tailTs (Leaf a)       = Leaf (Leaf a)
tailTs (s `Branch` t) =
  Branch (fmap (`Branch` t) (tailTs s)) (tailTs t)

Try it:

*T> t1
(1,(2,3))
*T> initTs t1
(1,((1,2),(1,(2,3))))
*T> tailTs t1
((1,(2,3)),((2,3),3))

*T> unflatten [1‥5]
((1,2),(3,(4,5)))
*T> initTs (unflatten [1‥5])
((1,(1,2)),(((1,2),3),(((1,2),(3,4)),((1,2),(3,(4,5))))))
*T> tailTs (unflatten [1‥5])
((((1,2),(3,(4,5))),(2,(3,(4,5)))),((3,(4,5)),((4,5),5)))

lastT ∘ initTs ≡ id
headT ∘ tailTs ≡ id

  lastT (initTs (Leaf a))
≡  {- initTs def -}
  lastT (Leaf (Leaf a))
≡  {- lastT def -}
  Leaf a

  lastT (initTs (s `Branch` t))
≡  {- initTs def -}
  lastT (Branch (⋯) (fmap (s `Branch`) (initTs t)))
≡  {- lastT def -}
  lastT (fmap (s `Branch`) (initTs t))
≡  {- lastT ∘ fmap f -}
  (s `Branch`) (lastT (initTs t))
≡  {- trivial -}
  s `Branch` lastT (initTs t)
≡  {- induction -}
  s `Branch` t

Scan specification

Now we can specify prefix & suffix scanning:

scanlT, scanrT ∷ Monoid a ⇒ T a → T a
scanlT = fmap fold ∘ initTs
scanrT = fmap fold ∘ tailTs

Try it out:

t3 ∷ T String
t3 = fmap (:[]) (unflatten "abcde")

*T> t3
(("a","b"),("c",("d","e")))
*T> scanlT t3
(("a","ab"),("abc",("abcd","abcde")))
*T> scanrT t3
(("abcde","bcde"),("cde",("de","e")))

To test on numbers, I’ll use a handy notation from Matt Hellige to add pre- and post-processing:

(↝) ∷ (a' → a) → (b → b') → ((a → b) → (a' → b'))
(f ↝ h) g = h ∘ g ∘ f

And a version specialized to functors:

(↝*) ∷ Functor f ⇒ (a' → a) → (b → b')
     → (f a → f b) → (f a' → f b')
f ↝* g = fmap f ↝ fmap g

t4 ∷ T Integer
t4 = unflatten [1‥6]

t5 ∷ T Integer
t5 = (Sum ↝* getSum) scanlT t4

Try it:

*T> t4
((1,(2,3)),(4,(5,6)))
*T> initTs t4
((1,((1,2),(1,(2,3)))),(((1,(2,3)),4),(((1,(2,3)),(4,5)),((1,(2,3)),(4,(5,6))))))
*T> t5
((1,(3,6)),(10,(15,21)))

fold ≡ lastT ∘ scanlT
fold ≡ headT ∘ scanrT

  lastT ∘ scanlT
≡  {- scanlT spec -}
  lastT ∘ fmap fold ∘ initTs
≡  {- lastT ∘ fmap f -}
  fold ∘ lastT ∘ initTs
≡  {- lastT ∘ initTs -}
  fold

  headT ∘ scanrT 
≡  {- scanrT def -}
  headT ∘ fmap fold ∘ tailTs
≡  {- headT ∘ fmap f -}
  fold ∘ headT ∘ tailTs
≡  {- headT ∘ tailTs -}
  fold

*T> fold t3
"abcde"
*T> (lastT ∘ scanlT) t3
"abcde"
*T> (headT ∘ scanrT) t3
"abcde"

Deriving faster scans

Recall the specifications:

scanlT = fmap fold ∘ initTs
scanrT = fmap fold ∘ tailTs

To derive more efficient implementations, proceed as in Deriving list scans. Start with prefix scan (scanlT), and consider the Leaf and Branch cases separately.

  scanlT (Leaf a)
≡  {- scanlT spec -}
  fmap fold (initTs (Leaf a))
≡  {- initTs def -}
  fmap fold (Leaf (Leaf a))
≡  {- fmap def -}
  Leaf (fold (Leaf a))
≡  {- fold def -}
  Leaf a

  scanlT (s `Branch` t)
≡  {- scanlT spec -}
  fmap fold (initTs (s `Branch` t))
≡  {- initTs def -}
  fmap fold (Branch (initTs s) (fmap (s `Branch`) (initTs t)))
≡  {- fmap def -}
   Branch (fmap fold (initTs s)) (fmap fold (fmap (s `Branch`) (initTs t)))
≡  {- scanlT spec -}
  Branch (scanlT s) (fmap fold (fmap (s `Branch`) (initTs t)))
≡  {- functor law -}
  Branch (scanlT s) (fmap (fold ∘ (s `Branch`)) (initTs t))
≡  {- rework as λ -}
  Branch (scanlT s) (fmap (λ t' → fold (s `Branch` t')) (initTs t))
≡  {- fold def -}
  Branch (scanlT s) (fmap (λ t' → fold s ⊕ fold t')) (initTs t))
≡  {- rework λ -}
  Branch (scanlT s) (fmap ((fold s ⊕) ∘ fold) (initTs t))
≡  {- functor law -}
  Branch (scanlT s) (fmap (fold s ⊕) (fmap fold (initTs t)))
≡  {- scanlT spec -}
  Branch (scanlT s) (fmap (fold s ⊕) (scanlT t))
≡  {- lastT ∘ scanlT ≡ fold -}
  Branch (scanlT s) (fmap (lastT (scanlT s) ⊕) (scanlT t))
≡  {- factor out defs -}
  Branch s' (fmap (lastT s' ⊕) t')
     where s' = scanlT s
           t' = scanlT t

Suffix scan has a similar derivation.

  scanrT (Leaf a)
≡  {- scanrT def -}
  fmap fold (tailTs (Leaf a))
≡  {- tailTs def -}
  fmap fold (Leaf (Leaf a))
≡  {- fmap on T -}
  Leaf (fold (Leaf a))
≡  {- fold def -}
  Leaf a

  scanrT (s `Branch` t)
≡  {- scanrT spec -}
  fmap fold (tailTs (s `Branch` t))
≡  {- tailTs def -}
  fmap fold (Branch (fmap (`Branch` t) (tailTs s)) (tailTs t))
≡  {- fmap def -}
  Branch (fmap fold (fmap (`Branch` t) (tailTs s))) (fmap fold (tailTs t))
≡  {- scanrT spec -}
  Branch (fmap fold (fmap (`Branch` t) (tailTs s))) (scanrT t)
≡  {- functor law -}
  Branch (fmap (fold ∘ (`Branch` t)) (tailTs s)) (scanrT t)
≡  {- rework as λ -}
  Branch (fmap (λ s' → fold (s' `Branch` t)) (tailTs s)) (scanrT t)
≡  {- functor law -}
  Branch (fmap (λ s' → fold s' ⊕ fold t) (tailTs s)) (scanrT t)
≡  {- rework λ -}
  Branch (fmap ((⊕ fold t) ∘ fold) (tailTs s)) (scanrT t)
≡  {- scanrT spec -}
  Branch (fmap (⊕ fold t) (scanrT s)) (scanrT t)
≡  {- headT ∘ scanrT -}
  Branch (fmap (⊕ headT (scanrT t)) (scanrT s)) (scanrT t)
≡  {- factor out defs -}
  Branch (fmap (⊕ headT t') s') t'
    where s' = scanrT s
          t' = scanrT t

Extract code from these derivations:

scanlT' ∷ Monoid a ⇒ T a → T a
scanlT' (Leaf a)       = Leaf a
scanlT' (s `Branch` t) =
  Branch s' (fmap (lastT s' ⊕) t')
     where s' = scanlT' s
           t' = scanlT' t

scanrT' ∷ Monoid a ⇒ T a → T a
scanrT' (Leaf a)       = Leaf a
scanrT' (s `Branch` t) =
  Branch (fmap (⊕ headT t') s') t'
    where s' = scanrT' s
          t' = scanrT' t

Try it:

*T> t3
(("a","b"),("c",("d","e")))
*T> scanlT' t3
(("a","ab"),("abc",("abcd","abcde")))
*T> scanrT' t3
(("abcde","bcde"),("cde",("de","e")))

Efficiency

Although I was just following my nose, without trying to get anywhere in particular, this result is exactly the algorithm I first thought of when considering how to parallelize tree scanning.

Let’s now consider the running time of this algorithm. Assume that the tree is balanced, to maximize parallelism. (I think balancing is optimal for parallelism here, but I’m not certain.)

For a tree with $n$ leaves, the work $Wn$ will be constant when $n=1$ and $2\cdot W(n/2)+n$ when $n>1$. Using the Master Theorem (explained more here), $Wn=\Theta (n\log n)$.

This result is disappointing, since scanning can be done with linear work by threading a single accumulator while traversing the input tree and building up the output tree.

I’m using the term "work" instead of "time" here, since I’m not assuming sequential execution.

We have a parallel algorithm that performs $n\log n$ work, and a sequential program that performs linear work. Can we construct a linear-parallel algorithm?

Yes. Guy Blelloch came up with a clever linear-work parallel algorithm, which I’ll derive in another post.

Generalizing head and last

Can we replace the ad hoc (tree-specific) headT and lastT functions with general versions that work on all foldables? I’d want the generalization to also generalize the list functions head and last or, rather, to total variants (ones that cannot error due to empty list). For totality, provide a default value for when there are no elements.

headF, lastF ∷ Foldable f ⇒ a → f a → a

I also want these functions to be as efficient on lists as head and last and as efficient on trees as headT and lastT.

The First and Last monoids provide left-biased and right-biased choice. They’re implemented as newtype wrappers around Maybe:

newtype First a = First { getFirst ∷ Maybe a }

instance Monoid (First a) where
  ∅ = First Nothing
  r@(First (Just _)) ⊕ _ = r
  First Nothing      ⊕ r = r

newtype Last a = Last { getLast ∷ Maybe a }

instance Monoid (Last a) where
  ∅ = Last Nothing
  _ ⊕ r@(Last (Just _)) = r
  r ⊕ Last Nothing      = r

For headF, embed all of the elements into the First monoid (via First ∘ Just), fold over the result, and extract the result, using the provided default value in case there are no elements. Similarly for lastF.

headF dflt = fromMaybe dflt ∘ getFirst ∘ foldMap (First ∘ Just)
lastF dflt = fromMaybe dflt ∘ getLast  ∘ foldMap (Last  ∘ Just)

For instance,

*T> headF 3 [1,2,4,8]
1
*T> headF 3 []
3

When our elements belong to a monoid, we can use ∅ as the default:

headFM ∷ (Foldable f, Monoid m) ⇒ f m → m
headFM = headF ∅

lastFM ∷ (Foldable f, Monoid m) ⇒ f m → m
lastFM = headF ∅

For instance,

*T> lastFM ([] ∷ [String])
""

Using headFM and lastFM in place of headT and lastT, we can easily handle addition of an Empty case to our tree functor in this post. The key choice is that fold Empty ≡ ∅ and fmap _ Empty ≡ Empty. Then headFM will choose the first leaf, and lastT

What about efficiency? Because headF and lastF are defined via foldMap, which is a composition of fold and fmap, one might think that we have to traverse the entire structure when used with functors like [] or T.

Laziness saves us, however, and we can even extract the head of an infinite list or a partially defined one. For instance,

  foldMap (First ∘ Just) [5 ‥]
≡ foldMap (First ∘ Just) (5 : [6 ‥])
≡ First (Just 5) ⊕ foldMap (First ∘ Just) [6 ‥]
≡ First (Just 5)

So

  headF d [5 ‥]
≡ fromMaybe d (getFirst (foldMap (First ∘ Just) [5 ‥]))
≡ fromMaybe d (getFirst (First (Just 5)))
≡ fromMaybe d (Just 5)
≡ 5

And, sure enough,

*T> foldMap (First ∘ Just) [5 ‥]
First {getFirst = Just 5}
*T> headF ⊥ [5 ‥]
5

Where to go from here?

• As mentioned above, the derived scanning implementations perform asymtotically more work than necessary. Future posts explore how to derive parallel-friendly, linear-work algorithms. Then we’ll see how to transform the parallel-friendly algorithms so that they work destructively, overwriting their input as they go, and hence suitably for execution entirely in CUDA or OpenCL.
• The functions initTs and tailTs are still tree-specific. To generalize the specification and derivation of list and tree scanning, find a way to generalize these two functions. The types of initTs and tailTs fit with the duplicate method on comonads. Moreover, tails is the usual definition of duplicate on lists, and I think inits would be extend for "snoc lists". For trees, however, I don’t think the correspondence holds. Am I missing something?
• In particular, I want to extend the derivation to depth-typed, perfectly balanced trees, of the sort I played with in A trie for length-typed vectors and From tries to trees. The functions initTs and tailTs make unbalanced trees out of balanced ones, so I don’t know how to adapt the specifications given here to the setting of depth-typed balanced trees. Maybe I could just fill up the to-be-ignored elements with ∅.

I’ve been playing with deriving efficient parallel, imperative implementations of "prefix sum" or more generally "left scan". Following posts will explore the parallel & imperative derivations, but as a warm-up, I’ll tackle the functional & sequential case here.

foldl ∷ (b → a → b) → b → [a] → b
foldl f z []     = z
foldl f z (x:xs) = foldl f (z `f` x) xs

foldr ∷ (a → b → b) → b → [a] → b
foldr f z []     = z
foldr f z (x:xs) = x `f` foldr f z xs

scanl ∷ (b → a → b) → b → [a] → [b]

scanl f z [x1, x2,  ⋯ ] ≡ [z, z `f` x1, (z `f` x1) `f` x2, ⋯]

scanr ∷ (a → b → b) → b → [a] → [b]

scanr f z [⋯, xn_1, xn] ≡ [⋯, xn_1 `f` (xn `f` z), xn `f` z, z]

last (scanl f z xs) ≡ foldl f z xs
head (scanr f z xs) ≡ foldr f z xs

last ∘ scanl f z ≡ foldl f z
head ∘ scanr f z ≡ foldr f z

scanl ∷ (b → a → b) → b → [a] → [b]
scanl f z ls = z : (case ls of
                     []   → []
                     x:xs → scanl f (z `f` x) xs)

scanr ∷ (a → b → b) → b → [a] → [b]
scanr _ z []     = [z]
scanr f z (x:xs) = (x `f` q) : qs
                   where qs@(q:_) = scanr f z xs

scanl f z xs = map (foldl f z) (inits xs)
scanr f z xs = map (foldr f z) (tails xs)

scanl f z = map (foldl f z) ∘ inits
scanr f z = map (foldr f z) ∘ tails

inits ∷ [a] → [[a]]
inits []     = [[]]
inits (x:xs) = [[]] ++ map (x:) (inits xs)

tails ∷ [a] → [[a]]
tails []         = [[]]
tails xs@(_:xs') = xs : tails xs'

inits xs = [] : case xs of
                  []     → []
                  (x:xs) → map (x:) (inits xs)

last ∘ scanl f z ≡ foldl f z
head ∘ scanr f z ≡ foldr f z

  scanl f z []
≡ {- scanl spec -}
  map (foldl f z) (inits [])
≡ {- inits def -}
  map (foldl f z) [[]]
≡ {- map def -}
  [foldl f z []]
≡ {- foldl def -}
  [z]

  scanl f z (x:xs)
≡ {- scanl spec -}
  map (foldl f z) (inits (x:xs))
≡ {- inits def -}
  map (foldl f z) ([] : map (x:) (inits xs))
≡ {- map def -}
  foldl f z [] : map (foldl f z) (map (x:) (inits xs))
≡ {- foldl def -}
  z : map (foldl f z) (map (x:) (inits xs))
≡ {- map g ∘ map f ≡ map (g ∘ f) -}
  z : map (foldl f z ∘ (x:)) (inits xs)
≡ {- (∘) def -}
  z : map (λ ys → foldl f z (x:ys)) (inits xs)
≡ {- foldl def -}
  z : map (λ ys → foldl f (z `f` x) ys)) (inits xs)
≡ {- η reduction -}
  z : map (foldl f (z `f` x))) (inits xs)
≡ {- scanl spec -}
  z : scanl f (z `f` x) xs

  scanr f z []
≡ {- scanr spec -}
  map (foldr f z) (tails [])
≡ {- tails def -}
  map (foldr f z) ([[]])
≡ {- map def -}
  [foldr f z []]
≡ {- foldr def -}
  [z]

  scanr f z (x:xs)
≡ {- scanr spec -}
  map (foldr f z) (tails (x:xs))
≡ {- tails def -}
  map (foldr f z) ((x:xs) : tails xs)
≡ {- map def -}
  foldr f z (x:xs) : map (foldr f z) (tails xs)
≡ {- scanr spec -}
  foldr f z (x:xs) : scanr f z xs
≡ {- foldr def -}
  (x `f` foldr f z xs) : scanr f z xs
≡ {- head/scanr property -}
  (x `f` head (scanr f z xs)) : scanr f z xs
≡ {- factor out shared expression -}
  (x `f` head qs) : qs where qs = scanr f z xs
≡ {- stylistic variation -}
  (x `f` q) : qs where qs@(q:_) = scanr f z xs

This post is the last of a series of six relating numbers, vectors, and trees, revolving around the themes of static size-typing and memo tries. We’ve seen that length-typed vectors form a trie for bounded numbers, and can handily represent numbers as well. We’ve also seen that n-dimensional vectors themselves have an elegant trie, which is the n-ary composition of the element type’s trie functor:

type VTrie n a = Trie a :^ n 

where for any functor f and natural number type n,

f :^ n ≅ f ∘ ⋯ ∘ f  -- (n times)

This final post in the series places this elegant mechanism of n-ary functor composition into a familiar & useful context, namely trees. Again, type-encoded Peano numbers are central. Just as BNat uses these number types to (statically) bound natural numbers (e.g., for a vector index or a numerical digit), and Vec uses number types to capture vector length, we'll next use number types to capture tree depth.

Edits:

• 2011-02-02: Changes thanks to comments from Sebastian Fischer
  • Added note about number representations and leading zeros (without size-typing).
  • Added pointer to Memoizing polymorphic functions via unmemoization for derivation of Tree d a ≅ [d] → a.
  • Fixed signatures for some Branch variants, bringing type parameter a into parens.
  • Clarification about number of VecTree vs pairing constructors in remarks on left- vs right-folded trees.
• 2011-02-06: Fixed link to From Fast Exponentiation to Square Matrices.

data BinTree a = BinTree a (BinTree a) (BinTree a)

data BinNat = Zero | Even BinNat | Odd BinNat

data BinTree a = BinTree a (Pair (BinTree a))

data Pair a = Pair a a

type Pair = Id × Id

data BinNat = Zero | NonZero Bool BinNat

type BinNat = [Bool]

data Tree d a = Tree a (d ↛ Tree d a)

Tree d a ≅ [d] → a

data Vec ∷ * → * → * where
  ZVec ∷                Vec Z     a
  (:<) ∷ a → Vec n a → Vec (S n) a

data BinTree ∷ * → * → * where
  Leaf   ∷ a                         → BinTree Z     a
  Branch ∷ BinTree n a → BinTree n a → BinTree (S n) a

  Branch ∷ Pair (BinTree n) a → BinTree (S n) a

data BinTree ∷ * → * → * where
  Leaf   ∷                    a → BinTree Z     a
  Branch ∷ (Bool ↛ BinTree n a) → BinTree (S n) a

  Branch ∷ (BNat TwoT ↛ BinTree n a) → BinTree (S n) a

data VecTree ∷ * → * → * → * where
  Leaf   ∷                        a → VecTree b Z     a
  Branch ∷ (BNat b ↛ VecTree b n a) → VecTree b (S n) a

  Branch ∷ Vec b (VecTree b n a) → VecTree b (S n) a

data FTree ∷ (* → *) → * → * → * where
  Leaf   ∷               a → FTree f Z     a
  Branch ∷ f (FTree f n) a → FTree f (S n) a

type VecTree b = FTree (Vec b)

type TrieTree i = FTree (Trie i)

type VecTree b = TrieTree (BNat b)

data (:^) ∷ (* → *) → * → (* → *) where
  ZeroC ∷             a  → (f :^ Z    ) a
  SuccC ∷ f ((f :^ n) a) → (f :^ (S n)) a

type TrieTree i = (:^) (Trie i)

data BinTree ∷ * → * → * where
  Leaf   ∷                  a → BinTree Z     a
  Branch ∷ Pair (BinTree n) a → BinTree (S n) a

data BinTree ∷ * → * → * where
  Leaf   ∷                         a → BinTree Z     a
  Branch ∷ (BNat TwoT ↛ BinTree n a) → BinTree (S n) a

  Branch ∷ BinTree n (Pair a) → BinTree (S n) a

  Branch ∷ BinTree n (BNat TwoT ↛ a) → BinTree (S n) a

data VecTree ∷ * → * → * → * where
  Leaf   ∷                        a → VecTree b Z     a
  Branch ∷ (BNat b ↛ VecTree b n a) → VecTree b (S n) a

data VecTree ∷ * → * → * → * where
  Leaf   ∷                    a  → VecTree b Z     a
  Branch ∷ Vec b (VecTree b n a) → VecTree b (S n) a

  Branch ∷ VecTree b n (BNat b ↛ a) → VecTree b (S n) a

  Branch ∷ VecTree b n (Vec b a) → VecTree b (S n) a

Branch (Branch (Leaf 0, Leaf 1), Branch (Leaf 2, Leaf 3))

Branch (Branch (Leaf ((0,1),(2,3))))

data BinTree a = BinTree a (BinTree a) (BinTree a)

data BinTree a = Leaf a | Branch (Pair (BinTree a))

data BinTree a = Leaf a | Branch (BinTree (Pair a))

As you might have noticed, I’ve been thinking and writing about memo tries lately. I don’t mean to; they just keep coming up.

Memoization is the conversion of functions to data structures. A simple, elegant, and purely functional form of memoization comes from applying three common type isomorphisms, which also correspond to three laws of exponents, familiar from high school math, as noted by Ralf Hinze in his paper Generalizing Generalized Tries.

In Haskell, one can neatly formulate memo tries via an associated functor, Trie, with a convenient synonym "k ↛ v" for Trie k v, as in Elegant memoization with higher-order types. (Note that I've changed my pretty-printing from "k :→: v" to "k ↛ v".) The key property is that the data structure encodes (is isomorphic to) a function, i.e.,

k ↛ a ≅ k → a

In most cases, we ignore non-strictness, though there is a delightful solution for memoizing non-strict functions correctly.

My previous four posts explored use of types to statically bound numbers and to determine lengths of vectors.

Just as (infinite-only) streams are the natural trie for unary natural numbers, we saw in Reverse-engineering length-typed vectors that length-typed vectors (one-dimensional arrays) are the natural trie for statically bounded natural numbers.

BNat n ↛ a ≡ Vec n a

and so

BNat n → a ≅ Vec n a

In retrospect, this relationship is completely unsurprising, since a vector of length n is a collection of values, indexed by 0, . . . , n - 1.

In that same post, I noted that vectors are not only a trie for bounded numbers, but when the elements are also bounded numbers, the vectors can also be thought of as numbers. Both the number of digits and the number base are captured statically, in types:

type Digits n b = Vec n (BNat b)

The type parameters n and b here are type-encodigs of unary numbers, i.e., built up from zero and successor (Z and S). For instance, when b ≡ S (S Z), we have n-bit binary numbers.

In this new post, I look at another question of tries and vectors. Given that Vec n is the trie for BNat n, is there also a trie for Vec n?

Edits:

• 2011-01-31: Switched trie notation to "k ↛ v" to avoid missing character on iPad.

Vec n a → b ≅ (a × ⋯ × a) → b
             ≅ a → ⋯ → a → b
             ≅ a ↛ ⋯ ↛ a ↛ b
             ≡ Trie a (⋯ (Trie a b)⋯)
             ≅ (Trie a ∘ ⋯ ∘ Trie a) b

type VTrie n a = Trie a :^ n 

f :^ n ≅ f ∘ ⋯ ∘ f  -- (n times)

f :^ Z   ≅ Id
f :^ S n ≅ f ∘ (f :^ n)

data (:^) ∷ (* → *) → * → (* → *) where
  ZeroC ∷                        a  → (f :^ Z) a
  SuccC ∷ IsNat n ⇒ f ((f :^ n) a) → (f :^ (S n)) a

instance Functor f ⇒ Functor (f :^ n) where
  fmap h (ZeroC a)  = ZeroC (h a)
  fmap h (SuccC fs) = SuccC ((fmap∘fmap) h fs)

instance (IsNat n, Applicative f) ⇒ Applicative (f :^ n) where
  pure = pureN nat
  ZeroC f  ⊛ ZeroC x  = ZeroC (f x)
  SuccC fs ⊛ SuccC xs = SuccC (liftA2 (⊛) fs xs)

pureN ∷ Applicative f ⇒ Nat n → a → (f :^ n) a
pureN Zero     a = ZeroC a
pureN (Succ _) a = SuccC ((pure ∘ pure) a)

pureN (Succ n) a = SuccC ((pure ∘ pureN n) a)

  fmap = inZeroC  ($) ⊔ inSuccC  (fmap∘fmap)

  (⊛)  = inZeroC2 ($) ⊔ inSuccC2 (liftA2 (⊛))

inZeroC  h (ZeroC a ) = ZeroC (h a )

inSuccC  h (SuccC as) = SuccC (h as)

inZeroC2 h (ZeroC a ) (ZeroC b ) = ZeroC (h a  b )

inSuccC2 h (SuccC as) (SuccC bs) = SuccC (h as bs)

f :^ Z   ≅ Id
f :^ S n ≅ (f :^ n) ∘ f

data (:^) ∷ (* → *) → * → (* → *) where
  ZeroC ∷                        a  → (f :^ Z) a
  SuccC ∷ IsNat n ⇒ (f :^ n) (f a) → (f :^ (S n)) a

instance (IsNat n, HasTrie a) ⇒ HasTrie (Vec n a) where
  type Trie (Vec n a) = Trie a :^ n

  ZeroC v `untrie` ZVec      = v
  SuccC t `untrie` (a :< as) = (t `untrie` a) `untrie` as

  trie = trieN nat

trieN ∷ HasTrie a ⇒ Nat n → (Vec n a → b) → (Trie a :^ n) b
trieN Zero     f = ZeroC (f ZVec)
trieN (Succ _) f = SuccC (trie (λ a → trie (f ∘ (a :<))))

instance (IsNat n, HasTrie a) ⇒ HasTrie (Vec n a) where
  type Trie (Vec n a) = Trie a :^ n

  ZeroC b `untrie` ZVec      = b
  SuccC t `untrie` (a :< as) = (t `untrie` as) `untrie` a
  
  trie = trieN nat

trieN ∷ HasTrie a ⇒ Nat n → (Vec n a → b) → (Trie a :^ n) b
trieN Zero     f = ZeroC (f ZVec)
trieN (Succ _) f = SuccC (trie (λ as → trie (f ∘ (:< as))))

Vec n a ≅ a × ⋯ × a
        ≅ (1 → a) × ⋯ × (1 → a)
        ≅ (1 + ⋯ + 1) → a

data BNat ∷ * → * where
  BZero ∷          BNat (S n)
  BSucc ∷ BNat n → BNat (S n)

data Vec ∷ * → * → * where
  ZVec ∷                Vec Z     a
  (:<) ∷ a → Vec n a → Vec (S n) a

data BNat ∷ * → * where
  BSucc ∷ BNat n → BNat (S n)
  BZero ∷          BNat (S n)

data Vec ∷ * → * → * where
  (:<) ∷ Vec n a → a → Vec (S n) a
  ZVec ∷                Vec Z     a

The last few posts posts followed a winding path toward a formulation of a type for length-typed vectors. In Fixing lists, I mused how something like lists could be a trie type. The Stream functor (necessarily infinite lists) is the natural trie for Peano numbers. The standard list functor [] (possibly finite lists) doesn't seem to be a trie, since it's built from sums. However, the functor Vec n of vectors ("fixed lists") of length n is built from (isomorphic to) products only (for any given n), and so might well be a trie.

Of what type is Vec n the corresponding trie? In other words, for what type q is Vec n a isomorphic to q → a (for all a).

Turning this question on its head, what simpler type gives rise to length-typed vectors in a standard fashion?

Edits:

• 2011-02-01: Define Digits n b as BNat n ↛ BNat b.

Vec n a ≅ a × ⋯ × a
        ≅ (1 → a) × ⋯ × (1 → a)
        ≅ (1 + ⋯ + 1) → a

data BNat ∷ * → * where
  BZero ∷          BNat (S n)
  BSucc ∷ BNat n → BNat (S n)

type Digits n b = Vec n (BNat b)  -- n-digit number in base b

type Digits n b = BNat n ↛ BNat b

littleEndianToZ, bigEndianToZ ∷ ∀ n b. IsNat b ⇒ Digits n b → Integer

littleEndianToZ = foldr' (λ x s → fromBNat x + b * s) 0
 where b = natToZ (nat ∷ Nat b)

bigEndianToZ    = foldl' (λ s x → fromBNat x + b * s) 0
 where b = natToZ (nat ∷ Nat b)

*Vec> let ds = map toBNat [3,5,7] ∷ [BNat TenT]
*Vec> let v3 = fromList ds ∷ Vec ThreeT (BNat TenT)
*Vec> v3
fromList [3,5,7]
*Vec> littleEndianToZ v3
753
*Vec> bigEndianToZ v3
357

The post Fixing lists defined a (commonly used) type of vectors, whose lengths are determined statically, by type. In Vec n a, the length is n, and the elements have type a, where n is a type-encoded unary number, built up from zero and successor (Z and S).

infixr 5 :<

data Vec ∷ * → * → * where
  ZVec ∷                Vec Z     a
  (:<) ∷ a → Vec n a → Vec (S n) a

It was fairly easy to define foldr for a Foldable instance, fmap for Functor, and (⊛) for Applicative. Completing the Applicative instance is tricky, however. Unlike foldr, fmap, and (⊛), pure doesn't have a vector structure to crawl over. It must create just the right structure anyway. I left this challenge as a question to amuse readers. In this post, I give a few solutions, including my current favorite.

You can find the code for this post and the two previous ones in a code repository.

instance Functor (Vec n) where
  fmap _ ZVec     = ZVec
  fmap f (a :< u) = f a :< fmap f u

instance Applicative (Vec n) where
  pure a = ??
  ZVec      ⊛ ZVec      = ZVec
  (f :< fs) ⊛ (x :< xs) = f x :< (fs ⊛ xs)

  pure ∷ a → Vec n a

fmap ∷ (a → b) → Vec n a → Vec n b

class IsNat n

instance IsNat Z
instance IsNat n ⇒ IsNat (S n)

instance IsNat n ⇒ Applicative (Vec n) where
  ⋯

data Vec ∷ * → * → * where
  ZVec ∷ Vec Z a
  (:<) ∷ IsNat n ⇒ a → Vec n a → Vec (S n) a

instance IsNat n ⇒ Applicative (Vec n) where
  pure = ???
  (⊛)  = applyV

applyV ∷ Vec n (a → b) → Vec n a → Vec n b
ZVec      `applyV` ZVec      = ZVec
(f :< fs) `applyV` (x :< xs) = f x :< (fs `applyV` xs)

class IsNat n where pureN ∷ a → Vec n a

instance IsNat Z                where pureN a = ZVec
instance IsNat n ⇒ IsNat (S n) where pureN a = a :< pureN a

class IsNat n where units ∷ Vec n ()

instance            IsNat Z     where units = ZVec
instance IsNat n ⇒ IsNat (S n) where units = () :< units