Static indices for speedy array-based algorithms

Posted on December 21, 2014 by Mike Ledger

GHC.TypeLits gives us some very powerful type-level functionality (with some caveats that will hopefully be worked out in future) that has not yet seen very widespread use. With some mildly ugly leg-work, I’ve made a small library that provides a n-dimensional statically bounded index type, and associated functions for their application with as minimal overhead as possible.

Motivation

Making plissken, I was unhappy with the state of linear algebra libraries on Hackage. The venerable hmatrix is excellent - once your input sizes outweigh the constant factors involved. For the matrices that my game engine lived and survived on – 4x4 matrices and 4-vectors – hmatrix was easily outperformed in simple 4x4 matrix multiplication by linear, which is a naiive Haskell implementation of the sort of “small” linear algebra that I need here. I ended up using hmatrix though, since it being built on Storable made for very easy interaction with OpenGL.

A trick I found in order to inline away “static recursion” was to abstract the use of the recursion parameter into a typeclass whose function is inlined. So long as you remember the all-powerful -O2 switch, GHC will happily inline away every recursive call, which will hopefully give way to further compile-time evaluation.

Here’s a little demonstration of such a class:

That’s pretty gross. We can improve it by wrapping it in a small helper:

Unfortunately GHC doesn’t see as we can intuitively, that since all Nat has an instance for ToPeano, For has the instance For (ToPeano (n :: Nat)), so we have to add that constraint.

Let’s try it out:

Some scrolling later…

Great! The loop is unrolled rather nicely.

Enter indices

indices provides a multi-dimensional (Int based) index type, with use for array-heavy code in mind. I hope Soon® to release on Hackage the package static (patches welcome and appreciated) which provides a ForiegnPtr based array type that leverages indices’s, ahem, indices, just about everywhere.

Indices in indices are these two types:

This is similar to the design seen in repa, except the “bound” of an index is its type. That’s crippling for code with arbitrarily-bounded arrays, but very nice otherwise. For instance, an index into a 4x4 matrix is 0:.0:.Z :: 4:.4:.Z .

For now, you can use indices for array-based code by leveraging its Ix instance. This is a little confusing due to the design of Ix, but it’s fairly simple: the type is always the upper bound, and zero is always the lower bound, not a value you give. That means that arrays constructed by array (_, _ :: t) are bounded by [0,t). As an aside, I’m leaning towards reworking static to simply use the array types found in arrays to simplify it greatly into the small linear-algebra package it yearns to be. Input is appreciated here. (At the moment I’m irked at having to store two superfluous lower/upper bounds – linked lists of Int – in the array constructors).

Here’s a demonstration, implementing vector dot product with a static, unrolled loop:

Note the slyness in me not writing the type signature to dot… Well, the important thing is that GHC happily infers the type without needing any hints. Right… guys?

You can contribute to indices on GitHub, find its documentation on Hackage, and install it with cabal.


blog comments powered by Disqus