http://people.csail.mit.edu/jaffer/Geometry/RFMDSFF  
Recursive Formulation of Multidimensional Hilbert SpaceFilling Curves 
The rigorous explanations of algorithms for multidimensional spacefilling functions tend to be cryptic and difficult to reconcile with their simple geometric interpretations. This article attempts to rectify that.
The properties we desire for spacefilling transforms are that they be:
Because the transforms are selfsimilar, the scalar can be considered to be composed of digits in some base, where digits correspond to positions in nested shapes. Because transforms are measurepreserving, the largest copies of the selfsimilar shape correspond to the highestorder digits of the scalar.
Any real scalar has a most significant digit; because the leastsignficant 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 highestorder digits; and information must pass from higherorder to lowerorder 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 tailrecursive, which is equalent to the interative loops in the extant algorithms.
A Gray code is an ordering of nonnegative integers in which exactly one bit differs between each pair of successive elements. There are multiple Gray codings. An nbit Gray code maps the integers from 0 to 2^{n} to a Hamiltonian cycle on an ndimensional hypercube. Graycodes are thus used to generate offsets at each hypercube level.
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 rotatebitfield
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 Graycode sequence. The
log2binaryfactors
machinations appear to compute the
additional rotation based on the next bit to be changed in the
Graycode sequence.
log2binaryfactors
is a simple but nonobvious function.
I encountered it in Doug Moore's
"hilbert.c"
([Moore 1999]). Here is the
scheme version:
(define (log2binaryfactors n) (+ 1 (integerlength (logand n ( n)))))
The log2binaryfactors
function returns the number of
factorsof2 of chnk, which is also the index of the
lowestorder 1 bit in chnk. But more importantly for this
application, log2binaryfactors
is related to the
nextbitchanged in incrementing Graycodes.
For Graycodes with odd numbers of 1 bits (which correspond to odd
integers), log2binaryfactors
returns one less than the
index of the bit different in the next Graycode in the sequence.
For Graycodes with even numbers of bits (which correspond to even
integers), log2binaryfactors
returns one less than
the index of the bit different in the previous Graycode in the
sequence.
sequence  Gray code  log2binaryfactors + 2  sequence  Gray code  log2binaryfactors + 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 loworder bit.
(define (integer>hilbertcoordinates scalar rank . nbits) (define igry (integer>graycode scalar)) (define rnkmsk (lognot (ash 1 rank))) (define rnkhib (ash 1 (+ 1 rank))) (define rank*nbits (if (null? nbits) (let ((rank^2 (* rank rank))) (* (quotient (+ 1 rank^2 (integerlength scalar)) rank^2) rank^2)) (* rank (car nbits)))) (define (loop bdxn chnk rotation flipbit lst) (if (positive? bdxn) (map graycode>integer (delaminatelist rank (reverse lst))) (loop (+ rank bdxn) (logxor rnkhib (logand (ash igry (+ rank bdxn))) rnkmsk) (modulo (+ rotation (log2binaryfactors chnk) 2) rank) (ash 1 rotation) (cons (logxor flipbit (rotatebitfield chnk rotation 0 rank)) lst)))) (loop ( rank rank*nbits) (ash igry ( rank rank*nbits)) 0 0 '())) 
integer>hilbertcoordinates
takes an optional argument,
nbits.
Flipping the highorder bit of a Graycoded number is the equivalent
of taking the onescomplement of the uncoded number.
For 0 <= k < 2*rnkhib:
(= ( (* 2 rnkhib) k 1) (graycode>integer (logxor rnkhib (integer>graycode k))))All but the highest order bdxn get their highorder bit flipped. With flipbit, all bdxn have an even number of flipped bits.
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 rank^{2} bits are processed.
loop
is the recursive function.
If bdxn is positive, then all the digits have been processed; convert the Graycoded coordinates to regular binary representation and return them.
If bdxn is not positive, then call loop
with
arguments:
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.
rank*nbits is 25.
(integer>graycode #b0011100001011010101011000) ==> #b0010010001110111111110100 xor rnkhib ==> #b0010000001010110111100100The high bit of each rankbit 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 exclusiveored 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  (01100 11001 00110 10001 00100) 
From left to right, lst gets reversed (top to bottom), delaminated (matrix transpose), and Graydecoded.
lst  reverse  delaminate  graycode>integer  

 ==> 
 ==> 
 ==> 

(integer>hilbertcoordinates #b0011100001011010101011000 5) ==> (#b01100 #b00010 #b11001 #b00111 #b01100)
graycode>integer
; is this related to item 1?
log2binaryfactors
ignore changes in the low order bit
of the Graycoded chnk?
I am a guest and not a member of the MIT Computer Science and Artificial Intelligence Laboratory.
My actions and comments do not reflect in any way on MIT.  
SpaceFilling Curves  
agj @ alum.mit.edu  Go Figure! 