/* Variation of lamxy that eliminates the branch in the loop. */
#include
#include
#include
/*
Given the "order" n of a Hilbert curve and a distance s along the curve,
this program computes the corresponding (x, y) coordinates. The square
that the Hilbert curve traverses is of size 2**n by 2**n.
The method is that given in [Lam&Shap], described by the following
table. Here i = n-1 for the most significant bit of x and y, and i = 0
for the least significant bits.
s[2i+1:2i] x[i] y[i] x[i-1:0] y[i-1:0]
-----------|-------------------------------
00 | 0 0 y[i-1:0] x[i-1:0]
01 | 0 1 x[i-1:0] y[i-1:0]
10 | 1 1 x[i-1:0] y[i-1:0]
11 | 1 0 ~y[i-1:0] ~x[i-1:0]
To use this table, start at the least significant two bits of s (i = 0).
If they are both 0 (first row), set the least significant bits of x and
y to 0 and 0 respectively, and interchange x and y. The last two
columns designate an interchange of the bits of x and y to the right of
bit i, but on this first iteration they are null, so there is no
interchange to do. If the least significant two bits of s are 10 (third
row), set the least significant bits of x and y to 1, and similarly for
the other rows.
Then, consider the next least significant two bits of s, and select
the appropriate row of the table to determine the next bits of x and y,
and how to change the bits of x and y to the right of i. Continue until
the most significant bits of x and y have been processed. */
// ------------------------------ cut ----------------------------------
void hil_xy_from_s(unsigned s, int n, unsigned *xp,
unsigned *yp) {
int i, sa, sb;
unsigned x, y, swap, cmpl;
for (i = 0; i < 2*n; i += 2) {
sa = (s >> (i+1)) & 1; // Get bit i+1 of s.
sb = (s >> i) & 1; // Get bit i of s.
swap = (sa ^ sb) - 1; // -1 if should swap, else 0.
cmpl = -(sa & sb); // -1 if should compl't, else 0.
x = x ^ y;
y = y ^ (x & swap) ^ cmpl;
x = x ^ y;
x = (x >> 1) | (sa << 31); // Prepend sa to x and
y = (y >> 1) | ((sa ^ sb) << 31); // (sa^sb) to y.
}
*xp = x >> (32 - n); // Right-adjust x and y
*yp = y >> (32 - n); // and return them to
} // the caller.
// ------------------------------ cut ----------------------------------
int main(int argc, char *argv[]) {
unsigned n, N, s, x, y;
if (argc <= 1) {
printf("Need one parameter, the order of the curve.\n");
exit(1);
}
n = atoi(argv[1]);
N = 1 << 2*n; // N = 2**2n.
printf("(x, y) coordinates along the Hilbert curve of order %d\n", n);
for (s = 0; s < N; s++) {
hil_xy_from_s(s, n, &x, &y); // Pass addresses of x and y.
printf("%5d %5d %5d\n", s, x, y);
}
return 0;
}