There is a surfeit of 2-dimensional space-filling curves, but
generalizing them to higher rank is not necessarily practical or even
possible. [Bongki 2001] attributes
generalization of space-filling curves to any rank to
[Butz 1969]; who gives algorithms
for the Hilbert and the Peano space-filling curves in separate papers.

The rigorous explanations of algorithms for multidimensional
space-filling functions tend to be cryptic and difficult to reconcile
with their simple geometric interpretations. This article attempts to
rectify that.

Properties

The properties we desire for space-filling transforms are that they
be:

space-filling

a continuous function whose range comes arbitrarily close to any
given point in the n-dimensional unit hypercube.

multidimensional

exist for all ranks greater than one.

self-similar

can be defined as a nonterminating recurrence.

Lebesgue measure-preserving.

equal length intervals map to equal volumes.

isotropic

the lengths of axis-aligned projections of the space-filling
curve should be equal.

Although there are many shapes which can tile the plane, the only
regular polyhedron which tiles
R3 is the cube. Isotropy rules
out parallelograms and their ilk. The general family of curves will
fill squares, cubes, and hypercubes.

Because the transforms are self-similar, the scalar can be considered
to be composed of digits in some base, where digits correspond to
positions in nested shapes. Because transforms are
measure-preserving, the largest copies of the self-similar shape
correspond to the highest-order digits of the scalar.

Any real scalar has a most significant digit; because the
least-signficant digits could be arbitrarily small, the same can't be
said of them. Thus the calculation of transformations between scalars
and multidimensional coordinates must start with the highest-order
digits; and information must pass from higher-order to lower-order
frames.

This suggests a recursive organization for conversion computations.
Were the task to sequentially enumerate points on the curve, then
recursive computation would have the familiar profile where the number
of calls increases exponentially with the depth of calls.
But that is only useful for drawing pictures. The other uses convert
or compare a few points sparsely scattered in the multidimensional
space. The computation will return a single value: a boolean, scalar
or vector (coordinates). This is still recursive; but tail-recursive,
which is equalent to the interative loops in the extant algorithms.

A Gray code is an ordering of non-negative integers in
which exactly one bit differs between each pair of successive
elements. There are multiple Gray codings. An n-bit Gray code
maps the integers from 0 to 2n
to a Hamiltonian cycle on an n-dimensional hypercube.
Gray-codes are thus used to generate offsets at each hypercube level.

Rotations

A cube's rotation is the combination of the rotations of all the cubes
it is part of. Thus the rotation is passed to the recursive function
and each call augments it.

Although any permutation which cycles all axes through all positions
could be used to "rotate" hypercubes in order to align them with the
corners they they are replacing, cyclic permuations are sufficient and
used here. The rotate-bit-field function is used to
rotate fields of rank bits.

But how much rotation should be applied? The coordinates of a vertex
alone do not tell which direction the next edge is to be drawn. That
direction is controlled by the Gray-code sequence. The
log2-binary-factors machinations appear to compute the
additional rotation based on the next bit to be changed in the
Gray-code sequence.

log2-binary-factors is a simple but non-obvious function.
I encountered it in Doug Moore's
"hilbert.c"
([Moore 1999]). Here is the
scheme version:

The log2-binary-factors function returns the number of
factors-of-2 of chnk, which is also the index of the
lowest-order 1 bit in chnk. But more importantly for this
application, log2-binary-factors is related to the
next-bit-changed in incrementing Gray-codes.

For Gray-codes with odd numbers of 1 bits (which correspond to odd
integers), log2-binary-factors returns one less than the
index of the bit different in the next Gray-code in the sequence.

For Gray-codes with even numbers of bits (which correspond to even
integers), log2-binary-factors returns one less than
the index of the bit different in the previous Gray-code in the
sequence.

sequence

Gray code

log2-binary-factors + 2

sequence

Gray code

log2-binary-factors + 2

0

00000

1

1

00001

2

2

00011

2

3

00010

3

4

00110

3

5

00111

2

6

00101

2

7

00100

4

8

01100

4

9

01101

2

10

01111

2

11

01110

3

12

01010

3

13

01011

2

14

01001

2

15

01000

5

16

11000

5

The other points in the sequence all are changes of the low-order bit.

All but the highest order bdxn get their high-order bit
flipped. With flipbit, all bdxn have an even
number of flipped bits.

rank*nbits is the number of bits to process from the
scalar. If the nbits argument is given, then it is just
nbits times rank; otherwise, the smallest
integer multiple of rank2
greater or equal to the number of bits necessary to represent
scalar.

The orientation of nested cubes precesses through all
rank axes. To guarantee match of orientation for varied
precision inputs, a multiple of rank cubes is processed.
There are rank bits in each cube; hence
rank2 bits are processed.

loop is the recursive function.

If bdxn is positive, then all the digits have been
processed; convert the Gray-coded coordinates to regular binary
representation and return them.

If bdxn is not positive, then call loop with
arguments:

bdxn is the (negative) number of bits to right shift
igry in order to bring the current set of rank bits to
lowest order.

chnk is the rank-bit Gray-coded number
currently being processed.

rotation is the number of bits to rotate the chunk by.
Notice that rotation accumulates the earlier rotations.

flipbit is the bit in the rotated chnk to
invert; it is the bit indexed by the calling function's
rotation. Note that while rotation is one
call delayed from the chnk which generated it,
flipbit is two calls delayed.

The list lst accumulates the results of processing
chnk.

The initial chnk is the highest order; it needs no
masking.

But if the rank (high) bit of the highest order group is
flipped, then zero maps to the other end of the curve from the
origin. If the initial flipbit is 1, then zero maps to a
point other than the origin.

The high bit of each rank-bit group except the highest order group are
flipped (by rnkhib) and are underlined. The leftmost
number in lst is derived from the chnk rotated
by rotation, then exclusive-ored with the
flipbit of the preceeding row.

bxdn

chnk

rotation

flipbit

lst

-20

00100

0

00000

()

-15

00001

4

00001

(00100)

-10

01011

1

10000

(10001 00100)

-5

01111

3

00010

(00110 10001 00100)

0

00100

0

01000

(11001 00110 10001 00100)

5

10000

4

00001

(01100 11001 00110 10001 00100)

From left to right, lst gets reversed (top to bottom),
delaminated (matrix transpose), and Gray-decoded.