A compact Hilbert curve routine

January 30, 2017

Here is a relatively compact Mathematica routine for generating the b-th iterate of a Hilbert curve in [0,1]^n. The algorithm is due to Skilling:

HilbertCurve =
     Compile[{{x, _Integer}, {b, _Integer}, {n, _Integer}},
             Module[{t = BitXor[x, Quotient[x, 2]], p = 2, k, q, r, xx},
                    xx = Total[Table[BitAnd[Quotient[t, 2^(r (n - 1) + k)], 2^r],
                                     {r, b - 1, 0, -1}, {k, n - 1, 0, -1}]];
                    Do[q = p - 1;
                       Do[t = xx[[k]];
                          If[BitAnd[t, p] != 0,
                             xx[[1]] = xx[[1]] ~BitXor~ q,
                             t = BitAnd[t ~BitXor~ xx[[1]], q];
                             xx[[{1, k}]] = BitXor[xx[[{1, k}]], t]],
                          {k, n, 2, -1}];
                       If[BitAnd[xx[[1]], p] != 0,
                          xx[[1]] = xx[[1]] ~BitXor~ q];
                       p *= 2, {r, b - 1}];
                    2 xx/(p - 1) - 1], RuntimeAttributes -> {Listable}]

I’ve seen a variety of routines for generating the Hilbert curve before, but none were as compact or as general as this one. If you want a version that generates exact rational coordinates, it is easy to take out the requisite parts from the Compile[] function.

Here are the 2D and 3D versions of the curve:

With[{p = 6}, 
     Graphics[{ColorData[97, 1],
               Line[HilbertCurve[Range[0, 2^(2 p) - 1], p, 2]]}, Frame -> True]]


With[{p = 3}, 
     Graphics3D[Tube[HilbertCurve[Range[0, 2^(3 p) - 1], p, 3]], Axes -> True]]


Here’s a perspective projection of a four-dimensional Hilbert curve:

With[{p = 3, k = 2, f = 4, d = 3/2, r = 1/8}, 
     Graphics3D[{Directive[EdgeForm[], CapForm[None], GrayLevel[1/5], 
                           Glow[ColorData["Legacy", "SteelBlue"]],
                           Specularity[ColorData["Legacy", "Gold"], 16]], 
                 Tube[(f Delete[#, k])/(d - Extract[#, k]) & /@ 
                      HilbertCurve[Range[0, 2^(4 p) - 1], p, 4], r]}, 
                Boxed -> False, Lighting -> "Neutral"]]


(This visualization is admittedly a bit messy-looking; if you have suggestions for making a nice display of the four-dimensional version, tell me about it in the comments!)


Some terms left unsimplified

January 6, 2017

I had been slowly going through some OEIS entries in an attempt to add or optimize the Mathematica code found there. One thing that struck me with respect to the sequences that are related to the elliptic integrals/elliptic functions was that they were often left in the form that was spit out by a CAS. I had railed about the deficiencies of current software for generating simple expressions for elliptic integrals many times before, so this was but a confluence with the fact that not many people these days have experience in the analytical manipulation of these functions, as classical and important as they are.

I certainly have my work cut out for me…