Hilbert Curves

The discrete Hilbert curve Hn of level n can be viewed as a function from a single 2n-bit unsigned integer to two n-bit unsigned integers. The straightforward method for computing a single point (x,y) = H(t) on this curve requires O(n) operations to compute one bit of the result at a time (see the Wikipedia article for details).

This can be sped up to O(log n) using parallel prefix ideas as follows:

  1. Uninterleave the bit string t into two n-bit strings (tx,ty), which takes O(log n) operations.
  2. Note that each pair of bits in the result (x,y) is given by a applying one of the isometries of the square to the corresponding pair of bits in (tx,ty). The particular isometry is determined by the product of the isometries produced by all higher bits.
  3. Use a parallel prefix circuit to compute the isometries for each bit pair. Conveniently, the isometries required form a subgroup of D4 isomophic to Z2 x Z2, so the group products are just xors. This takes O(log n) operations if n is within the word size of the machine.
  4. Compute all bits of the result using the isometries with O(1) operations.
The expensive parts of this are uninterleaving t and computing the parallel prefixes. Both of these operations can be incrementally updated in constant time, so this also gives an O(1) algorithm for incrementally walking along a Hilbert curve. Unfortunately, I do not know how to extend these ideas to compute the inverse of the Hilbert curve (i.e., t as a function of (x,y)).

The file hilbert.c contains random access and incremental versions of this algorithm, as well as testing code.

To see what one of these actually looks like, here's a Hilbert curve of level 8 mapping a smooth ramp from green to blue to red into the square:

Hilbert curve of level 8