3D (e.g.) Hilbert is easy with the finite state machine. Three bits of t in, one bit of x, y, and z out, plus state change. --rwg --well... actually figure out the state machine for arbitrary D, and do it in a way that we do not need no stinking huge 2^D entry table to store the state-rules, and then I'm sold. There's some disgusting computer programs out there I found with google to do multi-D hilbert curves, but there really ought to be a nice way that takes under 1 page of code. And they haven't found it. --wdsmith OK, but first, a vid of Gareth Conway <http://daftmusings.com/wp-content/uploads/2011/06/G4G9-052.jpg> bouncing on a steel one: http://blog.makezine.com/archive/2010/11/math_monday_3d_hilbert_curve_in_ste... Here's some relevant old funmail. The three pngs at the end have as usual vanished from gosper.org, so I wrote a three page pdf of the 208 KB notebook containing them. 78.8MB!! gosper.org/Trilbertloopscenterscorners.pdf If Adobe Reader ever finishes "optimizing" it, I'll rewrite it. I could write a png, but then you couldn't copy the text. On 2011-01-05 14:39, Bill Gosper wrote:
Last year I sent one named Treano and posted some graphs. It had three problems: 1) As Neil Bickford points out, it should have been named Trilbert. 2) It spanned a face diagonal instead of an edge, creating inconsistent textures when sampled at fixed frequencies. 3) It was so obscure that when I found the edge-spanning construction, I couldn't decipher and modify my own code.
Fortunately, Corey Ziegler Hunts was up to both tasks: I3 = IdentityMatrix[3];
In[720]:= Vects = {{0, 0, 0}, {0, 0, 1}, {0, 1, 1}, {0, 2, 1}, {1, 2, 1}, {2, 2, 1}, {2, 1, 1}, {2, 0, 1}, {2, 0, 0}}/2;
Mats = {{{0, 0, 1}, {0, 1, 0}, {-1, 0, 0}}.{{-1, 0, 0}, {0, 1, 0}, {0, 0, 1}}, {{0, -1, 0}, {1, 0, 0}, {0, 0, 1}}.{{1, 0, 0}, {0, 0, -1}, {0, 1, 0}}, {{0, -1, 0}, {1, 0, 0}, {0, 0, 1}}.{{1, 0, 0}, {0, 0, -1}, {0, 1, 0}}, {{1, 0, 0}, {0, 0, -1}, {0, 1, 0}}.{{1, 0, 0}, {0, -1, 0}, {0, 0, 1}}, {{1, 0, 0}, {0, 0, -1}, {0, 1, 0}}.{{1, 0, 0}, {0, -1, 0}, {0, 0, 1}}, {{0, 1, 0}, {-1, 0, 0}, {0, 0, 1}}.{{1, 0, 0}, {0, 0, -1}, {0, 1, 0}}, {{0, 1, 0}, {-1, 0, 0}, {0, 0, 1}}.{{1, 0, 0}, {0, 0, -1}, {0, 1, 0}}, {{0, 0, -1}, {0, 1, 0}, {1, 0, 0}}.{{-1, 0, 0}, {0, 1, 0}, {0, 0, 1}}, {{1, 0, 0}, {0, 1, 0}, {0, 0, 1}}}/2;
In[722]:= Trilbert[t_, mat_: I3, vect_: {0, 0, 0}] := (Trilbert[t, b1_: I3, b0_: {0, 0, 0}] = (Trilbert[t, matend_: I3, vectend_: {0, 0, 0}] := Inverse[matend - mat].(vect - vectend); Module[{t8 = 8*t, n}, n = Floor[t8] + 1; Vects[[n]] + Mats[[n]].Trilbert[FractionalPart[t8], mat.Mats[[n]], vect + mat.Vects[[n]]]]))
(style-tweaked in ways Corey might disavow).
The inverse images of the eight vertices are
In[724]:= Trilbert /@ {0, 5/28, 9/28, 3/7, 4/7, 19/28, 23/28, 1}
Out[724]= {{0, 0, 0}, {0, 0, 1}, {0, 1, 1}, {0, 1, 0}, {1, 1, 0}, {1, 1, 1}, {1, 0, 1}, {1, 0, 0}}
(which you couldn't directly test without the exact rational trick). There are eight(!) inverse[s] of the center point:
In[729]:= {19/28, 1 + 23/28, 2 + 5/28, 3 + 23/28, 4 + 5/28, 5 + 23/28, 6 + 5/28, 7 + 9/28}/8
Out[729]= 19 51 61 107 117 163 173 205 {---, ---, ---, ---, ---, ---, ---, ---} 224 224 224 224 224 224 224 224
Trilbert /@ % Out[728]= 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 {{-, -, -}, {-, -, -}, {-, -, -}, {-, -, -}, {-, -, -}, {-, -, -}, {-, -, -}, {-, -, -}} 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
Every interior triple of dyadic rationals has eight inverses.
Plots: http://gosper.org/Trilbertloops.png http://gosper.org/Trilbertcenters.png http://gosper.org/Trilbertcorners.png --rwg (While I was typing this, a popup said "Mathematica quit unexpectedly". When I clicked Restart, it said my license had expired, demanding a password, which I requested upon the first warning Christmas Eve. Why does WRI synchronize this with everyone's vacation?) It turns out the unexpected quit was due to a 7.0 bug eating the entire Apple VM while trying to convert to a .pdf the small .nb containing these three plots.
This version of Trilbert used the exact rational trick. Corey's faster but lengthier dyadic rational definition: ClearAll[Trilbert]; Trilbert[x_, funcs_: Table[{0 &, 1, 1}, {i, 1, 9, 1}]] := Block[{vects, mats, reflect, vect = {0, 0, 0}, mat = IdentityMatrix[3], vectr, matr, digits, digit, nonrep = 0, rep, special, period, i}, vects = {{0, 0, 0}, {0, 0, 1/2}, {0, 1/2, 1/2}, {0, 1, 1/2}, {1/2, 1, 1/2}, {1, 1, 1/2}, {1, 1/2, 1/2}, {1, 0, 1/2}, {1, 0, 0}}; mats = Table[0, {i, 1, 9, 1}]; mats[[1]] = (1/ 2)*{{0, 1, 0}, {-1, 0, 0}, {0, 0, 1}}.{{0, 0, -1}, {0, 1, 0}, {1, 0, 0}}; mats[[2]] = (1/ 2)*{{0, -1, 0}, {1, 0, 0}, {0, 0, 1}}.{{1, 0, 0}, {0, 0, -1}, {0, 1, 0}}; mats[[3]] = (1/ 2)*{{0, -1, 0}, {1, 0, 0}, {0, 0, 1}}.{{1, 0, 0}, {0, 0, -1}, {0, 1, 0}}; mats[[4]] = (1/ 2)*{{1, 0, 0}, {0, 0, -1}, {0, 1, 0}}.{{1, 0, 0}, {0, 0, -1}, {0, 1, 0}}; mats[[5]] = (1/ 2)*{{1, 0, 0}, {0, 0, -1}, {0, 1, 0}}.{{1, 0, 0}, {0, 0, -1}, {0, 1, 0}}; mats[[6]] = (1/ 2)*{{0, 1, 0}, {-1, 0, 0}, {0, 0, 1}}.{{1, 0, 0}, {0, 0, -1}, {0, 1, 0}}; mats[[7]] = (1/ 2)*{{0, 1, 0}, {-1, 0, 0}, {0, 0, 1}}.{{1, 0, 0}, {0, 0, -1}, {0, 1, 0}}; mats[[8]] = (1/ 2)*{{0, -1, 0}, {1, 0, 0}, {0, 0, 1}}.{{0, 0, 1}, {0, 1, 0}, {-1, 0, 0}}; mats[[9]] = (1/2)*{{1, 0, 0}, {0, 1, 0}, {0, 0, 1}}; reflect = {{1, 0, 0}, {0, 0, 1}, {0, 1, 0}}; If[And[0 < x, x < 1], digits = RealDigits[x + 1, 8][[1, 2 ;; -1]] + 1, If[x == 0, digits = {1}, digits = {9} ] ]; mat = MatrixPower[reflect, funcs[[1, 1]][0]]; For[{i = 1}, And[i <= Length[digits], ! ListQ[digits[[i]]]], i++, digit = digits[[i]]; vect = vect + mat.vects[[digit]]; mat = mat.mats[[digit]].MatrixPower[reflect, funcs[[digit, 1]][i]]; nonrep++ ]; If[ListQ[digits[[-1]]], digits = digits[[-1]]; rep = Length[digits]; special = Max[funcs[[Union[digits], 2]]] - nonrep; For[{i = 1}, i <= special, i++, digit = digits[[Mod[i - 1, rep] + 1]]; vect = vect + mat.vects[[digit]]; mat = mat.mats[[digit]].MatrixPower[reflect, funcs[[digit, 1]][i + nonrep]] ]; If[special < 0, special = 0 ]; vectr = {0, 0, 0}; matr = IdentityMatrix[3]; period = Apply[LCM, Append[funcs[[Union[digits], 3]], Length[digits]]]; For[{i = 1}, i <= period, i++, digit = digits[[Mod[i + special - 1, rep] + 1]]; vectr = vectr + matr.vects[[digit]]; matr = matr.mats[[digit]].MatrixPower[reflect, funcs[[digit, 1]][i + nonrep + special]] ]; vect + mat.Inverse[IdentityMatrix[3] - matr].vectr, vect ] ] Block[{order = 2, max = 85/126, funcs, data, str, cases, i}, funcs = Table[{0 &, 0, 1}, {i, 1, 9, 1}]; funcs[[1]] = {Piecewise[{{0, # == 0}, {1, # < order - 1}, {0, # == order - 1}}, 0] &, 2}; funcs[[2]] = {Piecewise[{{1, # < order - 1}, {0, # == order - 1}}, 0] &, order, 2}; funcs[[3]] = {Piecewise[{{1, # < order - 1}, {0, # == order - 1}}, 0] &, order, 2}; funcs[[4]] = {Piecewise[{{0, # < order - 1}, {1, # == order - 1}}, 0] &, order, 2}; funcs[[5]] = {Piecewise[{{0, # < order - 1}, {1, # == order - 1}}, 0] &, order, 2}; funcs[[6]] = {Piecewise[{{1, # < order - 1}, {0, # == order - 1}}, 0] &, order, 2}; funcs[[7]] = {Piecewise[{{1, # < order - 1}, {0, # == order - 1}}, 0] &, order, 2}; funcs[[8]] = {Piecewise[{{1, # < order - 1}, {0, # == order - 1}}, 0] &, order, 2}; data = Table[ Trilbert[ Block[{num = (i + max/8)/(8^order)}, If[And[0 <= num, num <= 1], num, Mod[num, 1]]], funcs], {i, 0, 8^order - 1}]; str = StringJoin @@ (Differences[data]*(2^order) /. {{1, 0, 0} -> "x", {0, 1, 0} -> "y", {0, 0, 1} -> "z", {-1, 0, 0} -> "a", {0, -1, 0} -> "b", {0, 0, -1} -> "c"}); cases = DeleteDuplicates[StringCases[str, a__ ~~ a_]]; Print[ If[Length[cases] < 1, cases, Table[{cases[[i]], StringPosition[str, cases[[i]]]}, {i, 1, Length[cases], 1}] // MatrixForm ] ]; str = StringJoin @@ (Abs[ Differences[data]*(2^order)] /. {{1, 0, 0} -> "x", {0, 1, 0} -> "y", {0, 0, 1} -> "z"}); cases = DeleteDuplicates[StringCases[str, a__ ~~ a_]]; Print[ If[Length[cases] < 1, cases, Table[{cases[[i]], StringPosition[str, cases[[i]]]}, {i, 1, Length[cases], 1}] // MatrixForm ] ]; Graphics3D[Tube[data, 1/50], Axes -> True, AxesLabel -> {"x", "y", "z"}, PlotRange -> {{-1/20, 21/20}, {-1/20, 21/20}, {-1/20, 21/20}}] ] That same notebook contains an isomorphic pair of statements defining and invoking "Spirilbert", which sweeps the cube in a sort of spiral, instead of by octants. --rwg