Find a point's position on the integral Hilbert curve

Question

The Hilbert curves are a class of space-filling fractal curves where every bend in the path is at a right angle and the curve as a whole fills up a square, with the property that sequences of consecutive segments are always displayed as contiguous blocks.

Traditionally, the curve is drawn as a series of very squiggly lines inside a square such that each segment becomes smaller and smaller until they fill up the whole square. However, like the Walsh matrix, it is also possible to represent the limit curve as a sequence of line segments connecting lattice points that grows infinitely until it fills up an entire quadrant of lattice points on the Cartesian plane.

Suppose we define such a curve as the integral Hilbert curve, as follows:

• Start at the origin.

• Move one unit right (R), one unit up (U), and one unit left (L). Down will be notated as D.

• Let RUL be H_1. Then define H_2 as H_1 U I_1 R I_1 D J_1 where:

  • I_1 is H_1 rotated 90 degrees counterclockwise and flipped horizontally.

  • J_1 is H_1 rotated 180 degrees.

  Perform the steps of H_2 after H_1 to complete H_2.

• Define H_3 as H_2 R I_2 U I_2 L J_2, where:

  • I_2 is H_2 rotated 90 degrees clockwise and flipped vertically.

  • J_2 is H_2 rotated 180 degrees.

  Perform the steps of H_3 after H_2 to complete H_3.

• Define H_4 and other H_n for even n the same way as H_2, but relative to H_{n-1}.

• Define H_5 and other H_n for odd n the same way as H_3, but relative to H_{n-1}.

• The integral Hilbert curve is H_infinity.

We get a function f(n) = (a, b) in this way, where (a, b) is the position of a point moving n units along the integral Hilbert curve starting from the origin. You may notice that if we draw this curve 2^{2n}-1 units long we get the nth iteration of the unit-square Hilbert curve magnified 2^n times.

Your task is to implement f(n). The shortest code to do so in any language wins.

---

Example inputs and outputs:

> 0
(0, 0)
> 1
(1, 0)
> 9
(2, 3)
> 17
(4, 1)
> 27
(6, 3)
> 41
(7, 6)
> 228
(9, 3)
> 325
(24, 3)

Note that your solution's algorithm must work in general; if you choose to implement your algorithm in a language where integers are of bounded size, that bound must not be a factor in your program's code.

Answer 1

GolfScript (62 bytes)

{1:a:b(.@4base{:?3%!:|b^?2<3?=.b^:b^^|!.a^:a^{\}*{@2*@+}2*}/\}

This is an anonymous function which takes n on the stack and leaves x y on the stack. Online demo

It's essentially the iterative code from the Wikipedia article reversed so that it processes the base-4 representation of n from most-significant to least-significant digit instead of least-significant to most-significant. This simplifies the bookkeeping required to output x and y in the right order.

Pseudocode translation (not tested, may contain errors):

a = b = true;
x = y = 0;
foreach (digit in n base 4) {
    p = (digit % 3 == 0) ^ b;
    q = (digit < 2) ^ b;
    b ^= (digit == 3);
    if (a) swap(p, q);
    a ^= (digit % 3 != 0);
    y = 2*y + (p ? 1 : 0);
    x = 2*x + (q ? 1 : 0);
    swap(x, y);
}
swap(x, y)

Answer 2

R: 178 169

Room for improvement in this, but it's a start

f=function(d){x=y=0;s=1;while(s<Inf){a=bitwAnd(1,d/2);b=bitwAnd(1,bitwXor(d,a));if(b==0){if(a==1){x=s-1-x;y=s-1-y};t=x;x=y;y=t};x=x+(s*a);y=y+(s*b);d=d/4;s=s*2};c(x,y)}

Test Run

> f(0)
[1] 0 0
> f(1)
[1] 1 0
> f(9)
[1] 2 3
> f(41)
[1] 7 6
> f(325)
[1] 24  3