oftenpaper.net

the sierpinski triangle page to end most sierpinski triangle pages ™


Constructing the Sierpinski triangle

Throughout my years playing around with fractals, the Sierpinski triangle has been a consistent staple. The triangle is named after Wacław Sierpiński and as fractals are wont the pattern appears in many places, so there are many different ways of constructing the triangle on a computer.

All of the methods are fundamentally iterative. The most obvious method is probably the triangle-in-triangle approach. We start with one triangle, and at every step we replace each triangle with 3 subtriangles:

  1. axiom = Polygon[{{0, 0}, {1, Sqrt[3]}/2, {1, 0}}];
    
    next[prev_] := prev /. Polygon[{p1_, p2_, p3_}] :> {
         Polygon[{p1, (p1 + p2)/2, (p1 + p3)/2}],
         Polygon[{p2, (p2 + p3)/2, (p1 + p2)/2}],
         Polygon[{p3, (p1 + p3)/2, (p2 + p3)/2}]};
    
    draw[n_] := Graphics[{EdgeForm[Black],
        Nest[next, N@axiom, n]}];
    

This triangle-in-triangle method strikes me as a disguised Lindenmayer system. L-systems are iterative symbol-based replacement mechanisms. There are a variety of more explicit L-system constructions for the triangle, such as the 'arrowhead' L-system (also see my L-systems program):

  1. axiom = {A};
    rules = {A -> {B, R, A, R, B}, B -> {A, L, B, L, A}};
    conversions = {A -> forward, B -> forward, L -> left, R -> right};
    
    (* state transformations *)
    forward[{z_, a_}] := {z + E^(I a), a};
    left[{z_, a_}] := {z, a + 2 Pi/6};
    right[{z_, a_}] := {z, a - 2 Pi/6};
    
    draw[n_] := Module[{program, zs},
      program = Flatten[Nest[# /. rules &, axiom, n]] /. conversions;
      zs = First /@ ComposeList[program, N@{0, 0}];
      Graphics[Line[{Re[#], Im[#]} & /@ First /@ Split[zs]]]];
    

There's the cellular automata approach, where the 'world' is a single array of bits and at each "instant" we alter a bit based on the state of it and its neighbors. If we plot the evolution of Rule 22 (and others), we get these patterns:

  1. draw[n_] := ArrayPlot[CellularAutomaton[22, {{1}, 0}, n]];
    

There are bound to be many elementary number-theoretic constructions of the Sierpinski triangle given that it looks like a percolation pattern (as in the cellular automata above). The Wikipedia article mentions that it appears in Pascal's Triangle when differentiating between even and odd numbers. Sure enough:

  1. draw[n_] := Module[{t},
        t = Table[Binomial[m, k], {m, 0, n}, {k, 0, m}];
    
        Column[Row[#, " "] & /@ t, Center] /. {
            x_?EvenQ :> Style[Framed[x], LightGray],
            x_?OddQ :> Framed[x]}];
    

If we look at these Pascal forms and reverse engineer the parity rules, we get Rule 22. Though it might depend on what exactly you're reverse engineering. We can generalize from even/odd to other moduli:

  1. Pascal's triangle mod 4

  2. Pascal's triangle $x\equiv 2\space(\text{mod }4)$

  3. draw[n_] := Module[{t},
       t = Table[Mod[Binomial[m, k], 4], {m, 0, n}, {k, 0, m}];
    
       Column[Row[#, " "] & /@ t, Center] /. x_?NumberQ :>
         Style[Framed["  ", FrameStyle -> None],
          Background -> ColorData[3][2 + x]]];
    

The Wikipedia article for Pascal's triangle mentions that we can construct a 'Pascal matrix' using the matrix exponential:

$$ e^{\left( \begin{array}{ccccc} 0 & 0 & 0 & 0 & 0 \\ 1 & 0 & 0 & 0 & 0 \\ 0 & 2 & 0 & 0 & 0 \\ 0 & 0 & 3 & 0 & 0 \\ 0 & 0 & 0 & \ddots & 0 \\ \end{array} \right)}=\left( \begin{array}{ccccc} 1 & 0 & 0 & 0 & 0 \\ 1 & 1 & 0 & 0 & 0 \\ 1 & 2 & 1 & 0 & 0 \\ 1 & 3 & 3 & 1 & 0 \\ 1 & 4 & 6 & 4 & \ddots \\ \end{array} \right) $$

"Ah, that makes sense." You say. Indeed, but what's cool is that we then have a pedantic way of specifying the Sierpinski triangle:

$$ \mathfrak{S}\equiv e^{\left( \begin{array}{ccccc} 0 & 0 & 0 & 0 & 0 \\ 1 & 0 & 0 & 0 & 0 \\ 0 & 2 & 0 & 0 & 0 \\ 0 & 0 & 3 & 0 & 0 \\ 0 & 0 & 0 & \ddots & 0 \\ \end{array} \right)}(\text{mod } 2) $$

This equation is in what's called "straight ballin'" form, and it gives us a fancy way of producing the triangle:

  1. draw[n_] := ArrayPlot[Mod[MatrixExp[DiagonalMatrix[Range[n], -1]], 2]];
    

Heawt deaowg /drawl. It's not very performant though. The following is faster and arguably more elegant:

  1. draw[n_] := ArrayPlot[Mod[Array[Binomial, {n, n}, 0], 2]];
    

Along these lines, it shouldn't be surprising that the Sierpinski pattern appears in other combinatorial expressions, such as the Stirling numbers:

  1. draw[n_] := Grid[Partition[#, 2]] &@
       Table[ArrayPlot[Mod[Array[f, {n, n}], 2],
         PlotLabel -> f, FrameStyle -> LightGray],
        {f, {Binomial, StirlingS1, StirlingS2, Multinomial}}];
    

If we treat the rows produced by these combinatorial functions as arrays of bits, what sequence of numbers do the bits represent? There's a variety of ways to interpret this question, but here's one assortment:

  1. $$ \left( \begin{array}{ccccccccccc} \text{Binomial} & 1 & 3 & 5 & 15 & 17 & 51 & 85 & 255 & 257 & \ldots \\ \text{StirlingS1} & 1 & 1 & 3 & 3 & 5 & 5 & 15 & 15 & 17 & \ldots \\ \text{StirlingS2} & 1 & 1 & 3 & 7 & 13 & 29 & 55 & 115 & 209 & \ldots \\ \text{Multinomial} & 511 & 341 & 409 & 273 & 481 & 321 & 385 & 257 & 255 & \ldots \\ \end{array} \right) $$
  2. 
    draw[n_] := With[{dropZeros = # /. {x__, 0 ..} :> {x} &},
       MatrixForm[Table[Flatten[
          {f, FromDigits[dropZeros[#], 2] & /@ Mod[Array[f, {n, n}, 0], 2], "\[Ellipsis]"}],
         {f, {Binomial, StirlingS1, StirlingS2, Multinomial}}]]];
    
    

The first, second, and fourth sequences are versions of each other, tautologically described in OEIS as A001317. The sequence for the Stirling numbers of the second kind doesn't seem to have any fame, but if you shift its bits around you can find A099901 and A099902.

The Wikipedia article for the Sierpinski triangle mentions its appearance in logic tables such as this one. If you stare blankly at that image long enough you'll notice it's a set-inclusion table. Take the subsets of a set and pair them against each other under set-inclusion (is subset A a subset of subset B?) and you will get that table.

Personally that's a more interesting interpretation than the binary logic one, though the apparent distinction between these subjects is likely just a matter of perspective. Another set-related Sierpinski pattern I found is set disjunction (when sets have no common elements):

  1. isSubset[a_, b_] := Union[a, b] == b;
    areDisjoint[a_, b_] := Intersection[a, b] == {};
    
    subs[0] = {{}};
    subs[n_] := Module[{s = subs[n - 1]},
       Join[s, Append[#, n] & /@ s]];
    
    draw[n_] := Grid[List[Table[
         ArrayPlot[Boole[Outer[f, subs[n], subs[n], 1]],
          PlotLabel -> f, FrameStyle -> LightGray],
         {f, {isSubset, areDisjoint}}]]]
    

One thing I noticed is that these set patterns depend on the order in which you place the subsets. It has to be the same order that you would get if you were constructing the subsets iteratively. I also wasn't able to find a straightforward ranking function that would order the sets into this iterative sequence. Mathematica's Combinatorica package refers to it as the binary ordering. I think I'm starting to understand what Gandalf meant when he said

" The Sierpinski triangle cannot-be wrought without heed to the creeping tendrils of recursion. Even the binomial coefficient has factorials which are recursively defined. "

MathWorld mentions a broader context for why binary logic can be used in the construction of the Sierpinski triangle. Namely the Lucas correspondence theorem which states that given two numbers written in a prime base,

$$n=n_mp^m+\cdots+n_1p^1+n_0p^0\space\space\space(0\le n_i\le p)$$ $$k=k_mp^m+\cdots+k_1p^1+k_0p^0\space\space\space(0\le k_i\le p)$$

We can get their binomial coefficient modulo that prime by performing binomial coefficients digit-wise and multiplying the results.

$$\binom{n}{k}=\prod _{i=0}^m \binom{n_i}{k_i}(\text{mod }p)$$

The binomial coefficient $\binom{n}{k}$ represents the number of $k$-element subsets of a set of $n$ elements. If we're using zeros and ones, then:

  1. $$ \begin{array}{cc} \binom{0}{0}=1 & \binom{0}{1}=0 \\ \binom{1}{0}=1 & \binom{1}{1}=1 \\ \end{array} $$
  2. 
    TraditionalForm[Grid[Outer[
        HoldForm[Binomial[##]] == Binomial[##] &,
      {0, 1}, {0, 1}]]]
    
    

The factorial definition is interesting in this case.

$$\binom{n}{k} = \frac {n!} {k!(n-k)!}$$

Notice that if we have $\binom 0 1$, we get the factorial of a negative number in the denominator. By sticking with the recursive definition of the factorial, the conclusion is that the denominator is some flavor of $\infty$, so you have $\frac 1 \infty=0$. ($0!$ is defined as 1).

The binary operation I found in our little binary binomial table was NOTing $n$, ANDing the result with $k$, and then NOTing that: $\neg(\neg n\land k)=n\lor \neg k$. Also notice it's equivalent to the greater than or equal to operation $n \ge k$.

If by some stroke of luck we happen to have the two numbers stored in binary on our computer, these operations can be performed atomically on the numbers as a whole. And since we're multiplying everything at the end, any presence of $\binom 0 1$ in the original numbers means the binomial is congruent to 2. The only trick would be tracking whatever the most significant bit of either number was.

  1. binaryBinomial[a_, b_] := Module[{bits},
       bits = IntegerDigits[{a, b}, 2];
       bits = PadLeft[#, Max[Length /@ bits]] & /@ bits;
    
       Boole[FreeQ[Transpose[bits], {0, 1}]]];
    
    draw[n_] := MatrixPlot[
       Array[binaryBinomial, {2^n, 2^n}, 0],
       Frame -> None];
    

There's a lot of related patterns:

  1. binaryWhoKnows[a_, b_] :=
      DigitCount[BitOr[a, BitNot[b]], 3, 1];
    
    draw[n_] := MatrixPlot[
       Array[binaryWhoKnows, {2^n, 2^n}, 0],
       Frame -> False];
    

And look what I found!

$$ 2 b\lor \neg 2 b = \text{true} $$

If we're looking for a one- or two-liner that's one- or two-linear in languages beside Mathematica, we'd have trouble doing better than the chaos game algorithm, which goes like this:

1 start at any point. call it p
2 pick one of the three vertices at random
3 find the point halfway between p and that vertex
4 call that point p and draw it
5 goto 2
  1. vertices = {{0, 0}, {1, Sqrt[3]}/2, {1, 0}};
    
    draw[numPoints_] := Graphics[{
        PointSize[0], Opacity[.1],
        Point[FoldList[(#1 + #2)/2 &, N@{0, 0},
          RandomChoice[N@vertices, numPoints]]]}];
    
  2. vertices = {{0, 0}, {1, Sqrt[3]}/2, {1, 0}};
    
    draw[numPoints_] := Graphics[{
        PointSize[0], Opacity[.01],
        Point[FoldList[(#1 + #2)/2 &, N@{0, 0},
          RandomChoice[N@vertices, numPoints]]]},
       ImageSize -> 2 1280];
    
    draw[50000000] // ImageAdjust // ImageResize[#, Scaled[1/2]] &
    

The chaos game doesn't render as crisply as a lot of the other methods, especially without transparency effects, but it has the advantage of being highly performant. It runs about one million points per second on my laptop. Mind you this is with Mathematica's RNG, which is not your everyday math.rand().

One thing I realized is that the randomness isn't actually a necessary aspect of the general algorithm. It's used as an approximating force (or perhaps something a bit more subtle than that). Otherwise with enough spacetime on your computer you can just perform all possible half-distancings:

  1. axiom = Polygon[{{0, 0}, {1, Sqrt[3]}/2, {1, 0}}];
    
    next[prev_] := prev /. Polygon[{p1_, p2_, p3_}] :> { 
         Polygon[{p1, (p1 + p2)/2, (p1 + p3)/2}],
         Polygon[{p2, (p2 + p3)/2, (p1 + p2)/2}],
         Polygon[{p3, (p1 + p3)/2, (p2 + p3)/2}]};
    
    draw[n_] := Graphics[{PointSize[0], Opacity[.1],
        Nest[next, N@axiom, n] /. Polygon -> Point}];
    
    
    
  2. axiom = Polygon[{{0, 0}, {1, Sqrt[3]}/2, {1, 0}}];
    
    next[prev_] := prev /. Polygon[{p1_, p2_, p3_}] :> {
         Polygon[{p1, (p1 + p2)/2, (p1 + p3)/2}],
         Polygon[{p2, (p2 + p3)/2, (p1 + p2)/2}],
         Polygon[{p3, (p1 + p3)/2, (p2 + p3)/2}]};
    
    points[n_] := DeleteDuplicates[Flatten[
        Nest[next, axiom, n] /. Polygon -> Sequence, n]];
    
    points[5]
    

These images look basically the same. Not surprising since they're both point-based. But I gander the distinction between these two algorithms may have been more than just an issue of curiousity 20 years ago. I still remember my first computer, the alien-processored TI-85, chugging away furiously for a good half a minute before the triangle became clear.

Notice that this specific algorithm is actually just a minor modification of the triangle-in-triangle algorithm. The difference is that polygon vertices are here rendered as points. This modification is possible because of Mathematica's symbolic semantics. The symbol Polygon is meaningless until it's processed by the Graphics function. Until then, we can perform structural operations such as replacing it by the Point symbol. In fact the following is completely valid:

axiom = triangle[{{0, 0}, {1, Sqrt[3]}/2, {1, 0}}];

next[prev_] := prev /. triangle[{p1_, p2_, p3_}] :> {
     triangle[{p1, (p1 + p2)/2, (p1 + p3)/2}],
     triangle[{p2, (p2 + p3)/2, (p1 + p2)/2}],
     triangle[{p3, (p1 + p3)/2, (p2 + p3)/2}]};

draw[n_] := Graphics[Nest[next, N@axiom, n] /.  triangle :> Polygon ];

triangle here doesn't have any meaning, ever, until we replace it:

  1.  triangle :> Polygon

  2.  triangle :> Line

  3. triangle[pts_] :> Line[RandomChoice[pts, RandomInteger[{2, 3}]]]

  4. triangle[pts_] :> Disk[Mean[pts], 1/2^(n + 1)]

  5. triangle[pts_] :> Sphere[Append[Mean[pts], 0], 1/2^(n + 1)]
    

Sidenote. What do you get when you methodically build a Lisp on top of symbolic replacement semantics? You get the Mathematica language, of which Mathematica and Mathics appear to be the only incarnations.

Let's say you forgot how to multiply matrices. Well, just type in some symbols and see the results empirically:


{{a, b}, {c, d}} . {{e, f}, {g, h}} // MatrixForm
$$ \left( \begin{array}{cc} a e+b g & a f+b h \\ c e+d g & c f+d h \\ \end{array} \right) $$

If that's still confusing, you can use strings, colored text, graphics, images, etc. instead of symbols. In fact if you have a Tron zapper you can even zap your cat into Mathematica and have him fill up one of those matrix slots, for the advancement of science.

  1. kitty = WolframAlpha["cat picture", "PodImages"][[2]];
    
    (* see http://mathematica.stackexchange.com/a/8291/950 *)
    text = First[First[ImportString[ExportString[
          Style["IM IN UR MATRIX...", FontFamily -> "Impact"], "PDF"]]]];
    
    sym = Framed[Overlay[{kitty,
         Graphics[{EdgeForm[Black], White, text}, ImageSize -> 150,
         PlotRangePadding -> 0]}], FrameStyle -> LightGray];
    
    {{a, b}, {Magnify[sym, 1/2], d}} . {{e, f}, {g, h}} // MatrixForm
    

There's poor Mr. Scruples. Our neighbor will miss him.

The exponential identity for the Pascal matrix is not difficult to understand based on the series definition of the exponential function:

$$ e^x=\frac{x^0}{0!}+\frac{x^1}{1!}+\frac{x^2}{2!}+\frac{x^3}{3!}+\frac{x^4}{4!}+\frac{x^5}{5!}+\cdots $$

You could work out the matrix arithmetic by hand, or you could do this:


power[n_, p_] := MatrixPower[
    DiagonalMatrix[ToString /@ Range[n], -1], p] // MatrixForm;

Grid[Partition[Table[power[6, p], {p, 1, 6}], 3]] /. 0 -> "\[CenterDot]"

$$ \begin{array}{ccc} \left( \begin{array}{ccccccc} \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ 1 & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & 2 & \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & 3 & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & 4 & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & 5 & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot & 6 & \cdot \\ \end{array} \right) & \left( \begin{array}{ccccccc} \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ 1 2 & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & 2 3 & \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & 3 4 & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & 4 5 & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & 5 6 & \cdot & \cdot \\ \end{array} \right) & \left( \begin{array}{ccccccc} \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ 1 2 3 & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & 2 3 4 & \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & 3 4 5 & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & 4 5 6 & \cdot & \cdot & \cdot \\ \end{array} \right) \\ \left( \begin{array}{ccccccc} \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ 1 2 3 4 & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & 2 3 4 5 & \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & 3 4 5 6 & \cdot & \cdot & \cdot & \cdot \\ \end{array} \right) & \left( \begin{array}{ccccccc} \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ 1 2 3 4 5 & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & 2 3 4 5 6 & \cdot & \cdot & \cdot & \cdot & \cdot \\ \end{array} \right) & \left( \begin{array}{ccccccc} \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ 1 2 3 4 5 6 & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ \end{array} \right) \\ \end{array} $$

These are the first 6 powers of the subdiagonal matrix. You can see that the diagonal gets multiplied by subsequently shifted versions of itself, so the calculation ends up creating factorial products. For example, 3x4x5x6 (in the fourth power) can be written in terms of factorials as $6!/2!$. If we factor in the denominator from the series for $e$, we have

$$\frac {6!} {4!2!}$$

From the factorial definition of the binomial coefficient:

$$\binom{n}{k} = \frac {n!} {k!(n-k)!}$$

We see that this particular slot in the matrix is $\binom 6 4$. The binomial coefficient itself is of course directly related to Pascal's triangle. Also notice that every power of the matrix has its numbers on a different diagonal, so when we sum up all the powers there is no interaction to account for. Every term in the series is a distinct diagonal of Pascal's triangle.

Powers of matrices have a well-known interpretation in terms of graph walks/probabilities. I didn't find anything interesting along this line though. What about graphs represented by the Sierpinski matrix itself?


$$\left( \begin{array}{cccccccccccccccc} 1 & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ 1 & 1 & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ 1 & \cdot & 1 & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ 1 & 1 & 1 & 1 & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ 1 & \cdot & \cdot & \cdot & 1 & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ 1 & 1 & \cdot & \cdot & 1 & 1 & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ 1 & \cdot & 1 & \cdot & 1 & \cdot & 1 & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ 1 & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & 1 & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ 1 & 1 & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & 1 & 1 & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ 1 & \cdot & 1 & \cdot & \cdot & \cdot & \cdot & \cdot & 1 & \cdot & 1 & \cdot & \cdot & \cdot & \cdot & \cdot \\ 1 & 1 & 1 & 1 & \cdot & \cdot & \cdot & \cdot & 1 & 1 & 1 & 1 & \cdot & \cdot & \cdot & \cdot \\ 1 & \cdot & \cdot & \cdot & 1 & \cdot & \cdot & \cdot & 1 & \cdot & \cdot & \cdot & 1 & \cdot & \cdot & \cdot \\ 1 & 1 & \cdot & \cdot & 1 & 1 & \cdot & \cdot & 1 & 1 & \cdot & \cdot & 1 & 1 & \cdot & \cdot \\ 1 & \cdot & 1 & \cdot & 1 & \cdot & 1 & \cdot & 1 & \cdot & 1 & \cdot & 1 & \cdot & 1 & \cdot \\ 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 \\ \end{array} \right)$$

Those were more interesting:

  1. graph[n_] := GraphPlot3D[
       Mod[Array[Binomial, {n, n}, 0], 2],
       Method -> "HighDimensionalEmbedding", Boxed -> False,
       VertexRenderingFunction -> ({White, Sphere[#, .05]} &),
       PlotStyle -> {Thick, Hue[2/3, 2/3, 2/3]}];
    

Note this is a 3D graph layout. It has some pretty symmetries. I did some tiresome work trying to figure out what polyhedron it might be.

Tooltip[PolyhedronData[#], #] & /@ Select[
  PolyhedronData[], PolyhedronData[#, "VertexCount"] == 14 &]

After much time, I find. It's the tetrakis hexahedron:

  1. Graphics3D[{Opacity[.94], FaceForm[Gray],
      PolyhedronData["TetrakisHexahedron", "Faces"]},
     Lighting -> "Neutral", Boxed -> False]
    

I'm certain it's this particular figure because we can just build a graph from its vertex data and then do a graph isomorphism check. And look, we can run this polyhedron grapherizer willy-nilly allabouts, like on the Archimedean solids:

  1. polyGraph[poly_, options___] :=
      Graph[UndirectedEdge @@@ PolyhedronData[poly, "EdgeIndices"], options];
    
    sierpinskiMatrixGraph[n_, options___] := Module[{a},
       (*symmetrize and remove self-loops to allow general isomorphism
       comparison. note that we remove the "inner" vertex since
       we're comparing the "external" geometry as rendered*)
       a = Mod[Array[Binomial, {n, n}, 1], 2];
       AdjacencyGraph[a + Transpose[a] /. 2 -> 0, options]];
    
    IsomorphicGraphQ[polyGraph["TetrakisHexahedron"], sierpinskiMatrixGraph[14]]
    IsomorphicGraphQ[polyGraph["CumulatedCube"], sierpinskiMatrixGraph[14]]
    
    Grid[Partition[#, 4]] &[
      polyGraph[#, PlotLabel -> PolyhedronData[#, "Name"]] & /@
       PolyhedronData["Archimedean"]]
    

Here are the first few powers of the Sierpinski matrix:

  1. $$ \begin{array}{cc} \left( \begin{array}{cccccccccccccccc} 1 &\cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ 1 & 1 & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ 1 & \cdot & 1 & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ 1 & 1 & 1 & 1 & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ 1 & \cdot & \cdot & \cdot & 1 & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ 1 & 1 & \cdot & \cdot & 1 & 1 & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ 1 & \cdot & 1 & \cdot & 1 & \cdot & 1 & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ 1 & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & 1 & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ 1 & 1 & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & 1 & 1 & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ 1 & \cdot & 1 & \cdot & \cdot & \cdot & \cdot & \cdot & 1 & \cdot & 1 & \cdot & \cdot & \cdot & \cdot & \cdot \\ 1 & 1 & 1 & 1 & \cdot & \cdot & \cdot & \cdot & 1 & 1 & 1 & 1 & \cdot & \cdot & \cdot & \cdot \\ 1 & \cdot & \cdot & \cdot & 1 & \cdot & \cdot & \cdot & 1 & \cdot & \cdot & \cdot & 1 & \cdot & \cdot & \cdot \\ 1 & 1 & \cdot & \cdot & 1 & 1 & \cdot & \cdot & 1 & 1 & \cdot & \cdot & 1 & 1 & \cdot & \cdot \\ 1 & \cdot & 1 & \cdot & 1 & \cdot & 1 & \cdot & 1 & \cdot & 1 & \cdot & 1 & \cdot & 1 & \cdot \\ 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 \\ \end{array} \right) & \left( \begin{array}{cccccccccccccccc} 1 & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ 2 & 1 & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ 2 & \cdot & 1 & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ 4 & 2 & 2 & 1 & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ 2 & \cdot & \cdot & \cdot & 1 & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ 4 & 2 & \cdot & \cdot & 2 & 1 & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ 4 & \cdot & 2 & \cdot & 2 & \cdot & 1 & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ 8 & 4 & 4 & 2 & 4 & 2 & 2 & 1 & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ 2 & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & 1 & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ 4 & 2 & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & 2 & 1 & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ 4 & \cdot & 2 & \cdot & \cdot & \cdot & \cdot & \cdot & 2 & \cdot & 1 & \cdot & \cdot & \cdot & \cdot & \cdot \\ 8 & 4 & 4 & 2 & \cdot & \cdot & \cdot & \cdot & 4 & 2 & 2 & 1 & \cdot & \cdot & \cdot & \cdot \\ 4 & \cdot & \cdot & \cdot & 2 & \cdot & \cdot & \cdot & 2 & \cdot & \cdot & \cdot & 1 & \cdot & \cdot & \cdot \\ 8 & 4 & \cdot & \cdot & 4 & 2 & \cdot & \cdot & 4 & 2 & \cdot & \cdot & 2 & 1 & \cdot & \cdot \\ 8 & \cdot & 4 & \cdot & 4 & \cdot & 2 & \cdot & 4 & \cdot & 2 & \cdot & 2 & \cdot & 1 & \cdot \\ 16 & 8 & 8 & 4 & 8 & 4 & 4 & 2 & 8 & 4 & 4 & 2 & 4 & 2 & 2 & 1 \\ \end{array} \right) \\ \left( \begin{array}{cccccccccccccccc} 1 & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ 3 & 1 & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ 3 & \cdot & 1 & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ 9 & 3 & 3 & 1 & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ 3 & \cdot & \cdot & \cdot & 1 & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ 9 & 3 & \cdot & \cdot & 3 & 1 & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ 9 & \cdot & 3 & \cdot & 3 & \cdot & 1 & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ 27 & 9 & 9 & 3 & 9 & 3 & 3 & 1 & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ 3 & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & 1 & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ 9 & 3 & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & 3 & 1 & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ 9 & \cdot & 3 & \cdot & \cdot & \cdot & \cdot & \cdot & 3 & \cdot & 1 & \cdot & \cdot & \cdot & \cdot & \cdot \\ 27 & 9 & 9 & 3 & \cdot & \cdot & \cdot & \cdot & 9 & 3 & 3 & 1 & \cdot & \cdot & \cdot & \cdot \\ 9 & \cdot & \cdot & \cdot & 3 & \cdot & \cdot & \cdot & 3 & \cdot & \cdot & \cdot & 1 & \cdot & \cdot & \cdot \\ 27 & 9 & \cdot & \cdot & 9 & 3 & \cdot & \cdot & 9 & 3 & \cdot & \cdot & 3 & 1 & \cdot & \cdot \\ 27 & \cdot & 9 & \cdot & 9 & \cdot & 3 & \cdot & 9 & \cdot & 3 & \cdot & 3 & \cdot & 1 & \cdot \\ 81 & 27 & 27 & 9 & 27 & 9 & 9 & 3 & 27 & 9 & 9 & 3 & 9 & 3 & 3 & 1 \\ \end{array} \right) & \left( \begin{array}{cccccccccccccccc} 1 & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ 4 & 1 & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ 4 & \cdot & 1 & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ 16 & 4 & 4 & 1 & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ 4 & \cdot & \cdot & \cdot & 1 & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ 16 & 4 & \cdot & \cdot & 4 & 1 & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ 16 & \cdot & 4 & \cdot & 4 & \cdot & 1 & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ 64 & 16 & 16 & 4 & 16 & 4 & 4 & 1 & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ 4 & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & 1 & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ 16 & 4 & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & 4 & 1 & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ 16 & \cdot & 4 & \cdot & \cdot & \cdot & \cdot & \cdot & 4 & \cdot & 1 & \cdot & \cdot & \cdot & \cdot & \cdot \\ 64 & 16 & 16 & 4 & \cdot & \cdot & \cdot & \cdot & 16 & 4 & 4 & 1 & \cdot & \cdot & \cdot & \cdot \\ 16 & \cdot & \cdot & \cdot & 4 & \cdot & \cdot & \cdot & 4 & \cdot & \cdot & \cdot & 1 & \cdot & \cdot & \cdot \\ 64 & 16 & \cdot & \cdot & 16 & 4 & \cdot & \cdot & 16 & 4 & \cdot & \cdot & 4 & 1 & \cdot & \cdot \\ 64 & \cdot & 16 & \cdot & 16 & \cdot & 4 & \cdot & 16 & \cdot & 4 & \cdot & 4 & \cdot & 1 & \cdot \\ 256 & 64 & 64 & 16 & 64 & 16 & 16 & 4 & 64 & 16 & 16 & 4 & 16 & 4 & 4 & 1 \\ \end{array} \right) \\ \end{array} $$
  2. power[n_, p_] := MatrixPower[
       Mod[Array[Binomial, {n, n}, 0], 2], p];
    
    Grid[Partition[#, 2]] &@
     Table[
      MatrixForm[power[16, p] /. 0 -> "\[CenterDot]"],
      {p, 1, 4}]
    

There's a lot of patterns here. For one, the powers of the Sierpinski matrix are Sierpinski matrices! This isn't necessarily interesting though. The powers of a triangular matrix are going to be triangular. But the numbers follow a curious sequence of powers. For example, in the third power we have the sequence {1, 3, 3, 9, 3, 9, 9, 27, 3, ... }. And this sequence occurs in every column and every row of the matrix, if you hop over the zeros. We can normalize the powers to find:

  1. $$ \left( \begin{array}{cccccccccccccccc} 0 & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ 1 & 0 & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ 1 & \cdot & 0 & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ 2 & 1 & 1 & 0 & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ 1 & \cdot & \cdot & \cdot & 0 & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ 2 & 1 & \cdot & \cdot & 1 & 0 & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ 2 & \cdot & 1 & \cdot & 1 & \cdot & 0 & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ 3 & 2 & 2 & 1 & 2 & 1 & 1 & 0 & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ 1 & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & 0 & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ 2 & 1 & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & 1 & 0 & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ 2 & \cdot & 1 & \cdot & \cdot & \cdot & \cdot & \cdot & 1 & \cdot & 0 & \cdot & \cdot & \cdot & \cdot & \cdot \\ 3 & 2 & 2 & 1 & \cdot & \cdot & \cdot & \cdot & 2 & 1 & 1 & 0 & \cdot & \cdot & \cdot & \cdot \\ 2 & \cdot & \cdot & \cdot & 1 & \cdot & \cdot & \cdot & 1 & \cdot & \cdot & \cdot & 0 & \cdot & \cdot & \cdot \\ 3 & 2 & \cdot & \cdot & 2 & 1 & \cdot & \cdot & 2 & 1 & \cdot & \cdot & 1 & 0 & \cdot & \cdot \\ 3 & \cdot & 2 & \cdot & 2 & \cdot & 1 & \cdot & 2 & \cdot & 1 & \cdot & 1 & \cdot & 0 & \cdot \\ 4 & 3 & 3 & 2 & 3 & 2 & 2 & 1 & 3 & 2 & 2 & 1 & 2 & 1 & 1 & 0 \\ \end{array} \right) $$
  2. 
    
    
    power[n_, p_] := MatrixPower[
       Mod[Array[Binomial, {n, n}, 0], 2], p];
    
    MatrixForm /@ Table[
      IntegerExponent[power[16, p], p] /. Infinity -> "\[CenterDot]",
      {p, 2, 4}]
    
    
    
    

This is the sequence in terms of the exponent, and it applies to each power of the Sierpinski matrix, including the first power. For example, 3 to the power of each of {0, 1, 1, 2, 1, 2, 2, 3, 1, ...} is {1, 3, 3, 9, 3, 9, 9, 27, 3, ...}. This power sequence appears in OEIS as the number of ones in the binary representation of n, among other descriptions.

Here is a totally practical application of all of this. A pretty array of buttons:

  1. power[n_, p_] := MatrixPower[Transpose@Reverse@
         Mod[Array[Binomial, {n, n}, 0], 2], p];
    
    Grid[Partition[#, 4]] &@
     With[{m = "you, are now infused, with, the power of, dot, dot, dot... "},
      Array[Function[p,
        Button[Rotate[#, -Pi/4], Speak[m <> ToString[p]]] &@
         Rasterize@MatrixPlot[IntegerExponent[power[2^p, 4], 10],
           ImageSize -> 94, Frame -> None, PlotRangePadding -> 0]],
       8]] 
    

The Towers of Hanoi is a variation on the sticks-in-holes game where instead of putting sticks in holes, you put holes around sticks. Thus the game is ultimately a quaint philosophical remark on the roles of the sexes. But for our purposes there is a claim on the internets that the states of the game form Sierpinski triangle-like graphs:

  1. validQ[s_state] := And @@ Less @@@ s;
    
    (*do all physically possible moves. remove invalid moves afterward.*)
    neighbors[states : {__state}] := Select[#, validQ] &@
       DeleteDuplicates@Flatten@
         Table[Module[{st2 = st},
    
           If[Length[st2[[from]]] > 0,
            PrependTo[st2[[to]], st2[[from, 1]]];
            st2[[from]] = Rest[st2[[from]]]];
    
           If[st2 =!= st && validQ[st2],
            Sow@UndirectedEdge[st, st2]];
    
           st2], {st, states}, {to, Length[st]}, {from, Length[st]}];
    
    toStyle[expr_] := expr /. s_state :> (
         Property[Tooltip[s, MatrixForm /@ List @@ s],
          VertexStyle -> {EdgeForm[None],
            ColorData[3][1 + Length[s] - Count[s, {}]]}]);
    
    hanoiGraph[s_, options___] := Module[{vertices, edges, n},
       n = Count[s, _Integer, Infinity];
       {vertices, {edges}} = Reap[Nest[neighbors, {s}, 2^n]];
    
       SetAttributes[UndirectedEdge, Orderless];
       Graph[toStyle[vertices], DeleteDuplicates[edges],
        options(*,GraphLayout->"SpringEmbedding"*)(*,
        VertexShapeFunction->(Style[#,7,Black]&@
        Text[Row[MatrixForm/@List@@#2],#1]&)*)]];
    
    hanoiGraph[state[{}, {}, Range[4]],
     Epilog -> Inset[Rotate[Style["F-", 300, Bold, Red, Opacity[.65]], Pi/7]]]
    

Which, as you can see, is a lie if I've ever seen one (internets, you are now on notice). Then again, if you fiddle with the layout and you squint a bit, you can kinda see it, but it's the sort of Sierpinski triangle that Maddox would stamp a huge red F over. To be clear, each vertex represents a single state of the game, and vertices are connected if there is a legal move between those states.

The nice thing about this algorithm is that at each step it just blindly constructs all possibilities, which is easy, and then afterwards removes the ones that aren't valid, which is also easy. Point being it works in broad strokes. And at the end of it we have a map to follow if we ever get stuck. You can do this sort of thing for all sorts of things, like say Rubik's cube. Though I don't know if the combinatorics are favorable in its case. The Towers of Hanoi can be played with more than three sticks:

  1. hanoiGraph[state[{}, {}, {}, Range[4]]]

  2. hanoiGraph[state[{}, {}, {}, {}, Range[4]]]
    

  3. hanoiGraph[state[{}, {}, Range[3], Range[3]]]
    

  4. hanoiGraph[state[{}, {}, {1}, Range[3]]]
    

  5. hanoiGraph[state[{}, {}, {2}, Range[3]]]
    

  6. hanoiGraph[state[{}, {}, {3}, Range[3]]]
    

  7. validQ[s_state] := And @@ Equal @@@ s;
    hanoiGraph[state[{}, {}, ConstantArray[1, 5]]]
    

  8. validQ[s_state] := And @@ LessEqual @@@ s;
    hanoiGraph[state[{}, ConstantArray[2, 3], ConstantArray[1, 3]]]

  9. validQ[s_state] := And @@ Equal @@@ s;
    hanoiGraph[state[{}, ConstantArray[2, 3], ConstantArray[1, 3]]]

"WHAT THE HELL IS THAT", you say. Indeed, it's messy because it's a low-D rendering. We can also play variations of the game that allow multiple holes of the same diameter, or variations where we adjust the rules a bit. In higher dimensions you can see the structure better:

    1. validQ[s_state] := And @@ Less @@@ s;
      
      (*do all physically possible moves. remove invalid moves afterward.*)
      neighbors[states : {__state}] := Select[#, validQ] &@
         DeleteDuplicates@Flatten@
           Table[Module[{st2 = st},
      
             If[Length[st2[[from]]] > 0,
              PrependTo[st2[[to]], st2[[from, 1]]];
              st2[[from]] = Rest[st2[[from]]]];
      
             If[st2 =!= st && validQ[st2],
              Sow@UndirectedEdge[st, st2]];
      
             st2], {st, states}, {to, Length[st]}, {from, Length[st]}];
      
      hanoiGraph[s_, options___] := Module[{vertices, edges, n},
          n = Count[s, _Integer, Infinity];
          {vertices, {edges}} = Reap[Nest[neighbors, {s}, 2^n]];
      
          SetAttributes[UndirectedEdge, Orderless];
          Graph[DeleteDuplicates[edges]]];
      
    2. toStyle3D[g_] := Module[{st = VertexList[g][[#2]]},
          Tooltip[{ColorData[3][1 + Length[st] - Count[st, {}]],
            Opacity[1], Sphere[#1, .045]}, MatrixForm /@ List @@ st]] &;
      
      hanoiGraph3D[s_, options___] := Module[{g = hanoiGraph[s]},
         GraphPlot3D[g,
          Method -> "SpringElectricalEmbedding",
          VertexRenderingFunction -> toStyle3D[g], options, Boxed -> False,
          PlotStyle -> {Lighter[Blue](*,Opacity[.5]*)}]];
      
      {vv, vp} = {{0, 0, 1}, {2, 0, 0}};
      Animate[
       hanoiGraph3D[state[{}, {}, {}, Range[4]],
        Lighting -> "Neutral", SphericalRegion -> True,
        ViewVertical -> Dynamic[vv], Boxed -> False,
        ViewPoint -> Dynamic[RotationTransform[\[Theta], vv][vp], (vp = #1) &]],
       {\[Theta], 2 Pi, 0}, SynchronousUpdating -> False]
      

Although the 3-stick Hanoi graphs merely resemble Sierpinski graphs, it would be folly to ignore that resemblance given the thread of recursion that runs through both. We can create Sierpinski graphs easily, by once again reusing our polygon-in-polygon approach and this time replacing the Polygon[{p1, p2, p3}] expression with {p1 <-> p2, p2 <-> p3, p3 <-> p1}:

  1. axiom = Polygon[{{0, 0}, {1, Sqrt[3]}/2, {1, 0}}];
    
    next[prev_] := prev /. Polygon[{p1_, p2_, p3_}] :> {
         Polygon[{p1, (p1 + p2)/2, (p1 + p3)/2}],
         Polygon[{p2, (p2 + p3)/2, (p1 + p2)/2}],
         Polygon[{p3, (p1 + p3)/2, (p2 + p3)/2}]};
    
    draw[n_] := Graph@Flatten@
        Nest[next, N@axiom, n] /. Polygon[pts_] :>
         Sequence @@ UndirectedEdge @@@ Partition[pts, 2, 1, 1];
    
  2. axiom = Polygon[{{0, 0}, {1, Sqrt[3]}/2, {1, 0}}];
    
    next[prev_] := prev /. Polygon[{p1_, p2_, p3_}] :> {
         Polygon[{p1, (p1 + p2)/2, (p1 + p3)/2}],
         Polygon[{p2, (p2 + p3)/2, (p1 + p2)/2}],
         Polygon[{p3, (p1 + p3)/2, (p2 + p3)/2}]};
    
    (*Orderless attribute not necessary here. but it causes the particular permutation
     of the edge list that results in the particular layout*)
    
    (*triple-click on "DynamicSetting" below, Right-click -> Evaluate in Place*)
    DynamicSetting[SetterBar[1, {SetAttributes, ClearAttributes}]][UndirectedEdge, Orderless];
    draw[n_] := Graph[#, VertexSize -> .05, GraphLayout -> "SpringEmbedding"] &@
        Flatten@Nest[next, N@axiom, n] /. Polygon[pts_] :>
         Sequence @@ UndirectedEdge @@@ Partition[pts, 2, 1, 1];
    
    g = draw[5]
    GraphPlot3D[g, VertexRenderingFunction -> None,
       PlotStyle -> Hue[2/3, 2/3, 2/3(*,1/2*)],
       Method -> "SpringEmbedding", Boxed -> False]
    

There's the Sierpinski triangle I know and love; the graph of. You might think it doesn't look good. But you don't realize it's a Sierpinski triangle wearing a cape made of Sierpinski triangles. Not only does it not not look good, it looks completely badass. Because we're using the coordinates of the points as vertices, we can straightforwardly recover the regular Sierpinski layout:

  1. axiom = Polygon[{{0, 0}, {1, Sqrt[3]}/2, {1, 0}}];
    
    next[prev_] := prev /. Polygon[{p1_, p2_, p3_}] :> {
         Polygon[{p1, (p1 + p2)/2, (p1 + p3)/2}],
         Polygon[{p2, (p2 + p3)/2, (p1 + p2)/2}],
         Polygon[{p3, (p1 + p3)/2, (p2 + p3)/2}]};
    
    draw[n_] := Module[{edges},
       edges = Flatten@Nest[next, N@axiom, n] /. Polygon[pts_] :>
          Sequence @@ UndirectedEdge @@@ Partition[pts, 2, 1, 1];
    
       Graph[edges, VertexCoordinates -> VertexList[Graph[edges]],
        VertexSize -> .25]];
    

The point of spending 1 or 2 LOC's worth of developer time to convert our geometric Sierpinski triangle into a graph is so that we can ask questions about the graph. Like for example, what are its Hamiltonicness and Eulerity quotients? What is the average degree of the graph, in Celsius? In Kelvin? Frankly most of these questions are boring, and I don't really know anything about graphs. But here is a picture of the line graphs of the first few Sierpinski iterations:

  1. axiom = Polygon[{{0, 0}, {1, Sqrt[3]}/2, {1, 0}}];
    
    next[prev_] := prev /. Polygon[{p1_, p2_, p3_}] :> {
         Polygon[{p1, (p1 + p2)/2, (p1 + p3)/2}],
         Polygon[{p2, (p2 + p3)/2, (p1 + p2)/2}],
         Polygon[{p3, (p1 + p3)/2, (p2 + p3)/2}]};
    
    draw[n_] := Module[{edges},
       edges = Flatten@Nest[next, N@axiom, n] /. Polygon[pts_] :>
          Sequence @@ UndirectedEdge @@@ Partition[pts, 2, 1, 1];
    
       Graph[edges, VertexCoordinates -> VertexList[Graph[edges]],
        VertexSize -> .25]];
    
    g = draw[2];
    LineGraph[g]
    
    cycle = RandomChoice[{FindHamiltonianCycle, FindEulerianCycle}][g][[1]];
    
    Animate[
     HighlightGraph[g, Graph[cycle[[1 ;; n]]],
      EdgeShapeFunction -> (Line[#1] &),
      VertexShapeFunction -> None,
      GraphHighlightStyle -> "DehighlightHide"],
     {n, 1, Length[cycle], 1}, AnimationRate -> 1]
    
  2. axiom = Polygon[{{0, 0}, {1, Sqrt[3]}/2, {1, 0}}];
    
    next[prev_] := prev /. Polygon[{p1_, p2_, p3_}] :> {
         Polygon[{p1, (p1 + p2)/2, (p1 + p3)/2}],
         Polygon[{p2, (p2 + p3)/2, (p1 + p2)/2}],
         Polygon[{p3, (p1 + p3)/2, (p2 + p3)/2}]};
    
    sierpinskiGraph[n_] := Graph@Flatten@
         Nest[next, N@axiom, n] /. Polygon[pts_] :>
        Sequence @@ UndirectedEdge @@@ Partition[pts, 2, 1, 1];
    
    dg = Mean[VertexDegree[sierpinskiGraph[8]]];
    
    (*if live in us, above comes out in farenheit, so have to convert*)
    us = Graphics[CountryData["UnitedStates", "Polygon"], ImageSize -> 8000];
    If[Rasterize[us] === Rasterize[Show[us, Graphics[{PointSize[0], Point[Reverse@FindGeoLocation[]]}]]],
     WolframAlpha[ToString[dg, InputForm] <> " degrees f in celcius", {{"Result", 1}, "NumberData"}],
     dg]
    

Also the minimal zig zag of the triangle, notable because it looks like a bunch of resistors (no doubt the inspiration for certain papers). And its minimal criss cross. I don't really see anything though. Do you see anything? I don't see anything. These graphs just vertices and edges to me.

They do raise a question though. What game (or what anything) does the Sierpinski graph represent? I wasn't able to produce the Sierpinski triangle from any variation of the Hanoi game beyond the first couple of trivial iterations. In any case, through the extensive research I've done here I've found that layered graph layouts are pretty:

  1. validQ[s_state] := And @@ Less @@@ s;
    neighbors[states : {__state}] := Select[#, validQ] &@
       DeleteDuplicates@Flatten@
         Table[Module[{st2 = st},
           If[Length[st2[[from]]] > 0,
            PrependTo[st2[[to]], st2[[from, 1]]];
            st2[[from]] = Rest[st2[[from]]]];
           If[st2 =!= st && validQ[st2],
            Sow@UndirectedEdge[st, st2]];
           st2], {st, states}, {to, Length[st]}, {from, Length[st]}];
    
    hanoiGraph[s_] := Module[{vertices, edges, n},
       n = Count[s, _Integer, Infinity];
       {vertices, {edges}} = Reap[Nest[neighbors, {s}, 2^n]];
       SetAttributes[UndirectedEdge, Orderless];
       Graph[DeleteDuplicates[edges]]];
    
    axiom = Polygon[{{0, 0}, {1, Sqrt[3]}/2, {1, 0}}];
    next[prev_] := prev /. Polygon[{p1_, p2_, p3_}] :> {
         Polygon[{p1, (p1 + p2)/2, (p1 + p3)/2}],
         Polygon[{p2, (p2 + p3)/2, (p1 + p2)/2}],
         Polygon[{p3, (p1 + p3)/2, (p2 + p3)/2}]};
    
    (*certain layout depends on this ordering*)
    (*next[prev_]:=prev/.Polygon[pts_]:>(
    Polygon[ScalingTransform[1/2{1,1},#][pts]]&/@pts);*)
    
    sierpinskiGraph[n_] := Graph@Flatten@
        Nest[next, N@axiom, n] /. Polygon[pts_] :>
         Sequence @@ UndirectedEdge @@@ Partition[pts, 2, 1, 1];
    
    draw[g_] := LayeredGraphPlot[g,
       EdgeRenderingFunction -> ({CapForm["Round"], Line[#]} &),
       VertexRenderingFunction -> None,
       PlotStyle -> {Thickness[.01], Black}];
    
    draw /@ {hanoiGraph[state[{}, {}, Range[3]]], sierpinskiGraph[3]}
    

Chaos

One of the nice things about the chaos game algorithm is that we can easily generalize it to more than three points. To begin with, we can place equiangular points on a circle using $\cos$ and $\sin$ (see also my screwing around with polygons).

  1. draw[v_, numPoints_] := Module[{vertices},
      vertices = {Cos[#], Sin[#]} & /@ (2 Pi Range[v]/v);
    
      Graphics[{PointSize[0], Opacity[.1],
        Point[FoldList[(#1 + #2)/2 &, N@{0, 0},
          RandomChoice[N@vertices, numPoints]]]}]];
    

These are drawn with 10 million points. The last two are drawn with 50 million points. The key to the quality here is giving the points transparency so that varying degrees of overlap/nearness form different shades. Higher vertex counts clearly have some structure, but it becomes blurry for one reason or another. You might be able to pull out the structure better with a more methodical approach and some image trickery.

If you play around with pentagons in a vector editor (Mathematica itself has basic vector editing capabilities), you will find this figure:

I've highlighted one of the inner pentagons. You can see that this figure reproduces the faded stellation pattern in the center of the chaos game rendition. So the chaos game algorithm remains consistent in this geometric fashion: At each vertex of the figure, attach a copy of the larger figure, but with sidelength one-half of the original (note the red edge in the above image).

This also explains why the 4-vertex rendering is a block. And since we now have the geometric rule, we can turn to an explicit geometric construction to see if we can make the structure of these chaos games clearer. After some hiccups, I was able to get something working:

    1. SetAttributes[toXY, Listable];
      toXY[z_] := {Re[z], Im[z]};
      
      ring[c_, r_, 0] := c;
      ring[c_, r_, depth_] := Module[{zs},
         zs = c + r E^(I 2 Pi Range[3]/3);
         ring[c + # Normalize[# - c], r/2., depth - 1] & /@ zs];
      
      Graphics[Rotate[{Opacity[.95], LightGray, EdgeForm[Black],
         Polygon /@ toXY /@ Level[ring[0, 1, 5], {-2}]}, Pi]]
      
  1. draw[v_, n_] := Module[{ring},
       ring[c_, r_, depth_] := Module[{ps},
         ps = c + r {Cos[#], Sin[#]} & /@ (2. Pi Range[v]/v);
    
         If[depth == 0, Polygon[ps],
          ring[(c + #)/2, r/2, depth - 1] & /@ ps]];
    
       Graphics[{Transparent,
         EdgeForm[{Opacity[.28], Black}],
         ring[0., 1., n]}]];
    

The snowflake has all sorts of symmetries, probably because $6=2 \times 3$. It even has 3D grids and cubes. It's an infinite cubic matryoshka snowflake. And there is a lot of amazing detail in these drawings.

At this point I should mention that all of the code snippets on this page are self-contained. If you have Mathematica you can copy-paste this and start producing these figures.

The chaos game has another generalization. Instead of moving halfway between the active point and the randomly-chosen vertex, we can move 1/3rd of the way, or 3/2 of the way, etc:

  1. draw[v_, df_, numPoints_: 10000] := Module[{vertices},
       vertices = {Cos[#], Sin[#]} & /@ (2 Pi Range[v]/v);
    
       Graphics[{PointSize[0], Opacity[.5],
         Point[FoldList[df, N@{0, 0},
           RandomChoice[N@vertices, numPoints]]]}]];
    
    functions = Function[r, (#1 + #2) r &] /@ {1, .96, .7, .6, .5, .2};
    
    Grid[Join[
      {TraditionalForm[#[a, b]] & /@ functions},
      Table[draw[v, df], {v, 3, 6}, {df, functions}]]]
    

In the case where we're just adding the numbers, we get a normal n-directional random walk. Of course, the geometric approach has its own similar generalization:

  1. draw[v_, df_, n_] := Module[{ring},
       ring[c_, r_, depth_] := Module[{ps},
         ps = c + r {Cos[#], Sin[#]} & /@ (2. Pi Range[v]/v);
    
         If[depth == 0, Polygon[ps],
          ring[(c + #)/2, df[0, r], depth - 1] & /@ ps]];
    
       Graphics[{Transparent,
         EdgeForm[{Opacity[.26], Black}],
         ring[0., 1., n]}]];
    
    functions = Function[r, (#1 + #2) r &] /@ {1, .7, .6, .5, .35, .2};
    
    Grid[Join[
      {TraditionalForm[#[a, b]] & /@ functions},
      Table[draw[v, df, 4], {v, 3, 6}, {df, functions}]]]
    

One of the things you might try to do, if you're me, is adjust the ratio until the corners match up:

  1. drawGeom[v_, ratio_, n_] := Module[{ring},
       ring[c_, r_, depth_] := Module[{ps},
         ps = c + r {Cos[#], Sin[#]} & /@ (2. Pi Range[v]/v);
    
         If[depth == 0, Polygon[ps],
          ring[(c + #)/2, r ratio, depth - 1] & /@ ps]];
    
       Graphics[{Transparent,
         EdgeForm[{Opacity[.28], Black}],
         ring[0., 1., n]}]];
    
    drawChaos[v_, ratio_, numPoints_] := Module[{vertices},
       vertices = {Cos[#], Sin[#]} & /@ (2 Pi Range[v]/v);
    
       Graphics[{PointSize[0], Opacity[.5],
         Point[FoldList[(#1 + #2) ratio &, N@{0, 0},
           RandomChoice[N@vertices, numPoints]]]}]];
    
    With[{
      verticesC = Control[{vertices, 3, 8, 1, ImageSize -> Tiny}],
      iterationsC = Control[{iterations, 0, 8, 1, ImageSize -> Tiny}],
      numPointsC = Control[{numpoints, 0, 100000, 1, ImageSize -> Tiny}]},
    
     Manipulate[Overlay[{
        drawGeom[vertices, ratio, iterations],
        drawChaos[vertices, ratio, numpoints]}],
      Row[{verticesC, iterationsC, numPointsC}, " "],
      {{ratio, .5}, 0, 1, Appearance -> "Labeled"},
      Alignment -> Center]]
    

Look! There are Koch snowflake figures that form in the negative space. The boundary becomes snowflaked. A Koch snowflake can easily be made with an L-system construction:

    1. axiom = {F, right[2], F, right[2], F};
      rules = F -> {F, left[1], F, right[2], F, left[1], F};
      conversions = {F -> forward, dir_[n_] :> ConstantArray[dir, n]};
      
      (*state transformations*)
      forward[{z_, theta_}] := {z + E^(I theta), theta};
      left[{z_, theta_}] := {z, theta + 2 Pi/6};
      right[{z_, theta_}] := {z, theta - 2 Pi/6};
      
      draw[n_] :=
        Module[{program, zs},
         program = Flatten[Nest[# /. rules &, axiom, n] /. conversions];
         zs = First /@ ComposeList[program, {0., 0.}];
         Graphics[{Thin, Line[{Re[#], Im[#]} & /@ zs]}]];
      
      Grid[Partition[#, 3]] &[draw /@ Range[0, 5]]
      
    1. axiom = {F, right[1], F, right[1], F, right[1], F, right[1], F};
      rules = F -> {F, left[1], F, right[2], F, left[1], F};
      conversions = {F -> forward, dir_[n_] :> ConstantArray[dir, n]};
      
      (*state transformations*)
      forward[{z_, theta_}] := {z + E^(I theta), theta};
      left[{z_, theta_}] := {z, theta + 2 Pi/5};
      right[{z_, theta_}] := {z, theta - 2 Pi/5};
      
      draw[n_] :=
        Module[{program, zs},
         program = Flatten[Nest[# /. rules &, axiom, n] /. conversions];
         zs = First /@ ComposeList[program, {0., -Pi/10.}];
         Graphics[{Thin, Line[{Re[#], Im[#]} & /@ zs]}]];
      
      Grid[Partition[#, 3]] &[draw /@ Range[0, 5]]
      

With some minoradjustments we get our pentagonal snowflake. If we do the same procedure for the hexagonal chaos game we get the familiar triangular snowflake. All of the geometries seem to create Koch snowflakes, which makes sense given that indentations are triangles.

Of course, there are much more interesting generalizations we can come up with than simple ratios:

  1. draw[v_, df_, numPoints_: 1000] := Module[{vertices},
       vertices = {Cos[#], Sin[#]} & /@ (2 Pi Range[v]/v);
    
       Graphics[{PointSize[0], Opacity[.15],
         Point[FoldList[df, N@{0, 0},
           RandomChoice[N@vertices, numPoints]]]}]];
    
    rotate = RotationTransform;
    functions = {
       (#1 + #2)/RandomChoice[Prime[Range[3]]] &,
       (#1 + #2)/RandomChoice[Prime[Range[3]]!] &,
       (#1 + #2)/RandomChoice[Prime[Range[10]]] &,
       #1 + .5 rotate[10. Degree, #1][#2 - #1] &};
    
    Grid[Join[
      {TraditionalForm[Trace[#[a, b]][[2]]] & /@ functions},
      ParallelTable[draw[v, df], {v, 3, 5}, {df, functions}]]]
    

Some of these drawings remind me of the kind of fractal scattering found in the more deterministic algorithms. I wonder what kind of relation there is. The best distance function I found was logarithm-based:

  1. game = Compile[{{v, _Integer}, {wowzerz, _Real}, {numPoints, _Integer}},
       Module[{diff, vertices},
        vertices = {Cos[#], Sin[#]} & /@ (2 Pi Range[v]/v);
    
        (*distance function is Clip[(a+b)Log[EuclideanDistance[a,b]+wowzerz]*)
        FoldList[(
           diff = #2 - #1;(*note each of these is an x-y pair*)
           Clip[(#1 + #2) Log[Sqrt[diff.diff] + wowzerz], 1.1 {-2, 2}]) &,
         {0, 0}, RandomChoice[vertices, numPoints]]]];
    
    Graphics[{PointSize[0], Opacity[.08],
      Point[game[5, .8, 300000]]}(*,PlotRange->1.15*)]
    

All of these images are from the same distance function. The 'holes' on the inward-folded leaves of this one are interesting. It's like a fractal Klein bottle thing goin on there. If my computer was worth more than my car, as it some day will be, I would burn a lot of lightning-sequestered power in my mad scientist laboratory in the process of rendering different distance functions. There's a lot of pretty pictures in these simple chaos games. As it stands all this lightning is going to waste.

The geometric approach, not one to have been served, decides to go Tron:

    1. 
      
      draw[v_, df_, n_] := Module[{ring},
         ring[c_, r_, depth_] := Module[{zs},
           zs = c + r E^(I 2. Pi Range[v]/v);
      
           If[depth == 0, Polygon[{Re[#], Im[#]} & /@ zs],
            ring[(c + #)/2, df[0, r], depth - 1] & /@ zs]];
      
         Framed[Graphics[{Transparent,
            EdgeForm[{Thick, LightBlue}],
            ring[0., 1., n]}], FrameStyle -> LightBlue]];
      
      functions = {
         (#1 + #2)/RandomChoice[Prime[Range[3]]] &,
         (#1 + #2)/RandomChoice[Prime[Range[3]]!] &,
         (#1 + #2)/RandomChoice[Prime[Range[10]]] &,
         #1 + (1/2) (#2 - #1) E^(I 10. Degree) &};
      
      Framed[Grid[Join[
         {TraditionalForm[Trace[#[a, b]][[2]]] & /@ functions},
         Table[draw[v, df, 3], {v, 3, 5}, {df, functions}]]],
       Background -> Black, BaseStyle -> LightBlue]
      
      
      
    1. 
      
      
      
      
      draw[v_, df_, n_] := Module[{ring},
         ring[c_, r_, depth_] := Module[{ps},
           ps = c + r {Cos[#], Sin[#]} & /@ (2. Pi Range[v]/v);
      
           If[depth == 0, Polygon[ps],
            ring[(c + #)/2, df[0, r], depth - 1] & /@ ps]];
      
         Graphics[{EdgeForm[White], Opacity[.4],
           RGBColor[.4, 1, 1], ring[0., 1., n]}]];
      
      Show[(*repeatedly draw to cover more possibilities*)
       draw[4, (#1 + #2)/RandomChoice[Prime[Range[4]]] &,
          RandomChoice[{.1, 1.5} -> {2, 3}]] & /@ Range[20],
       Background -> Black, ImageSize -> 600]
      
      
      
      
      
      
      

Or Asteroids. Same thing.

The most interesting place I've seen the chaos game is in genetics. The idea is that instead of randomly picking the vertex at each step, you let the letters of the genetic code pick for you. There are 4 letters in DNA: A, T, G, C. So you run a chaos game with 4 vertices. If some sequence of DNA is AAAATC, your active point will approach the point labeled A 4 times, then it will approach the T point, then the C point.

If the DNA sequence is completely random, you will just recreate our beautiful block, which I have named the Charcoal Diamond:

But what you get is not random, as this chaos plot shows:

  1. coords = N@{"A" -> {1, -1}, "T" -> {-1, -1}, "G" -> {-1, 1}, "C" -> {1, 1}};
    dat = Characters[GenomeData[{"ChromosomeX", {28000000, 36000000}}]] /. coords;
    
    draw[data_] := Graphics[{PointSize[0], Opacity[.01],
        Point[FoldList[(#1 + #2)/2 &, N@{0, 0}, data]]},
       Epilog -> Style[Text @@@ coords, Red, Background -> White]];
    
    draw[dat]
    

This is a chaos game plot of an arbitrarily-chosen 8 million basepair sequence from our chromosome X (for scale, a typical protein is encoded in only a few hundred basepairs). You might insensibly think this happens because the letters occur with different frequencies, but that's not the case. The following is a chaos game plot of a sequence that was randomly generated over the same statistical frequencies as the above sequence:




  1. chars = Characters[GenomeData[{"ChromosomeX", {28000000, 36000000}}]];
    tallies = Tally[chars];
    BarChart[Last /@ tallies, ChartLabels -> First /@ tallies]
    
    1. coords = N@{"A" -> {1, 1}, "T" -> {-1, -1}, "G" -> {-1, 1}, "C" -> {1, -1}};
      dat = Characters[GenomeData[{"ChromosomeX", {28000000, 36000000}}]] /. coords;
      
      scale = 2;
      grid = Tuples[Range[{-1, -1}, {1, 1} - 2^-scale, 2^-scale]];
      color = If[Mod[Plus @@ #*2^scale, 2] == 1,
          Blend[{Lighter@Purple, Yellow}],
          Blend[{Lighter@Blue, Red}]] &;
      
      overunder = Graphics[{Opacity[.25],
          {color[#], Rectangle[#, # + 2^-scale]} & /@ grid}];
      
      draw[data_] := Show[overunder,
         Graphics[{PointSize[0], Opacity[.06],
           Point[FoldList[(#1 + #2)/2 &, N@{0, 0},
             RandomChoice[data, Length[data]]]]}],
         overunder, PlotRange -> 1, ImageSize -> {600, 600}];
      
      draw[dat] // Rasterize
      
  2. coords = N@{"A" -> {1, -1}, "T" -> {-1, -1}, "G" -> {-1, 1}, "C" -> {1, 1}};
    dat = Characters[GenomeData[{"ChromosomeX", {28000000, 36000000}}]] /. coords;
    
    draw[data_] := Graphics[{PointSize[0], Opacity[.01],
        Point[FoldList[(#1 + #2)/2 &, N@{0, 0}, RandomChoice[data, Length[data]]]]},
       Epilog -> Style[Text @@@ coords, Red, Background -> White]];
    
    draw[dat]
    

The letters do occur in different frequencies, but that doesn't make any interesting patterns. If you move the letters around you get a pattern related in this case to the fact that the frequencies are bilateral, but otherwise it's just a glorified chessboard. And look what happens when we do the same vertex movearounding for our genetic code:

  1. coords = N@{"A" -> {1, 1}, "T" -> {-1, -1}, "G" -> {-1, 1}, "C" -> {1, -1}};
    dat = Characters[GenomeData[{"ChromosomeX", {28000000, 36000000}}]] /. coords;
    
    draw[data_] := Graphics[{PointSize[0], Opacity[.01],
        Point[FoldList[(#1 + #2)/2 &, N@{0, 0}, data]]},
       Epilog -> Style[Text @@@ coords, Red, Background -> White]];
    
    draw[dat]
    

The original paper does a good job explaining how the chaos plot is something of a "fractal subsequence histogram." Assume your active point is anywhere in the entire square, and the next move is toward the bottom-left corner. Because you move halfway toward that corner (instead of, say, only one third or one fifth of the way), you will land inside that corner's quadrant regardless of where your point was to begin with.

Furthermore, you can apply this argument to subquadrants. It's easy to see this if you "work backwards." The formula for going from the current active point toward the next vertex is

$$p_{i+1}=\frac 1 2 (p_i + v)$$

By the DeLorean transform, we can go backwards like this:

$$p_{i-1}=2p_i-v$$

So, reversing all the points in a particular subquadrant:

  1. Manipulate[Module[{delorean, preimage, pt1, pt2},
      pt1 = ptc - width;
      pt2 = ptc + width;
    
      delorean[x_] := 2 x - coord;
      preimage = delorean /@ Tuples[Range[pt1, pt2, .025]];
    
      Graphics[{
        EdgeForm[{Thickness[.01], Darker[Gray, .4]}],
        {Transparent, Rectangle[{-1, -1}, {1, 1}]},
        {LightGray, Rectangle[pt1, pt2]},
        {Gray, Point[preimage]}},
       PlotRange -> 1.2, GridLines -> Automatic,
       GridLinesStyle -> Lighter[Gray, .8]]],
    
     {{width, .25}, 0, .5, 2.^-3},
     {{ptc, {-.75, -.25}}, Locator},
     {{coord, {-1, -1}}, Locator,
      Appearance -> Style["\[FilledSquare]", Red]}]
    

You can see here that for all the points which were just ordered to move toward the little red dot in the bottom-left corner, those that landed in the gray square had to have come from the top-left quadrant of the main square (the region with gray dots). So for the points in that gray region, we not only know that they were ordered to move toward the bottom-left vertex in the last step, but also that they were ordered to move toward the top-left vertex in the step before that.

And so on. The points in this little gray region were ordered to move previously toward the bottom-left, and before that the top-left, and before that the bottom-right. All points in that square have that history. So going back to our genetic chaos plot:

What those big holes mean is that CG is a rare sequence. As we just saw, a point can only get to that big empty square by coming from the bottom-right quadrant and going toward the top-left vertex. And since that square is so empty, there are rarely any points that are available to go toward other subsquares, such as this one, and so on.

This accounts for the texture of the first chaos plot as well. It just looks more wacked out because the CG vertices are adjacent, so the empty squares touch each other and create those staggered serrations. A simple histogram confirms our suspicions of CG Paucity — a.k.a. biology's Dark Energy:

  1. chars = Characters[GenomeData[{"ChromosomeX", {28000000, 36000000}}]];
    tallies = Sort[Tally[Partition[chars, 2, 1]]];
    BarChart[Last /@ tallies, ChartLabels -> CenterDot @@@ First /@ tallies]
    

If you sample subsequences instead of individual letters, and use those samples to simulate a genetic sequence, what's the smallest subsequence-sampling size you can get away with while still faithfully reproducing the texture of the chaos plot?

Asked differently, what length of subsequence is it that accounts for the texture of the chaos plot? Here is a graph of our DNA letters with pair-wise sequences labeled by probability:

  1. chars = Characters[GenomeData[{"ChromosomeX", {28000000, 36000000}}]];
    tallies = Sort[Tally[Partition[chars, 2, 1]]];
    sum = Total[Last /@ tallies];
    stats = {Rule @@ #1, #2/sum} & @@@ tallies;
    
    With[{r = .05},
     edgeF[pts_List, e_] := Arrow[pts, r];
     edgeF[pts_List, h_[a_, a_]] := Scale[Arrow[pts, r/.3], .3, pts[[1]]];
     vertexF = {EdgeForm[Opacity[.5]], Disk[#1, r], Darker[Gray, .7],
        Style[Text[#2, #1], 13, Bold, FontFamily -> "Comic Sans MS"]} &;
     edgeLabels = #1 -> Style[Round[#2, .01], 12, Bold] & @@@ stats;]
    
    Graph[First /@ stats, EdgeLabels -> edgeLabels,
     EdgeStyle -> Directive[{Thick, Opacity[.56]}],
     VertexShapeFunction -> vertexF, VertexStyle -> Orange,
     EdgeShapeFunction -> edgeF, PlotRangePadding -> .1]
    

This is a graph of what's called a Markov chain, but don't quote me on the formalities. (Mathematica 9 has built-in Markov whatitswhats, but I'm using version 8). The point is we can generate a sequence whose letter-to-letter statistics are the same as those of our original DNA by following the graph probaballistically:

  1. getStats[data_] := Module[{tallies, sum},
       tallies = Tally[Partition[data, 2, 1]];
       sum = Total[Last /@ tallies];
       {Rule @@ #1, #2/sum} & @@@ tallies];
    
    draw[data_, options___] := Graphics[{PointSize[0], Opacity[.01],
        Point[FoldList[(#1 + #2)/2 &, N@{0, 0}, data]]}, options,
       Epilog -> Style[Text @@@ coords, Red, Background -> White]];
    
    chars = Characters[GenomeData[{"ChromosomeX", {28000000, 36000000}}]];
    stats = getStats[chars];
    
    Do[With[{weights =
        Rule @@ Transpose@Cases[stats, {letter -> to_, p_} :> {N[p], to}]},
      next[letter] := RandomChoice[weights]],
     {letter, DeleteDuplicates[stats[[All, 1, 1]]]}]
    
    coords = N@{"A" -> {1, 1}, "T" -> {-1, -1}, "G" -> {-1, 1}, "C" -> {1, -1}};
    pseudoDat = NestList[next, "A", Length[chars]] /. coords;
    
    draw[pseudoDat]
    

You can see that, while similar, the fake plot immediately stands out as too Hollywood compared to the verisimilous beauty of the real data. The most notable distinction between them, besides the grain, is the dark diagonal that crosses A and T in the real plot, presumably because those two letters have a lot of interplay. That it's not replicated by our pseudosequence may mean there are a relatively large amount of ATA, TAT subsequences.

So it looks like subsequences of length 2 aren't sufficient. We could generalize our Markovizer, but what I think is actually interesting here is the grain. We can do some image processing to see if we can bring it out:




  1. ListPlot[FoldList[(#1 + #2)/2 &, 1,
      Mod[Range[1000], 2]], PlotRange -> 1,
     Ticks -> {Automatic, Range[0, 1, 1/3]}]
  2. getStats[data_] := Module[{tallies, sum},
       tallies = Tally[Partition[data, 2, 1]];
       sum = Total[Last /@ tallies];
       {Rule @@ #1, #2/sum} & @@@ tallies];
    
    draw[data_, options___] := Graphics[{PointSize[0], Opacity[.01],
        Point[FoldList[(#1 + #2)/2 &, N@{0, 0}, data]]}, options];
    
    chars = Characters[GenomeData[{"ChromosomeX", {28000000, 36000000}}]];
    stats = getStats[chars];
    
    Do[With[{weights =
        Rule @@ Transpose@Cases[stats, {letter -> to_, p_} :> {N[p], to}]},
      next[letter] := RandomChoice[weights]],
     {letter, DeleteDuplicates[stats[[All, 1, 1]]]}]
    
    coords = N@{"A" -> {1, 1}, "T" -> {-1, -1}, "G" -> {-1, 1}, "C" -> {1, -1}};
    pseudoDat = NestList[next, "A", Length[chars] - 1] /. coords;
    realDat = chars /. coords;
    
    With[{upsc = 2},
     pseudo = draw[pseudoDat, ImageSize -> upsc 600] // Rasterize;
     real = draw[realDat, ImageSize -> upsc 600] // Rasterize;
    
     With[{\[Theta] = ColorNegate},
       (ImageSubtract[\[Theta][real], \[Theta][pseudo]] // \[Theta])
         ~MinFilter~1
         ~ImageMultiply~1.1
        (*~ImageAdjust~(9!)*)
        // ImageAdjust]
      ~ImageResize~Scaled[1/upsc]]
    

My intuition here is to subtract the Hollywood plot from the real plot (as images) in order to highlight artifacts that are due to longer subsequence patterns. The two horizontal streaks are at the one-third and two-thirds marks of the square as a whole, which I think implies a lot of AT/TA. Note that a point at $\small{\frac 1 3}$ is halfway between zero and $\small{\frac 2 3}$, and vice versa. $\small{\{\frac 1 3, \frac 2 3\}}$ is the fixed point of alternating T/A, so to speak. Here's a histogram of 3-sequences:

    1. chars = Characters[GenomeData[{"ChromosomeX", {28000000, 36000000}}]]; 
      tallies = Sort[Tally[Partition[chars, 3, 1]]];
      BarChart[Last /@ tallies, ChartLabels -> Column /@ First /@ tallies]
      
      
    1. chars = Characters[GenomeData[{"ChromosomeX", {28000000, 36000000}}]];
      tallies = Sort[Tally[Partition[chars, 5, 1]]];
      ListPlot[Cases[tallies, {seq_, count_} :> Tooltip[count, Column[seq]]],
       Axes -> None, Filling -> Axis, PlotRange -> Full]
      

Which actually just shows a lot of TTT and AAA. Longer subsequence statistics show a similar picture. And of course we can always just do this:

  1. $$ \begin{array}{llllllll} \text{A} & 1119676 & \text{T} & 1120986 & \text{G} & 882866 & \text{C} & 882887 \\ \text{AA} & 316677 & \text{TT} & 322527 & \text{GG} & 195417 & \text{CC} & 199940 \\ \text{AAA} & 127806 & \text{TTT} & 130977 & \text{GGG} & 45560 & \text{CCC} & 46320 \\ \text{AAAA} & 51308 & \text{TTTT} & 52674 & \text{GGGG} & 10022 & \text{CCCC} & 9962 \\ \text{AAAAA} & 18837 & \text{TTTTT} & 19401 & \text{GGGGG} & 2177 & \text{CCCCC} & 2199 \\ \text{AAAAAA} & 5617 & \text{TTTTTT} & 5772 & \text{GGGGGG} & 468 & \text{CCCCCC} & 409 \\ \text{AAAAAAA} & 2308 & \text{TTTTTTT} & 2326 & \text{GGGGGGG} & 69 & \text{CCCCCCC} & 59 \\ \text{AAAAAAAA} & 803 & \text{TTTTTTTT} & 797 & \text{GGGGGGGG} & 18 & \text{CCCCCCCC} & 14 \\ \text{AAAAAAAAA} & 419 & \text{TTTTTTTTT} & 458 & \text{GGGGGGGGG} & 4 & \text{CCCCCCCCC} & 8 \\ \text{AAAAAAAAAA} & 259 & \text{TTTTTTTTTT} & 265 & \text{GGGGGGGGGG} & 3 & \text{CCCCCCCCCC} & 8 \\ \text{AAAAAAAAAAA} & 141 & \text{TTTTTTTTTTT} & 173 & \text{GGGGGGGGGGG} & 2 & \text{CCCCCCCCCCC} & 12 \\ \text{AAAAAAAAAAAA} & 125 & \text{TTTTTTTTTTTT} & 114 & \text{GGGGGGGGGGGG} & 1 & \text{CCCCCCCCCCCC} & 2 \\ \text{AAAAAAAAAAAAA} & 83 & \text{TTTTTTTTTTTTT} & 112 & \text{GGGGGGGGGGGGG} & 1 & \text{CCCCCCCCCCCCC} & 1 \\ \text{AAAAAAAAAAAAAA} & 59 & \text{TTTTTTTTTTTTTT} & 110 & \text{GGGGGGGGGGGGGG} & 1 & \text{CCCCCCCCCCCCCC} & 2 \\ \text{AAAAAAAAAAAAAAA} & 63 & \text{TTTTTTTTTTTTTTT} & 82 & \text{GGGGGGGGGGGGGGG} & 1 & \text{CCCCCCCCCCCCCCC} & 1 \\ \text{AAAAAAAAAAAAAAAA} & 62 & \text{TTTTTTTTTTTTTTTT} & 66 & \text{GGGGGGGGGGGGGGGG} & 2 & \text{CCCCCCCCCCCCCCCCC} & 1 \\ \text{AAAAAAAAAAAAAAAAA} & 34 & \text{TTTTTTTTTTTTTTTTT} & 52 & \text{GGGGGGGGGGGGGGGGGGGGG} & 1 & \text{} & \text{} \\ \text{AAAAAAAAAAAAAAAAAA} & 43 & \text{TTTTTTTTTTTTTTTTTT} & 49 & \text{} & \text{} & \text{} & \text{} \\ \text{AAAAAAAAAAAAAAAAAAA} & 28 & \text{TTTTTTTTTTTTTTTTTTT} & 38 & \text{} & \text{} & \text{} & \text{} \\ \text{AAAAAAAAAAAAAAAAAAAA} & 28 & \text{TTTTTTTTTTTTTTTTTTTT} & 31 & \text{} & \text{} & \text{} & \text{} \\ \text{AAAAAAAAAAAAAAAAAAAAA} & 34 & \text{TTTTTTTTTTTTTTTTTTTTT} & 28 & \text{} & \text{} & \text{} & \text{} \\ \text{AAAAAAAAAAAAAAAAAAAAAA} & 31 & \text{TTTTTTTTTTTTTTTTTTTTTT} & 23 & \text{} & \text{} & \text{} & \text{} \\ \text{AAAAAAAAAAAAAAAAAAAAAAA} & 26 & \text{TTTTTTTTTTTTTTTTTTTTTTT} & 14 & \text{} & \text{} & \text{} & \text{} \\ \text{AAAAAAAAAAAAAAAAAAAAAAAA} & 15 & \text{TTTTTTTTTTTTTTTTTTTTTTTT} & 23 & \text{} & \text{} & \text{} & \text{} \\ \text{AAAAAAAAAAAAAAAAAAAAAAAAA} & 11 & \text{TTTTTTTTTTTTTTTTTTTTTTTTT} & 12 & \text{} & \text{} & \text{} & \text{} \\ \text{AAAAAAAAAAAAAAAAAAAAAAAAAA} & 12 & \text{TTTTTTTTTTTTTTTTTTTTTTTTTT} & 15 & \text{} & \text{} & \text{} & \text{} \\ \text{AAAAAAAAAAAAAAAAAAAAAAAAAAA} & 15 & \text{TTTTTTTTTTTTTTTTTTTTTTTTTTT} & 8 & \text{} & \text{} & \text{} & \text{} \\ \text{AAAAAAAAAAAAAAAAAAAAAAAAAAAA} & 11 & \text{TTTTTTTTTTTTTTTTTTTTTTTTTTTT} & 7 & \text{} & \text{} & \text{} & \text{} \\ \text{AAAAAAAAAAAAAAAAAAAAAAAAAAAAA} & 7 & \text{TTTTTTTTTTTTTTTTTTTTTTTTTTTTT} & 3 & \text{} & \text{} & \text{} & \text{} \\ \text{AAAAAAAAAAAAAAAAAAAAAAAAAAAAAA} & 6 & \text{TTTTTTTTTTTTTTTTTTTTTTTTTTTTTT} & 5 & \text{} & \text{} & \text{} & \text{} \\ \text{AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA} & 5 & \text{TTTTTTTTTTTTTTTTTTTTTTTTTTTTTTT} & 2 & \text{} & \text{} & \text{} & \text{} \\ \text{AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA} & 3 & \text{TTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTT} & 3 & \text{} & \text{} & \text{} & \text{} \\ \text{AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA} & 1 & \text{TTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTT} & 3 & \text{} & \text{} & \text{} & \text{} \\ \text{AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA} & 4 & \text{TTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTT} & 3 & \text{} & \text{} & \text{} & \text{} \\ \text{AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA} & 5 & \text{TTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTT} & 1 & \text{} & \text{} & \text{} & \text{} \\ \text{AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA} & 1 & \text{TTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTT} & 2 & \text{} & \text{} & \text{} & \text{} \\ \text{AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA} & 2 & \text{TTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTT} & 2 & \text{} & \text{} & \text{} & \text{} \\ \text{AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA} & 2 & \text{TTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTT} & 1 & \text{} & \text{} & \text{} & \text{} \\ \text{AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA} & 1 & \text{TTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTT} & 1 & \text{} & \text{} & \text{} & \text{} \\ \text{AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA} & 1 & \text{} & \text{} & \text{} & \text{} & \text{} & \text{} \\ \end{array} $$
  2. 
    
    
    string = GenomeData[{"ChromosomeX", {28000000, 36000000}}];
    tallies = Sort[Tally[StringCases[string, ("A" .. | "T" .. | "C" .. | "G" ..)]]];
    
    Grid[Join[
       Sequence @@ Reverse@Sort@
          SplitBy[tallies, StringTake[First[#], 1] &], 2],
      Alignment -> Left] // Magnify[#, 1/2] &
    
    
    
    
  3. 
    
    
    string = GenomeData[{"ChromosomeX", {28000000, 36000000}}];
    strings = ParallelTable[Module[{cases =
          StringCases[string, Alternatives @@ cs .., Overlaps -> True]},
        Last@SortBy[cases, StringLength]],
       {cs, Subsets[{"A", "C", "T", "G"}, {2}]}];
    
    Grid[{Tooltip[Short[#], #], StringLength[#]} & /@
      Reverse@SortBy[strings, StringLength], Alignment -> Left]
    
    
    
    

The longest single-letter run length in this section of DNA is 48 As. The longest string of A-or-T is 222 basepairs long. Quite long, but the longest pairing is actually T/C which has a sequence of length 231. C/G's longest sequence is 34 basepairs long. I wonder what it is about CG. Maybe an unusually (un)useful amino acid or some hydrophobilia issue. I wonder too if these are blanket statistical patterns or if certain quirks are only present, say, in non-coding regions.

You might be wondering why we don't just ask a biologist about these mysteries. The reason is because you're inside a car right now, I'm driving, we're lost, both of us are tourists, and I'm one of those people that would sooner burn hours of gasoline/diesel than ask for directions. You also suspect I might be some kind of criminal, so you're afraid of bringing up the issue. All around it's pretty awkward in here.

We can do a lot better than these static diagrams by giving ourselves the ability to manually movearound the vertices to see if we can find interesting patterns:

  1. coords = N@{"A" -> {1, 1}, "T" -> {-1, -1}, "G" -> {-1, 1}, "C" -> {1, -1}};
    chars = Characters[GenomeData[{"ChromosomeX", {28000000, 36000000}}]];
    
    draw[data_, options___] := Graphics[{PointSize[0], Opacity[.1],
        Point[FoldList[(#1 + #2)/2 &, N@{0, 0}, data]]}, options];
    
    Manipulate[
     draw[
      chars[[1 ;; 100000]] /. Thread[First /@ coords -> pts],
      PlotRange -> 1.1],
     {{pts, Last /@ coords}, Locator,
      Appearance -> (Framed[#, BaseStyle -> Red,
           FrameStyle -> None, FrameMargins -> 0,
           Background -> White] & /@ First /@ coords)}]
    

And a tool that repeatedly applies the DeLorean transform to rebuild the sequence leading up to a region:

  1. coords = N@{"A" -> {1, 1}, "T" -> {-1, -1}, "G" -> {-1, 1}, "C" -> {1, -1}};
    dat = Characters[GenomeData[{"ChromosomeX", {28000000, 36000000}}]] /. coords;
    nfLetter = Module[{nf = Nearest[Reverse /@ coords]}, nf[#, 1][[1]] &];
    
    draw[data_, options___] := Graphics[{PointSize[0], Opacity[.01],
        Point[FoldList[(#1 + #2)/2 &, N@{0, 0}, data]]}, options];
    
    background = Raster[Reverse@ImageData@Rasterize@
          draw[dat, PlotRange -> 1, ImageSize -> 600 {1, 1}],
       {{-1, -1}, {1, 1}}];
    
    delorean[p_] := 2 p - (nfLetter[p] /. coords);
    Manipulate[Module[{seq, r = 2^radius},
      seq = Most[NestList[delorean, pt, Floor[1/(2 r)]]];
    
      Graphics[{background, {Darker@Gray, Point[seq]},
        {Orange, Thick, Opacity[.8], Arrow[Reverse[seq]],
         MapThread[Circle, {seq, 2^(radius + 1) r 2^Range[Length[seq]]}]}},
       PlotRange -> 1.015, PlotLabel -> (nfLetter /@ Reverse[seq]),
       Epilog -> Style[Text @@@ coords, Red, Background -> White]]],
    
     {{radius, -2}, -4, -2, 1},
     {{pt, {-.75, -.25}}, Locator}]
    

I'm not actually sure how legit the maths of the program are, but there it be. Let's return to our charcoal diamond, here rotated:

  1. drawDiamond[numPoints_] := Graphics[{PointSize[0], Opacity[.01],
        Point[FoldList[(#1 + #2)/2 &, N@{0, 0}, RandomInteger[{0, 1}, {numPoints, 2}]]]},
       ImageSize -> 600, PlotRangePadding -> 0];
    
    diamond = drawDiamond[15000000] // Rasterize;
    
    axiom = {{Transparent, Rectangle[Scaled[{0, 0}], Scaled[{2, 2}]]}, White,
       If[hl, EdgeForm[{Opacity[.1], Green}]], Rectangle[Scaled[{1, 1}], Scaled[{2, 2}]]};
    
    next[prev_] := Translate[Scale[prev, .5], {{-1, 1}, {-1, -1}, {1, -1}}];
    
    Control[{hl, {True, False}}] Control[{n, 0, 10, 1}]
    Dynamic[Overlay[{diamond, Graphics[NestList[next, axiom, n],
        ImageSize -> (ImageSize /. AbsoluteOptions[diamond]), PlotRange -> 1]}]]
    

Imagine we suddenly removed one vertex. That would mean that points can no longer land in that quadrant. Which would mean that no points could go from that quadrant to these subquadrants. Which would mean no points going to these subquadrants. And soon and soforth, until.

So that explains the holes in the Sierpinski triangle. I call this the "Sierpinski triangle by infinite quadrilateral descent" method of construction. It seems very natural to me, but it raises the question of what these regions in the various deterministic constructions have to do with each other:

  1. draw[n_] := Array[Tooltip[Mod[Binomial[##], 2],
         TraditionalForm[HoldForm[Binomial[##]] == Binomial[##]]] &, {2^n, 2^n}, 0];
    proc[a_ /; Length[a] == 2] := a;
    proc[arr_] := Module[{l = Length[arr]/2},
       ArrayFlatten@Map[Function[square,
          If[FreeQ[square, Tooltip[1, _]],
           (**)Map[Style[#, Bold, ColorData[3][l]] &, square, {2}],
           (**)proc[square]]], Partition[arr, {l, l}], {2}]];
    Style[MatrixForm[proc[draw[5]]], Background -> GrayLevel[.98]]
    

(To be clear, the chaos game is just an algorithmic tradeoff vs the geometric approach. It is not necessarily doing anything non-deterministic in the larger scheme.) In this case I think the parity/binary explanations are going to be the simplest, though I'm a math noob and I don't see an immediately obvious way of approaching this, if the question even makes sense in the way I seem to be implying. However with some inspiration we can find an iterative angle that seems to me like a kind of multiplication:


$$ \begin{array}{l} \left( \begin{array}{c} 1 \\ \end{array} \right)\to \alpha \\ \alpha \left( \begin{array}{cc} 1 & 0 \\ 1 & 1 \\ \end{array} \right)=\left( \begin{array}{cc} \alpha & 0 \\ \alpha & \alpha \\ \end{array} \right)=\left( \begin{array}{cc} \left( \begin{array}{c} 1 \\ \end{array} \right) & \left( \begin{array}{c} 0 \\ \end{array} \right) \\ \left( \begin{array}{c} 1 \\ \end{array} \right) & \left( \begin{array}{c} 1 \\ \end{array} \right) \\ \end{array} \right)=\left( \begin{array}{cc} 1 & 0 \\ 1 & 1 \\ \end{array} \right)\to \alpha \\ \alpha \left( \begin{array}{cc} 1 & 0 \\ 1 & 1 \\ \end{array} \right)=\left( \begin{array}{cc} \alpha & 0 \\ \alpha & \alpha \\ \end{array} \right)=\left( \begin{array}{cc} \left( \begin{array}{cc} 1 & 0 \\ 1 & 1 \\ \end{array} \right) & \left( \begin{array}{cc} 0 & 0 \\ 0 & 0 \\ \end{array} \right) \\ \left( \begin{array}{cc} 1 & 0 \\ 1 & 1 \\ \end{array} \right) & \left( \begin{array}{cc} 1 & 0 \\ 1 & 1 \\ \end{array} \right) \\ \end{array} \right)=\left( \begin{array}{cccc} 1 & 0 & 0 & 0 \\ 1 & 1 & 0 & 0 \\ 1 & 0 & 1 & 0 \\ 1 & 1 & 1 & 1 \\ \end{array} \right)\to \alpha \\ \end{array} $$

So the Sierpinski triangle is the infinith power of ${ \left( \begin{array}{cc} 1 & 0 \\ 1 & 1 \\ \end{array} \right)} $ under this 'multiplication.' Someone who knows enough group theory might recognize what's going on here. Unfortunately I don't, but one thing we can do is investigate this 'multiplication' in general. I was going to make a simple program to do that, but I got carried away and made this:

    1. (* minimal *)
      
      iterate[matrix_, power_] := Nest[ArrayFlatten[
           ConstantArray[#, Dimensions[matrix]] matrix] &, 1, power];
      
      draw[matrix_, power_] :=
        ArrayPlot[iterate[matrix, power],
         Frame -> False, PixelConstrained -> 1];
      
      draw[{{1, 0}, {1, 1}}, 10]
      
    2. matrixInput[Dynamic[m_], Dynamic[rot_]] :=
        Dynamic[Rotate[Deploy[MatrixForm[#, TableSpacing -> {0, 0}]], rot] &@
          Array[(*(*better performance*)Rotate[#,-rot]&@*)
           Checkbox[Dynamic[m[[##]]], {0, 1}] &, Dimensions[m]]];
      
      bg = White;
      dims = # -> If[# > 4, Style[#, Red], #] & /@ Range[7];
      
      iterate[matrix_, power_] := Nest[ArrayFlatten[
           ConstantArray[#, Dimensions[matrix]] matrix] &, 1, power];
      
      controls = With[{
          mC = Control[{{m, 2, ""}, dims, ControlType -> PopupMenu}],
          nC = Control[{{n, 2, ""}, dims, ControlType -> PopupMenu}],
          matrixInputC = matrixInput[Dynamic[matrix], Dynamic[rot]],
          colorC = Control[{{color, Black}, ColorSlider}],
          rotC = Control[{{rot, 0, "\[Theta]"}, Pi, -Pi, Pi/16}],
          powerC = Control[{{power, 3}, 1, 4, 1, Appearance -> "Labeled"}],
          opacityC = Control[{{opacity, 1}, 0, 1, ImageSize -> Small}],
          primitiveC = Control[{{primitive, Rectangle[]},
             (# -> Graphics[{color, #}, ImageSize -> 20] &) /@ {
               {PointSize[Tiny], Point[N@{0, 0}]},
               {EdgeForm[None], Disk[N@{0, 0}, .5]},
               Rotate[Scale[Rectangle[], 1./Sqrt[2]], Pi/4],
               Rectangle[]}, SetterBar}],
          backgroundC = Row[{"background   ",
             Framed[ColorSlider[Dynamic[background, (bg = background = #) &],
               AppearanceElements -> "Swatch"], FrameStyle -> Darker[Gray]],
              " ",
             ColorSlider[Dynamic[background, (bg = background = #) &],
              AppearanceElements -> "Spectrum", ImageSize -> Small]}]},
      
         Row[{
           Column[{
             Row[{mC, "   \[Times]", nC}],
             Row[{"    ", matrixInputC}]}],
           Spacer[40],
           Column[{colorC, rotC, powerC}],
           Column[{backgroundC, opacityC, primitiveC}]}]];
      
    3. Panel[#, Background -> Dynamic[bg]] &@
       Manipulate[
        If[{m, n} =!= Dimensions[matrix], matrix = PadRight[matrix, {m, n}]];
      
        With[{primitives = Rotate[#, rot - Pi/2] &@
            Translate[primitive, Position[
              iterate[matrix /. 0 matrix -> {{1}}, power], 1]]},
      
         Graphics[{Dynamic[EdgeForm[{Opacity[opacity], color}]],
           Dynamic[color], Dynamic[Opacity[opacity]], primitives},
          ImageSize -> {{400, Large}, {400, Large}},
          Background -> Dynamic[background]]],
      
        Evaluate[controls],
      
        (*declare variables here for persistence*)
        {{background, bg = White}, ControlType -> None},
        {{matrix, {{1, 0}, {1, 1}}}, ControlType -> None},
      
        Bookmarks :> {
          "Random" :> (matrix = RandomChoice[{.4, .6} -> {0, 1}, Dimensions[matrix]]),
          "Invert" :> (matrix = BitXor[matrix, 1]),
          "Array Print" :> (With[{p = power, m = matrix, c = color, o = opacity, bg = background},
            CellPrint[ExpressionCell[Defer[
               ArrayPlot[iterate[m, p], Frame -> False, PixelConstrained -> 1,
                ColorRules -> {0 -> bg, 1 -> c /. RGBColor[r_, g_, b_] :> RGBColor[r, g, b, o]}]],
              "Input"]]]),
          "Clear" :> (matrix = 0 matrix)},
      
        Paneled -> False, SynchronousUpdating -> Automatic,
        SaveDefinitions -> True, LabelStyle -> Darker[Gray], Alignment -> Center]
      

I know what this image reminds you of. Those little candle chandoliers that you hit in Castlevania to make hearts and morning stars come out. I also found The Sierpinski Scream, a letter H that would definitely beat you up if it was human, the up arrow and its Hot Topic-donning offspring, pink infinities made of pink infinities, even the vaunted Sierpinski Chronobracket.

Essentially what we have here in these little matrices is a notation for specifying translations. It's yet another algorithm with different tradeoffs for doing more or less the same thing that our chaos game and geometric algorithms are doing. We can bring this characteristic out by allowing arbitrary rules:

    1. (* minimal *)
      
      iterate[matrix_, power_, matrix1_: {{1}}] := Module[{rules =
           {0 -> (0 # &), 1 -> (# &), T -> Transpose,
            R -> (Transpose[Reverse[#]] &), L -> (Reverse[Transpose[#]] &)}},
      
         Nest[Function[prev,
           ArrayFlatten[Map[#[prev] &, matrix /. rules, {2}]]],
          matrix1, power]];
      
      draw[matrix_, power_] :=
        ArrayPlot[iterate[matrix, power],
         Frame -> False, PixelConstrained -> 1];
      
      draw[{{1, 0}, {T, R}}, 10]
      
    2. matrixInput1[Dynamic[m_], Dynamic[rot_]] :=
        Dynamic[Rotate[Deploy[MatrixForm[#, TableSpacing -> {0, 0}]], rot] &@
          Array[(*(*better performance*)Rotate[#,-rot]&@*)
           EventHandler[Checkbox[Dynamic[m[[##]]], {0, 1}],
             {"MouseDown", 2} :> (m[[##]] = 0)] &,
           Dimensions[m]], 0];
      
      matrixInput2[Dynamic[m_], Dynamic[rules_], Dynamic[color_],
         Dynamic[rot_]] :=
        With[{
          tooltip = Tooltip[#, "Click to cycle\nRight-click to zero", TooltipDelay -> .8] &,
          eatRightClick = EventHandler[#, {"MouseDown", 2} :> {}] &,
          matrixForm = MatrixForm[#, TableSpacing -> {1, 1}] &},
      
         Dynamic[
          eatRightClick@Style[#, color] &@
             Rotate[#, rot] &@tooltip@Deploy@matrixForm@
              Array[
               EventHandler[Toggler[Dynamic[m[[##]]], First /@ rules],
                 {"MouseDown", 2} :> (m[[##]] = 0)] &,
               Dimensions[m]]]];
      
      bg = White;
      dims = # -> If[# > 4, Style[#, Red], #] & /@ Range[4];
      defaultRules = {0 -> (0 # &), 1 -> (# &), T -> Transpose,
         R -> (Transpose[Reverse[#]] &), L -> (Reverse[Transpose[#]] &)};
      
      iterate[matrix_, matrix1_, rules_, power_] :=
        Nest[
         Function[prev, ArrayFlatten[Map[#[prev] &, matrix /. rules, {2}]]],
         matrix1, power];
      
    3. controls = With[{
          m1C = Control[{{m1, 2, ""}, dims, ControlType -> PopupMenu}],
          m2C = Control[{{m2, 2, ""}, dims, ControlType -> PopupMenu}],
          matrixInput1C = matrixInput1[Dynamic[matrix1], Dynamic[rot]],
          matrixInput2C = matrixInput2[Dynamic[matrix], Dynamic[rules], Dynamic[color], Dynamic[rot]],
          rulesC = OpenerView[{"Rules", Control[{{rules, defaultRules, ""}, InputField,
               Background -> Dynamic[Lighter[background, .65]], FieldSize -> {45, 5}}]}],
          colorC = Control[{{color, Black}, ColorSlider}],
          rotC = Control[{{rot, 0, "\[Theta]"}, Pi, -Pi, Pi/16}],
          powerC = Control[{{power, 3}, 1, 4, 1, Appearance -> "Labeled"}],
          backgroundC = Row[{"background   ",
             Framed[ColorSlider[Dynamic[background, (bg = background = #) &],
               AppearanceElements -> "Swatch"], FrameStyle -> Darker[Gray]], " ",
             ColorSlider[Dynamic[background, (bg = background = #) &],
              AppearanceElements -> "Spectrum", ImageSize -> Small]}],
          opacityC = Control[{{opacity, 1}, 0, 1, ImageSize -> Small}],
          primitiveC = Control[{{primitive, Rectangle[]},
             (# -> Graphics[{color, #}, ImageSize -> 20] &) /@ {
               {PointSize[Tiny], Point[{0, 0}]},
               {EdgeForm[None], Disk[{0, 0}, .5]},
               Rotate[Scale[Rectangle[], 1/Sqrt[2]], Pi/4],
               Rectangle[]}, SetterBar}]},
      
         Row[{
           Column[{
             Row[{m1C, "   |", m2C}],
             Row[{"    ", matrixInput1C, "  ", matrixInput2C}]}],
           Spacer[40],
           Column[{rulesC, OpenerView[#, True] &@
              {"Style", Row[{Column[{colorC, rotC, powerC}],
                 Column[{backgroundC, opacityC, primitiveC}]}]}}]}]];
      
      bookmarks = {
         "Random" :> (
           matrix1 = RandomChoice[{0, 1}, Dimensions[matrix1]];
           matrix = RandomChoice[First /@ rules, Dimensions[matrix]]),
      
         "Array Print" :> With[
           {m1 = matrix1, m = matrix, r = rules, p = power, c = color, o = opacity, bg = background},
           CellPrint[ExpressionCell[Defer[
              ArrayPlot[
               iterate[m /. 0 m -> {{1}}, m1 /. 0 m1 -> {{1}}, r, p], PixelConstrained -> 1, Frame -> False,
               ColorRules -> {0 -> bg, 1 -> c /. RGBColor[r_, g_, b_] :> RGBColor[r, g, b, o]}]],
             "Input"]]],
      
         "Clear" :> (matrix = 0 matrix)};
      
    4. Panel[#, Background -> Dynamic[bg]] &@
       Manipulate[
        If[{m1, m1} =!= Dimensions[matrix1], matrix1 = PadRight[matrix1, {m1, m1}]];
        If[{m2, m2} =!= Dimensions[matrix], matrix = PadRight[matrix, {m2, m2}]];
      
        (*remove rules from matrix that no longer exist*)
        Module[{matrixP, default = rules[[1, 1]]},
         matrixP = Replace[matrix, a_ /; ! MemberQ[First /@ rules, a] -> default, {2}];
         If[matrix =!= matrixP, matrix = matrixP]];
      
        With[{primitives =
           Rotate[Translate[primitive, Position[#, 1]], rot - Pi/2] &@
            iterate[
             matrix /. 0 matrix -> {{1}},
             matrix1 /. 0 matrix -> {{1}}, rules,
             ControlActive[Max[power - 2, 2], power]]},
      
         Graphics[{Dynamic[EdgeForm[{Opacity[opacity], color}]],
           Dynamic[color], Dynamic[Opacity[opacity]], primitives},
          ImageSize -> {{400, Large}, {400, Large}},
          Background -> Dynamic[background]]],
      
        (*declare variables here for persistence*)
        {{background, bg = White}, ControlType -> None},
        {{matrix1, {{1, 0}, {1, 1}}}, ControlType -> None},
        {{matrix, {{1, 0}, {1, 1}}}, ControlType -> None},
      
        Evaluate[controls],
        Bookmarks :> Evaluate[bookmarks],
      
        LabelStyle -> Darker[Gray], SynchronousUpdating -> Automatic,
        Paneled -> False, SaveDefinitions -> True, Alignment -> Center]
      

Cha-ching baby. If Snoop Dogg ever used Mathematica, that's what square brackets in his custom font would look like. And I know what this one reminds you of. The folds of the brain. And check out the Black Riddler's Question Mark.


Inversion

What happens if we turn the Sierpinski triangle "inside out"? This is easy to answer because all we have to do is invert. You may be familiar with this plot of $\sin \small{\frac 1 x}$:



It's typically given as an example of a function that isn't differentiable at a point (at 0 in this case). It can be seen as a composition of two functions:

$$\frac 1 x \longrightarrow \sin x$$

The important function is the first one. It inverts the entire number line around 1, mapping $[1,\infty)$ to $[1,0)$ and vice versa. The reason the plot of $\sin \small{\frac 1 x}$ looks like that is because it's essentially the regular sine function with its values from 1 to infinity all crammed between 0 and 1. In this sense, $\small{\frac 1 x \rightarrow f}$ is like a logarithmic plot on hypercrack for $\small f$.

So to turn our Sierpinski triangle inside out, we can do the same thing. For each point, we invert its distance but keep it at the same angle, using $\small{\frac z {|z|}}$ to normalize the point:

$$z_\text{inv} = \frac 1 {|z|}\frac z {|z|} = \frac z {|z|^2}$$
  1. axiom = Polygon[{Cos[#], Sin[#]} & /@ (Pi/2 + 2 Pi Range[3]/3)];
    
    next[prev_] := prev /. Polygon[{p1_, p2_, p3_}] :> {
         Polygon[{p1, (p1 + p2)/2, (p1 + p3)/2}],
         Polygon[{p2, (p2 + p3)/2, (p1 + p2)/2}],
         Polygon[{p3, (p1 + p3)/2, (p2 + p3)/2}]};
    
    invert[p_] := p/Norm[p]^2;
    
    draw[n_] := Graphics[{EdgeForm[Black], Nest[next, N@axiom, n]}];
    
    g = draw[2];
    Show[g, g /. Polygon[pts_] :> Polygon[invert /@ pts]]
    

The radius of inversion is right at the corners of the triangle, and I've left the univerted triangle in the center. Here's what the first few construction steps of the triangle look like if we invert them:

  1. axiom = Polygon[{Cos[#], Sin[#]} & /@ (Pi/2 + 2 Pi Range[3]/3)];
    
    next[prev_] := prev /. Polygon[{p1_, p2_, p3_}] :> {
         Polygon[{p1, (p1 + p2)/2, (p1 + p3)/2}],
         Polygon[{p2, (p2 + p3)/2, (p1 + p2)/2}],
         Polygon[{p3, (p1 + p3)/2, (p2 + p3)/2}]};
    
    invert[p_] := p/Norm[p]^2;
    
    draw[n_] := Module[{ps = Nest[next, N@axiom, n]},
       Graphics[{EdgeForm[Black], Transparent,
         ps, ps /. Polygon[pts_] :> Polygon[invert /@ pts]}]];
    
    Grid[Partition[draw /@ Range[0, 8], 3]]
    

Notice that we're just inverting the endpoints of the lines, not the lines-as-curves. Visually this doesn't make a difference at higher iterations:

  1. axiom = Polygon[{Cos[#], Sin[#]} & /@ (Pi/2 + 2 Pi Range[3]/3)];
    
    next[prev_] := prev /. Polygon[{p1_, p2_, p3_}] :> {
         Polygon[{p1, (p1 + p2)/2, (p1 + p3)/2}],
         Polygon[{p2, (p2 + p3)/2, (p1 + p2)/2}],
         Polygon[{p3, (p1 + p3)/2, (p2 + p3)/2}]};
    
    invert[p_] := p/Norm[p]^2;
    
    draw[n_] := Graphics[{
        EdgeForm[Black], Transparent,
        Nest[next, N@axiom, n]}];
    
    Show[draw[10], draw[12] /. Polygon[pts_] :> Polygon[invert /@ pts],
       Method -> {"ShrinkWrap" -> True}, ImageSize -> 4 750] //
      Rasterize // ImageResize[#, Scaled[1/4]] &
    

What about varying the radius of inversion? You first perform the same inversion as before, but with respect to the radius:

$$\frac r {|z|}$$

Normally this would be enough to get the inverted distance, but division is performed through the lens of the unit $1$. When we perform that $\small{\frac r {|z|}}$, we "lose" the information about what the radius may have been. (This innocent-sounding thing strikes me as much more involved than it appears). If you work out what an inversion should look like on the number line you'll find that you have to scale the result back by multiplying by the radius:

$$\frac r {|z|} r = \frac {r^2} {|z|}$$

It took me a while, but eventually I realized that the edges of the triangle were being mapped to curves, and that if you continued those curves they would form circles that intersected the origin, like this:

What this means is that by this inversion, infinite lines and circles that cross the origin are inverses of eachother. This realization almost punched my brain in its face, but apparently this is well-known. In fact it's called inversive geometry. I found myself quite disappoint, however, that online descriptions presented our $r^2$ factor as part of the definition, rather than as the arithmetic consequence of the inversion operation. Son.

Let's not forget we have a bountiful cornucopia from which to invert:

    1. draw[v_, n_] := Module[{ring},
         ring[c_, r_, depth_] := Module[{ps},
           ps = c + r {Cos[#], Sin[#]} & /@ (2 Pi Range[v]/v);
      
           If[depth == 0, Polygon[ps],
            ring[(c + #)/2, r/2, depth - 1] & /@ ps]];
      
         Graphics[{Transparent,
           EdgeForm[Black],
           ring[0., 1., n]}]];
      
      Clear[invert];
      (*invert[p_/;Norm[p]<.0001]:={0,0};*)
      invert[p_] := p/Norm[p]^2;
      
      Column[Panel[Row[#]] & /@
        Table[
         draw[v, n] /. Polygon[pts_] :> Polygon[invert /@ pts],
         {v, 3, 6}, {n, 0, 4}]]
                  
    2. draw[v_, n_] := Module[{ring},
         ring[c_, r_, depth_] := Module[{ps},
           ps = c + r {Cos[#], Sin[#]} & /@ (2 Pi Range[v]/v);
      
           If[depth == 0, Polygon[ps],
            ring[(c + #)/2, r/2, depth - 1] & /@ ps]];
      
         Graphics[ring[0., 1., n]]];
      
      invert[p_] := p/Norm[p]^2;
      
      Column[Panel[Row[#]] & /@
        Table[
         draw[v, n]
           /. Polygon[pts_] :>
            Line /@ Quiet@Partition[invert /@ pts, 2, 1, 1]
          /. l_Line /; MemberQ[l, Indeterminate, Infinity] :> {},
         {v, 3, 6}, {n, 0, 4}]]
                  

Those light red shades are Mathematica QQing about plotting points at infinity. I thought Mathematica can do anything??? The problem is that the inverse of the origin under this scheme is essentially "everything at infinity" (from division by 0) and algebraically this inverse doesn't even have any specific 'direction' like $\pm\space\infty$ do. The easiest solution is to just leave points at (0, 0) alone or remove the incident lines.

Sidenote. Notice that in this program we aren't even touching our original uninverted geometric renderer, because we don't need to. Our original renderer returns a Graphics structure. This structure (which you might call an M-expression) is to us a set of straightforward vector graphics directives, but is to Mathematica meaningless until the frontend gets ahold of it. Until then (and even afterwards) we can perform the same kinds of structural slicing and dicing that we can perform on any other structure. In this case, replacing points by their inverses.

A more complete solution to our point at infinity/division by 0 problem is to put the inverse of (0, 0) not at infinity, but really far. This doesn't come out of the algebra, but we can do it in a well-behaved way because we know from which direction our lines are coming from since we're defining things as polygons:

  1. invert[p_] := p/Norm[p]^2;
    
    (*you ever see that show Long Ago and Far Away? that show was awesome*)
    invertPoly[Polygon[pts_], farAway_: 20000] :=
      With[{indQ = MemberQ[#, Indeterminate] &},
       Line /@ Quiet@Partition[invert /@ pts, 2, 1, 1] /.
        {_?indQ, p_} | {p_, _?indQ} :> {p, farAway Normalize[p]}];
    
    plotRangeInv[g_Graphics] := PlotRange /.
       AbsoluteOptions[g /. Polygon[pts_] :> Quiet@Polygon[invert /@ pts]];
    
    draw[v_, n_] := Module[{ring},
       ring[c_, r_, depth_] := Module[{ps},
         ps = c + r {Cos[#], Sin[#]} & /@ (2 Pi Range[v]/v);
    
         If[depth == 0, Polygon[ps],
          ring[(c + #)/2, r/2, depth - 1] & /@ ps]];
    
       Graphics[ring[0., 1., n]]];
    
    (*if it ever came out on DVD i'd buy it like 100 times*)
    drawInv[v_, n_] := Module[{g = draw[v, n]},
       Show[g /. poly_Polygon :> invertPoly[poly],
        PlotRange -> 1.1 plotRangeInv[g]]];
    
    (*lines=Cases[drawInv[6,4],Line[ps_]/;
    EuclideanDistance@@ps<10000:>Line[Sort[ps]],Infinity];
    Graphics[DeleteDuplicates@lines]*)
    

What if, for no particular reason, we vary the exponents of the inversion formula?

$$\frac {z^{\lt {\text{insert number here}} \gt} \leftarrow {\small\text{(elementwise exponent)}}} {|z|^{\lt {\text{insert some other, not necessarily distinct, number here} \gt}}}$$

Most of the results of this were boring, but the one for $\small z^3 / {|z|^2}$ was cool:

  1. axiom = Polygon[{Cos[#], Sin[#]} & /@ (Pi/2 + 2 Pi Range[3]/3)];
    
    next[prev_] := prev /. Polygon[{p1_, p2_, p3_}] :> {
         Polygon[{p1, (p1 + p2)/2, (p1 + p3)/2}],
         Polygon[{p2, (p2 + p3)/2, (p1 + p2)/2}],
         Polygon[{p3, (p1 + p3)/2, (p2 + p3)/2}]};
    
    invert[p_] := p^3/Norm[p]^2;
    
    draw[n_] := Graphics[{EdgeForm[Black],
        Nest[next, N@axiom, n]}];
    
    Grid[Partition[draw /@ Range[0, 8], 3]] /.
     Polygon[pts_] :> Polygon[invert /@ pts]
    

See this one. One day you're going to be driving home from work. It's going to be dark. Pitch black. All a sudden out the corner your eye you're gonna see a flash in your rear view mirror. And when you look, you're gonna see that same Black Cobra Grill on my car speeding towards you at some unspeakable number of kilometers per hour. And then I'll disappear into the night. Like an episode off an MJ's Thriller$\small\times$Knight Rider mashup.

If we mangle the formula every which way we can find a lot of interesting effects:

  1. axiom = Polygon[{Cos[#], Sin[#]} & /@ (Pi/4 + 2 Pi Range[3]/3)];
    
    next[prev_] := prev /. Polygon[{p1_, p2_, p3_}] :> {
         Polygon[{p1, (p1 + p2)/2, (p1 + p3)/2}],
         Polygon[{p2, (p2 + p3)/2, (p1 + p2)/2}],
         Polygon[{p3, (p1 + p3)/2, (p2 + p3)/2}]};
    
    transform[p_] := (Reverse[p].p) p^2/Norm[p]^2;
    
    drawFishie1[n_] := Graphics[{Red,
        EdgeForm[{Thickness[.01], Opacity[.3],
          JoinForm["Round"], Lighter[Blue, .6]}],
        Rotate[Nest[next, N@axiom, n] /.
          Polygon[pts_] :> Polygon[transform /@ pts], 3 Pi/4]},
       PlotRange -> .8 {{-.85, 1.51}, .4 {-1, 1.1}}];
    
    drawFishie2[n_] := Graphics[{
        Transparent, EdgeForm[Black],
        Rotate[Nest[next, N@axiom, n] /.
          Polygon[pts_] :> Polygon[transform /@ pts], -Pi/4]}];
    

The self-crossings form hexagonal figures. And American iconography? Here's another nifty one:

  1. next[prev_] := prev /. Polygon[{p1_, p2_, p3_}] :> {
         Polygon[{p1, (p1 + p2)/2, (p1 + p3)/2}],
         Polygon[{p2, (p2 + p3)/2, (p1 + p2)/2}],
         Polygon[{p3, (p1 + p3)/2, (p2 + p3)/2}]};
    
    transform[p_] := (Reverse[p].p) p/Norm[p]^2;
    
    notFishieAwwwSadFace[n_, a1_: 0, a2_: 0, options___] :=
      Module[{axiom},
       axiom = Polygon[{Cos[#], Sin[#]} & /@ (a1 + 2 Pi Range[3]/3)];
    
       Graphics[{Transparent, EdgeForm[Black],
         Rotate[Nest[next, N@axiom, n] /.
           Polygon[pts_] :> Polygon[transform /@ pts], a2]},
        options]];
    
    notFishieAwwwSadFace[6, Pi/4]
    

I call it the Sierpinski Stiletto of Triangular Destruction. Hell yea. Also pay heed to the Sierpinski Butterfly of Poisonous Death, lest yee regret it. We can also move the circle of inversion around. I was going to write a program to do only that, but before I realized it I had accidentally built this:

    1. polys[v_, n_, offset_: {0, 0}, size_: 1, rot_: 0] := Module[{ring},
         ring[c_, r_, depth_] := Module[{ps},
           ps = c + r {Cos[#], Sin[#]} & /@ (rot + 2 Pi Range[v]/v);
      
           If[depth == 0, Polygon[(offset + # &) /@ ps],
            ring[(c + #)/2, r/2, depth - 1] & /@ ps]];
      
         ring[0., size, n]];
      
      (*Polygon in, Lines out /statictyping*)
      invertPoly[Polygon[pts_], tf_, farAway_: 20000] :=
        With[{indQ = MemberQ[#, Indeterminate] &},
         Line /@ Quiet@Partition[tf /@ pts, 2, 1, 1] /.
          {_?indQ, p_} | {p_, _?indQ} :> {p, farAway Normalize[p]}];
      
      d = Norm;
      (*d=ChessboardDistance[{0, 0},#]&;*)
      transformationFunctions = {
         #/d[#]^2 &,
         #^3/d[#]^2 &,
         (Reverse[#].#) #/d[#] &,
         (# + #^3/d[#]^2)/2 &,
         (# + #/d[#]^2)/2 &,
         Round[#, .05] &};
      
    2. With[{
        verticesC = Control[{{vertices, 3}, 3, 8, 1, Appearance -> "Labeled", ImageSize -> Small}],
        iterationsC = Control[{{iterations, 5}, 0, 10, 1, Appearance -> "Labeled", ImageSize -> Small}],
        rangeC = Control[{{range, 4}, 1, 12, Appearance -> "Labeled"}],
        rotC = Control[{{rot, -Pi/6}, -Pi, Pi, Appearance -> "Labeled", ImageSize -> Tiny}],
        functionC = {{tf, transformationFunctions[[1]], "function"},
          (# -> TraditionalForm@Quiet@Trace[#[z]][[2]] &) /@ transformationFunctions,
             ControlType -> SetterBar},
        originalC = Control[{original, {True, False}}],
        circleC = Control[{circle, {True, False}}],
        opacityC = Control[{{opacity, 1}, 0, 1, ImageSize -> Small}],
        sizeC = Control[{{size, 1}, .001, 4, ImageSize -> Small}],
        offsetC = {{offset, {0, 0}}, Locator}},
      
       Manipulate[Module[{polygons},
         polygons = polys[vertices,
           ControlActive[Min[3, iterations - 2], iterations],
           offset, size, rot];
      
         Graphics[{
           If[circle, {LightGray, Circle[]}],
           If[original, {EdgeForm[LightGray], {Transparent, polygons}}],
           {Black, polygons /. p_Polygon :> invertPoly[p, tf]}},
          PlotRange -> range]],
      
        Row[{verticesC, iterationsC}],
        Row[{rangeC, rotC}], functionC,
        Row[{originalC, circleC, sizeC, opacityC}, "  "], offsetC,
        SynchronousUpdating -> False]]
      

Oops. This hilarious function doesn't allow anything inside the unit disk. It's just waiting for someone to make a Yakety Sax movie about shapes crashing into the circle and crawling around it.

Someone on the internet asked an interesting question: Are there "zoom out" fractals? We know that if we zoom in on the Sierpinski triangle, we'll continue seeing detail endlessly. But are there fractals that no matter how far you zoom out, you can't get out of them?

Of course there are. We can just take a quote-unquote "zoom in" fractal and place one of its points of detail right at the origin, and then invert the fractal. Because the inverse of the origin is some kind of crazy infinity, we know that no matter how far we zoom out, we won't reach the end of the fractal. This example is really a formality though. You have a lot of liberty to make things up in math.

Cornucopia.

    1. invert[p_ /; Norm[p] < .0001] := {0, 0};
      invert[p_] := p/Norm[p]^2;
      
      draw[v_, df_, n_] := Module[{ring},
         ring[c_, r_, depth_] := Module[{ps},
           ps = c + r {Cos[#], Sin[#]} & /@ (2. Pi Range[v]/v);
      
           If[depth == 0, Polygon[invert /@ ps],
            ring[(c + #)/2, df[0, r], depth - 1] & /@ ps]];
      
         Graphics[{EdgeForm[White], Opacity[.4],
           RGBColor[.4, 1(*;.6*), 1], ring[0., 1., n]}]];
      
      candy[g_, res_: 600, upsc_: 1, style_: EdgeForm[Thick]] :=
        Module[{a = Show[g, ImageSize -> upsc res, Background -> White]},
      
         a = a /. p_Polygon :>
            {RGBColor[.6, RandomReal[], RandomReal[]], style, p};
         a = Rasterize[a];
      
         (*move this downscale to end for different quality*)
         a = ImageResize[a, Scaled[1/upsc]];
         ImageDifference[a, ImageReflect[a, Left]] // ColorNegate];
      
      g = With[{f = (#1 + #2)/RandomChoice[Prime[Range[4]]] &},
        Show[
         Table[draw[(*repeatedly draw to cover more possibilities*)
            RandomChoice[{1, 1, .25} -> {3, 4, 5}], f,
            RandomChoice[{.1, 1.5} -> {2, 3}]]
           /. p_Polygon :> Rotate[p, 0(*;Pi/4*), {0, 0}],
          {12}], Background -> Black, ImageSize -> Medium]];
      
      (*note you can edit g in-place*)
      Defer[candy][g, 1280, 4]
      
    1. ring[c_, r_, depth_] := Module[{ps},
         ps = c + r {Cos[#], Sin[#]} & /@ (-Pi/10 + 2. Pi Range[5]/5);
      
         If[depth == 0, Polygon[ps],
          ring[#, r/2, depth - 1] & /@ ps]];
      
      invert[p_] := p/Norm[p]^2;
      
      Graphics[{Opacity[.15], Black,
         ring[0., 1., #] & /@ Range[5(*;8*)]}] /.
       Polygon[pts_] :> Polygon[invert /@ pts, VertexColors ->
         (ColorData["AvocadoColors"] /@ (#^1.7 &) /@ Norm /@ pts)]
      

From what I can tell, one of the settings used to deal with division by 0 is the so-called Riemann sphere, which is where we take a space shuttle and use it to fly over and drop a cow on top of a biodome, and then have the cow indiscriminately fire laser beams at the grass inside and around the biodome. That's my intuitive understanding of it anyway.

  1. shuttle = ExampleData[{"Geometry3D", "SpaceShuttle"}, "GraphicsComplex"];
    cow = ExampleData[{"Geometry3D", "Cow"}, "GraphicsComplex"];
    grass = ExampleData[{"Texture", "Grass4"}];
    drop = Sound[Play[Sin[1000 (2 - t)^2.3], {t, 0, .9}]];
    lazer[] := EmitSound[Sound[Play[{TriangleWave[10 (2 - t)^5],
          Sin[500 (2 - t)^5] + SawtoothWave[20 (2 - t)^4.99]},
         {t, 0, .2}, SampleRate -> 4000]]];
    
    Animate[Module[{dir, angle},
      (*xs animates from +30 to -60*)
      If[xs > 25, dropTrigger = False];
      If[dropTrigger == False && xs < 10, dropTrigger = True; EmitSound[drop]];
      If[xs < -10,
       lazer[];
       dir = RandomReal[{-15, 15}, 2]~Join~{-20};
       angle = ArcTan @@ Take[dir, 2], angle = 0];
    
      Labeled[#, Style["Understanding the Riemann Sphere", FontFamily -> "Verdana"], Top] &@
       With[{
         cowLoc = {Clip[xs, {1.3, 50}], 0, Clip[8 + (-(.25 xs - 5)^2), {-13, -.25}]},
         greenLight = Lighting -> {{"Directional", Green, {{5, 5, 5}, {0, 0, 0}}}}},
    
        Graphics3D[{EdgeForm[None],
          (*shuttle*){Specularity[White, 7], Translate[shuttle, {xs, 0, 0}]},
          (*cow*) Rotate[Translate[Scale[cow, 5], cowLoc], angle, {0, 0, 1}],
    
          (*grass*) Translate[{If[False, Sequence @@ {greenLight, Texture[grass]}, Green],
            Polygon[50 {{-1, -1, 0}, {1, -1, 0}, {1, 1, 0}, {-1, 1, 0}},
             VertexTextureCoordinates -> 1.1 {{0, 0, 0}, {1, 0, 0}, {1, 1, 0}, {0, 1, 0}}]},
            {0, 0, -19}],
    
          (*biodome*){Specularity[White, 3], Opacity[.5], Sphere[.95 {0, 0, -20}, 5]},
          (*sun*){Glow[White], Sphere[{18, 18, 18}, 2]},
          (*laser*) If[xs < -10, {Red, Glow[Red], Opacity[.5], Tube[{{0, 0, -12.4}, dir}]}]},
         (*sky*)Background -> LightBlue, PlotRange -> 20, Boxed -> False]]],
     {xs, 30, -60, 3.2}, AnimationRate -> 10, DisplayAllSteps -> True]
    

(Note the cow cannot be spherical or it will roll off). Personally I don't have any beef with Riemann or any of his manifolds, but for our purposes the Riemann sphere is inadequate since it maps our inverses vertically. One interesting consequence of this is that in the 2D cross section where the imaginary component is zero (essentially the 'Weierstrass circle'), it maps multiplicative inverses vertically and additive inverses horizontally. This all seems mathematically expedient, but it's otherwise boring.

The Riemann sphere does give one explanation though about 'why' our circles and lines are inverses. In the Riemann sphere, the inverse of a circle that crosses the origin is a circle that crosses the North Pole, and since the lasers are being shot from the North Pole, they're limited to tracing out a line as they follow the circle. I was going to make a simple 3D diagram demonstrating this, but I accidentally made this:

    1. toRiemann = Compile[{{pts, _Real, 2}}, Module[{k},
          Map[(k = 2/(1 + #.#);
             {k #[[1]], k #[[2]], 1 - k}) &, pts]]];
      
      fromRiemann[pts_] := (-1/(#3 - 1)) {#1, #2, 0} & @@@ pts;
      
      shuttle = ExampleData[{"Geometry3D", "SpaceShuttle"}, "GraphicsComplex"];
      cow = Append[ExampleData[{"Geometry3D", "Cow"}, "GraphicsComplex"],
         VertexTextureCoordinates -> 1/500 ExampleData[{"Geometry3D", "Cow"}, "PolygonData"]];
      sphere = SphericalPlot3D[1, {\[Theta], 0, Pi}, {\[Phi], 0, 2 Pi},
         MeshStyle -> Opacity[.05], PlotStyle -> Opacity[.1]];
      plane = {Lighting -> "Neutral", Opacity[.5], LightGray, EdgeForm[None],
         Line[{{-50, 0, 0}, {50, 0, 0}}], Line[{{0, -50, 0}, {0, 50, 0}}],
         Polygon[15 {{-1, -1, 0}, {1, -1, 0}, {1, 1, 0}, {-1, 1, 0}}]};
      
      numPoints = 160/4;
      circle = .25 {1 + Cos[#], Sin[#], 0} & /@ (-Pi + 2 Pi Range[numPoints]/numPoints);
      rCircle = toRiemann[circle];
      invRCircle = {#1, #2, -#3} & @@@ rCircle;
      invCircle = Quiet[fromRiemann[invRCircle]];
      
      ind = Indeterminate {1, 1, 1};
      transform[\[Theta]_] := Composition[
         RotationTransform[\[Theta], {0, 0, 1}],
         TranslationTransform[{0, 0, 1.2}]];
      {l, r} = {{.335, -.044, .1894}, {.335, .044, .1894}};
      
    2. slides = ParallelTable[Module[{angle},
          angle = Quiet[ArcTan @@ (invCircle /. ind -> {0, 0})[[pti, 1 ;; 2]] /.  ArcTan[0, 0] -> 0];
          Show[sphere,
           Graphics3D[{
             plane, Opacity[.7], Sphere[{0, 0, 1}, .01],
             (*cow*) Rotate[#, angle, {0, 0, 1}] &@Translate[#, {0, 0, 1.2}] &@
              {EdgeForm[None],(*Opacity[.999],*)Texture[Graphics[Disk[]]],
               Lighting -> "Neutral", Lighting -> {{"Point", Darker[Red], l}},
               cow, Red, Glow[Red], Sphere[{l, r}, .01]},
             (* keep shuttle in orbit in case need more cows *)
             {EdgeForm[None], Translate[shuttle, {0, 0, 100}]},
             (*lazerz*) If[invCircle[[pti]] =!= ind, {Red,
               Line[{transform[angle][l], invCircle[[pti]]}],
               Line[{transform[angle][r], invCircle[[pti]]}],
               {Red, Glow[Red], Opacity[.1], Sphere[invCircle[[pti]], .02 RandomReal[]]}}],
             (*etc*) {Lighter[Gray], Dashed, Line[{rCircle[[pti]], invRCircle[[pti]]}]},
             {Opacity[.1], Lighter[Blue], Line[circle]},
             {Opacity[.5], Lighter[Blue], Line[Take[circle, pti]]},
             {Lighter[Blue], Line[{{0, 0, 1}, 50 (-{0, 0, 1} + circle[[pti]])}]},
             {Lighter[Blue], Line[Take[rCircle, pti]]},
             {Red, Line[Take[invRCircle, pti]]},
             (*Purple,Line[{{0,0,1},50(-{0,0,1}+invCircle[[pti]])}/.tride->{0,0,1}],*)
             (*burn mark*) Thick, Red, Line[DeleteCases[Take[invCircle, pti], ind]]}],
           ImageSize -> 1/4 {16, 9} (1080/9),
           ViewAngle -> 4 Degree, PlotRange -> 10, Boxed -> False, Axes -> False]],
         {pti, 1, Length[circle]}];
      ListAnimate[slides]
      

Oops. But since we now have this tool, let's see what other plots look like:

  1. toRiemann = Compile[{{pts, _Real, 2}}, Module[{k},
        Map[(k = 2/(1 + #.#);
           {k #[[1]], k #[[2]], 1 - k}) &, pts]]];
    
    shuttle = ExampleData[{"Geometry3D", "SpaceShuttle"}, "GraphicsComplex"];
    cow = Append[ExampleData[{"Geometry3D", "Cow"}, "GraphicsComplex"],
       VertexTextureCoordinates -> 1/500 ExampleData[{"Geometry3D", "Cow"}, "PolygonData"]];
    sphere = SphericalPlot3D[1, {\[Theta], 0, Pi}, {\[Phi], 0, 2 Pi},
       MeshStyle -> Opacity[.05], PlotStyle -> Opacity[.1]];
    plane = {Lighting -> "Neutral", Opacity[.5], LightGray, EdgeForm[None],
       Line[{{-50, 0, 0}, {50, 0, 0}}], Line[{{0, -50, 0}, {0, 50, 0}}],
       Polygon[15 {{-1, -1, 0}, {1, -1, 0}, {1, 1, 0}, {-1, 1, 0}}]};
    
    listRiemannPlot[pts_] := Module[{rPts = toRiemann[pts]},
       Show[
        sphere,
        Graphics3D[{
          plane, Opacity[.7], Sphere[{0, 0, 1}, .01],
    
          (*cow*)Translate[#, {0, 0, 1.2}] &@
           {EdgeForm[None],(*Opacity[.999],*)
            Texture[Graphics[Disk[]]], Lighting -> "Neutral", cow},
    
          (*shuttle*){EdgeForm[None], Translate[shuttle, {0, 0, 100}]},
          (*original*){Opacity[2 .2], Lighter[Blue], Line[pts]},
          (*riemannized*){Opacity[.8], Blue, Line[rPts]}}],
    
        ViewAngle -> 4 Degree, PlotRange -> 10, Boxed -> False, Axes -> False]];
    
    listRiemannPlot[Table[{x, Sin[2 x], 0}, {x, -40 Pi, 40 Pi, .01}]]
    

I call this one the Riemann cowsine, sort of like cosine but with "cow" instead of "co". This one is $\small{\frac 1 {2 \sin(2x)}}$. If you've seen the tangent function, you know it has a lot of infinities, which means the cowtangent is going to have a lot of circles. I suppose the fact that the circles enter and exit the north pole without being deflected is a result of the asymptotic behavior of the function being the same going up as going down. Come to think of it, those circles are more like a continous infinity symbol that goes in infinitely on itself. And look:

  1. toRiemann = Compile[{{pts, _Real, 2}}, Module[{k},
        Map[(k = 2/(1 + #.#);
           {k #[[1]], k #[[2]], 1 - k}) &, pts]]];
    
    shuttle = ExampleData[{"Geometry3D", "SpaceShuttle"}, "GraphicsComplex"];
    cow = Append[ExampleData[{"Geometry3D", "Cow"}, "GraphicsComplex"],
       VertexTextureCoordinates -> 1/500 ExampleData[{"Geometry3D", "Cow"}, "PolygonData"]];
    sphere = SphericalPlot3D[1, {\[Theta], 0, Pi}, {\[Phi], 0, 2 Pi},
       MeshStyle -> Opacity[.05], PlotStyle -> Opacity[.1]];
    plane = {Lighting -> "Neutral", Opacity[.5], LightGray, EdgeForm[None],
       Line[{{-50, 0, 0}, {50, 0, 0}}], Line[{{0, -50, 0}, {0, 50, 0}}],
       Polygon[15 {{-1, -1, 0}, {1, -1, 0}, {1, 1, 0}, {-1, 1, 0}}]};
    
    riemannize[Graphics[g_, rest___], options___] :=
      Show[sphere,
       Graphics3D[{plane,
         (*cow*)Translate[#, {0, 0, 1.2}] &@
          {EdgeForm[None],
           Texture[Graphics[Disk[]]], Lighting -> "Neutral", cow},
         (*shuttle*)Translate[shuttle, {0, 0, 100}],
         (*curves*)g /. Line[pts_] :> {Line[toRiemann[pts]],
            Opacity[.25], Line[{#1, #2, 0} & @@@ pts]}}],
       Boxed -> False, Axes -> None, PlotRange -> 5, options];
    
    plot = Plot[
       Evaluate[y /. Solve[y^2 == x^3 - 3 x + 1, y]],
       {x, -100, 100}, PlotPoints -> 10000];
    
    riemannize[plot, ViewAngle -> 15 Degree]
    

This is one of those fangled elliptic curves. Apparently they do form pairs of circle things on the Riemann sphere. I thought that was just an old wive's tale. The nice thing about this program is that it works on Graphics structures, such as those returned by Plot. That means you can plug arbitrary 2D plots and graphics into this function and have them automatically Riemannized. Like say you're trying to educe from without your incorrigible students' crania some particular factoid:

  1. toRiemann = Compile[{{pts, _Real, 2}}, Module[{k},
        Map[(k = 2/(1 + #.#); {k #[[1]], k #[[2]], 1 - k}) &, pts]]];
    
    shuttle = ExampleData[{"Geometry3D", "SpaceShuttle"}, "GraphicsComplex"];
    cow = Append[ExampleData[{"Geometry3D", "Cow"}, "GraphicsComplex"],
       VertexTextureCoordinates -> 1/500 ExampleData[{"Geometry3D", "Cow"}, "PolygonData"]];
    sphere = SphericalPlot3D[1, {\[Theta], 0, Pi}, {\[Phi], 0, 2 Pi},
       MeshStyle -> Opacity[.05], PlotStyle -> Opacity[.1]];
    plane = {Lighting -> "Neutral", Opacity[.5], LightGray, EdgeForm[None],
       Line[{{-50, 0, 0}, {50, 0, 0}}], Line[{{0, -50, 0}, {0, 50, 0}}],
       Polygon[15 {{-1, -1, 0}, {1, -1, 0}, {1, 1, 0}, {-1, 1, 0}}]};
    
    riemannize[Graphics[g_, ___], options___] := Show[sphere,
       Graphics3D[{plane,
         (*cow*)Translate[#, {0, 0, 1.2}] &@
          {EdgeForm[None],(*Opacity[.999],*)
           Texture[Graphics[Disk[]]], Lighting -> "Neutral", cow},
         (*shuttle*)Translate[shuttle, {0, 0, 100}],
         (*curves*)g /. Line[pts_] :> {Line[toRiemann[pts]],
            Opacity[.25], Line[{#1, #2, 0} & @@@ pts]}}],
       Boxed -> False, Axes -> None, PlotRange -> 5, options];
    
    vint = Sqrt[1 - 3/2 (1 + Sqrt[21]) + 1/8 (1 + Sqrt[21])^3];
    plot = Show[Plot[y /. Solve[y^2 == x^3 - 3 x + 1, y], {x, -100, 100}, PlotPoints -> 10000],
       Graphics[{Orange, Line[Table[{x, x + 1}, {x, -500, 500, .1}]],
         Black, Dashing[1/10 {0.08, 0.02}],
         Line[Table[{1/2 (1 + Sqrt[21]), y}, {y, -vint, vint, .1}]]}]];
    
    riemannize[plot, ViewAngle -> 8 Degree]
    

You've set this up with 2D graphics. But you can just plug the output of this into our Riemannizer to get this. In fact in Mathematica you can even copy/paste the 2D plot (itself an interactwithable vector object) like this. And you can vector-edit that plot in-place and when you re-evaluate the expression, the differences will appear in the Riemannization. Not bad for what essentially amounts to one line of code:

g /. Line[pts_] :> Line[toRiemann[pts]]

This is the power of Mathematica's macro-at-will symbolic semantics and well-curated architecture. Specifically in this case, it's the fact that the built-in plotting functions return the same laid-bare Graphics vector structures that your own versions of those functions would return. This Riemannizer only does a direct endpoint conversion of lines, but you can easily have it 3Dify whatever you want in a more thorough fashion.

After I made my Cyclotron 4000 masterpiece, I considered what a version 2 might be. Now I know. With some adjustments to the contraption, we now have the Cycowtron 4800 Deluxe (pronounced psy-cow-tron forty-eight-hundred de-lux):

    1. toRiemann = Compile[{{pts, _Real, 2}}, Module[{k},
          Map[(k = 2/(1 + #.#); {k #[[1]], k #[[2]], 1 - k}) &, pts]]];
      
      invertRiemann[pts_] := {#1, #2, -#3} & @@@ pts;
      (*fromRiemann[pts_]:=(-1/(#3-1)) {#1,#2,0}&@@@pts;*)
      
      shuttle = With[{shuttleGC = ExampleData[{"Geometry3D", "SpaceShuttle"}, "GraphicsComplex"]},
         Translate[#, {0, 0, 100}] &@
          {EdgeForm[None], Lighting -> "Neutral",
           Append[shuttleGC, VertexColors -> RandomReal[.65 + {0, 1}, Length[shuttleGC[[1]]]]^2]}];
      
      cow = Translate[#, {0, 0, 1.2}] &@
         {EdgeForm[None],(*Opacity[.999],*)Texture[Graphics[Disk[]]], Lighting -> "Neutral",
          Append[ExampleData[{"Geometry3D", "Cow"}, "GraphicsComplex"],
           VertexTextureCoordinates -> 1/500 ExampleData[{"Geometry3D", "Cow"}, "PolygonData"]]};
      
      allSettings = {"sphere", "cow", "original", "riemann", "inverse", "shuttle"};
      
      bookmarks = {
         "Random" :> (
           With[{R := RandomReal[]},
            color = RGBColor[R, R, R];
            inflection = RandomChoice[{1, -1}];
            effect = .25*R^8; thickness = 10*R^8;
            spirality = R^3;
            If[spirality < .035, spirality = 0];
            x = 12 R; y = 2 R];
      
           Module[{weights = If[x > 3, {.1, .03, .87}, {.1, .07, .83}]},
            \[Zeta] = RandomChoice[{-1, 1}]*
              RandomChoice[weights -> {
                 RandomReal[{0, 7.5}],
                 RandomChoice[{0, .5, .5}],
                 Round[RandomReal[{1, 7.5}], .5]}]]),
      
         "Rose" :> {color = Red, effect = 0, inflection = -1, x = 0, y = 1.1343, spirality = 0, \[Zeta] = 3},
         "Glyph" :> {color = Black, effect = 0.198, inflection = 1, x = 5.2, y = 0, spirality = 0, \[Zeta] = 2},
         "Mass Atomic" :> { effect = 0, inflection = -1, x = 5.84, y = 0.412, spirality = 0, \[Zeta] = -4.2504},
         "Jello" :> {color = Red, effect = 0, inflection = -1, x = 12, y = 0.846, spirality = 1, \[Zeta] = -1},
         "Grim" :> {thickness = 3.35, color = Black, effect = 0.0675, inflection = 1, x = 8, y = 0.296, spirality = 1, \[Zeta] = -1},
         "Angelwings" :> {color = RGBColor[.07694, .39046, 1], effect = 0, inflection = 1, spirality = 1, x = 12, y = 0, \[Zeta] = -5.4947},
         "Rollers" :> {color = Black, effect = 0, inflection = -1, spiral3ity = 0, x = 5.5, y = .1, \[Zeta] = -0.984032039033508},
         "Lifespark" :> {color = RGBColor[.1026, .9878, .0201], effect = 0, inflection = 1, spirality = .0995, x = 3.2757, y = .2002, \[Zeta] = -5.5}};       
      
    2. With[{
        colorC = Control[{{color, Black, "line color"}, ColorSlider}],
        backgroundC = Control[{{background, White}, ColorSlider}],
        thicknessC = Control[{{thickness, .001, "line thickness"}, .001, 10, Appearance -> "Labeled", ImageSize -> Small}],
        effectC = Control[{{effect, 0., "charcoal effect"}, 0, .25, Appearance -> "Labeled", ImageSize -> Small}],
        inflectionC = Control[{{inflection, 1}, {1 -> " concave ", -1 -> " convex "}, Appearance -> "Vertical"}],
        angularityC = Control[{{\[Zeta], 2., "angularity"}, -7.5, 7.5, .5, Appearance -> "Labeled", ImageSize -> Small}],
        tensionC = Control[{{x, 8., "tension"}, 0, 12, Appearance -> "Labeled", ImageSize -> Small}],
        yC = Control[{{y, 2., "cycle width"}, 0, 2, ImageSize -> Tiny}],
        spiralityC = Control[{{spirality, 0.}, 0, 1, ImageSize -> Tiny}],
        scaleC = Control[{{scale, 3.157, "sphere size"}, .00001, 15, ImageSize -> Medium}],
        settingsC = Control[{{settings, Take[allSettings, 4], "view"}, allSettings, ControlType -> TogglerBar}],
        opacityC = Control[{{opacity, .43}, 0, 1, Appearance -> "Labeled", ImageSize -> Small}],
        resetC = DynamicWrapper[
          Tooltip[Setter[Dynamic[reset], "reset"], "reset perspective", TooltipDelay -> .3],
          If[reset === "reset", (reset = False; vp = {1.3, -2.4, 2}; vv = {0, 0, 1}; {va, vc} = Automatic {1, 1})]]},
      
       With[{
         controls = Sequence[
           OpenerView[{"Style",
             Column[{
               Row[{
                 Column[{backgroundC, colorC}, Alignment -> Right],
                 Column[{effectC, thicknessC, opacityC}, Alignment -> Right]},
                Spacer[30]],
               Style[\[HorizontalLine], Lighter[LightGray]]}, Spacings -> 0]}],
           Row[{scaleC, Spacer[30], settingsC}],
           Row[{inflectionC, Spacer[30],
             Column[{angularityC, tensionC}],
             Column[{yC, spiralityC}], Spacer[30], resetC}]],
      
         storedVars = Sequence @@ ({{#, Automatic}, ControlType -> None} & /@ {vp, vv, va, vc}),
         dynamicView = Sequence[
           ViewPoint -> Dynamic[vp], ViewVertical -> Dynamic[vv],
           ViewAngle -> Dynamic[va], ViewCenter -> Dynamic[vc]]},
      
        (# /. switch[a_, b_] :> (*macro*)
             Unevaluated[Dynamic[If[MemberQ[settings, a], b, {}]]] &)@
         Manipulate[
          DynamicModule[{g, lines, riemannLines, invertedLines,
            \[Psi] = Round[Abs[FractionalPart[\[Zeta]]]*1., .25] /. {
               0. -> y, .5 -> y/2, .25 | .75 -> y/4}},
      
           g = ParametricPlot[
             1/(scale^1.4) (1. + spirality*(Log[\[Theta] + 1.] - 1.))*
              {\[Psi] Cos[\[Theta]] + x  Cos[64. \[Theta]] + (1 - effect*RandomReal[])*\[Zeta]* (Cos[512. \[Theta]] + Cos[64. \[Zeta] \[Theta]]),
               \[Psi] Sin[\[Theta]] - x Sin[64. \[Theta]] + (1 - effect*RandomReal[])* inflection*\[Zeta]* (Sin[512. \[Theta]] + Sin[64. \[Zeta] \[Theta]])}
      
             , {\[Theta], 0, 2 \[Pi]}, ImageSize -> {640, 480}, PerformanceGoal -> "Quality",
             Epilog -> {Gray, Thick, Circle[{0, 0}, 1]},
             PlotStyle -> Dynamic[{{color, Opacity[.43]}}], PlotRange -> Full,
             Background -> Dynamic[background], PlotPoints -> 270, Axes -> None];
      
           lines = Cases[g, Line[pts_] :> pts, Infinity];
           riemannLines = toRiemann /@ lines;
           invertedLines = invertRiemann /@ riemannLines;
           lines = Map[{##, 0} & @@@ # &, lines];
      
           ControlActive[g,
            Graphics3D[{
              switch[ "sphere", {Lighting -> "Neutral", Opacity[.1], Sphere[]}],
              switch["shuttle", shuttle], switch["cow", cow],
              Dynamic[color], Dynamic[Opacity[opacity]],
              Dynamic[AbsoluteThickness[thickness]],
              switch["original", Line /@ lines],
              switch["riemann", Line /@ riemannLines],
              switch["inverse", Line /@ invertedLines]}, Boxed -> False,
             dynamicView, Background -> Dynamic[background],
             ImageSize -> {640, 480}]]]
      
          , controls,
          storedVars,
          {{reset, "reset"}, ControlType -> None},
      
          Bookmarks -> bookmarks, Alignment -> Center]]]
      

This thing's almost as curly as my hair. Note that in Mathematica these aren't static renderings. They're regular Graphics3D panes that you can spin and move around every which way. But let's not forget why we're here:

  1. toRiemann = Compile[{{pts, _Real, 2}}, Module[{k},
        Map[(k = 2/(1 + #.#); {k #[[1]], k #[[2]], 1 - k}) &, pts]]];
    
    cow = {EdgeForm[None], Texture[Graphics[Disk[]]], Lighting -> "Neutral",
       Append[ExampleData[{"Geometry3D", "Cow"}, "GraphicsComplex"],
        VertexTextureCoordinates -> 1/500 ExampleData[{"Geometry3D", "Cow"}, "PolygonData"]]};
    
    riemannize[Graphics[g_, ___], options___] := Graphics3D[{
        (*sphere*){Lighting -> "Neutral", Opacity[.07], Sphere[{0, 0, 0}, 1]},
        (*cow*)Rotate[Rotate[Translate[cow, {0, 0, 1.2}], -Pi/2, {0, 0, 1}], -Pi/ 2, {1, 0, 0}],
        (*curves*)g /. (h : Line | Polygon)[pts_] :> {h[{##, 0} & @@@ pts], h[toRiemann[pts]]}},
       options, Boxed -> False];
    
    axiom = Polygon[2 {Cos[#], Sin[#]} & /@ (Pi/2 - 2 Pi Range[3]/3)];
    
    next[prev_] := prev /. Polygon[{p1_, p2_, p3_}] :> {
         Polygon[{p1, (p1 + p2)/2, (p1 + p3)/2}],
         Polygon[{p2, (p2 + p3)/2, (p1 + p2)/2}],
         Polygon[{p3, (p1 + p3)/2, (p2 + p3)/2}]};
    
    draw[n_] := Graphics[{Transparent, EdgeForm[Black], Nest[next, N@axiom, n]}];
    
    riemannize[draw[5], ViewPoint -> Top]
    

.


Look at the symmetry of this inversion:

  1. toRiemann = Compile[{{pts, _Real, 2}}, Module[{k},
        Map[(k = 2/(1 + #.#); {k #[[1]], k #[[2]], 1 - k}) &, pts]]];
    
    fromRiemann = Compile[{{pts, _Real, 2}}, Module[{k},
        Map[(-1/(#[[3]] - 1)) {#[[1]], #[[2]], 0} &, pts]]];
    
    cow = {EdgeForm[None], Texture[Graphics[Disk[]]], Lighting -> "Neutral",
       Append[ExampleData[{"Geometry3D", "Cow"}, "GraphicsComplex"],
        VertexTextureCoordinates -> 1/500 ExampleData[{"Geometry3D", "Cow"}, "PolygonData"]]};
    
    riemannTableau[Graphics[g_, ___], options___] := Module[{tmp},
       Graphics3D[{
         Translate[cow, {0, 0, 1.2}],
         {Opacity[.07], Sphere[{0, 0, 0}, 1]},
         g /. (h : Line | Polygon)[pts_] :> {
            EdgeForm[Opacity[.3]],
            (*original*)EdgeForm[Purple], Purple, h[{#1, #2, 0} & @@@ pts],
            (*riemann*)EdgeForm[Blue], Blue, h[tmp = toRiemann[pts]],
            (*riemann inverse*)EdgeForm[Red], Red, h[tmp = {#1, #2, -#3} & @@@ tmp],
            (*inverse*)h[fromRiemann[tmp]]}},
        options, Lighting -> "Neutral", Boxed -> False,
        Axes -> None, PlotRange -> All]];
    
    axiom = Polygon[{Cos[#], Sin[#]} & /@ (Pi/2 - 2 Pi Range[3]/3)];
    
    next[prev_] := prev /. Polygon[{p1_, p2_, p3_}] :> {
         Polygon[{p1, (p1 + p2)/2, (p1 + p3)/2}],
         Polygon[{p2, (p2 + p3)/2, (p1 + p2)/2}],
         Polygon[{p3, (p1 + p3)/2, (p2 + p3)/2}]};
    
    draw[n_] := Graphics[Nest[next, N@axiom, n]];
    
    riemannTableau[draw[5], ViewPoint -> {Top, Left}]
    

We have the original in the middle in purple, its Riemann mapping in blue, their inverses in red, and the cow in black and white. And except for the cow they all meet at the same three points. How gangster is that. (My friend's six-month old informs me that it's "substantially gangster" (paraphrasing)). But witness the scene of a zoom-out fractal:

  1. toRiemann = Compile[{{pts, _Real, 2}}, Module[{k},
        Map[(k = 2/(1 + #.#); {k #[[1]], k #[[2]], 1 - k}) &, pts]]];
    
    fromRiemann[pts_] :=
      Quiet@DeleteCases[(-1/(#3 - 1)) {#1, #2, 0} & @@@ pts,
        x_ /; MemberQ[x, Indeterminate]];
    
    cow = {EdgeForm[None], Texture[Graphics[Disk[]]],
       Append[ExampleData[{"Geometry3D", "Cow"}, "GraphicsComplex"],
        VertexTextureCoordinates -> 1/500 ExampleData[{"Geometry3D", "Cow"}, "PolygonData"]]};
    
    riemannTableau[Graphics[g_, ___], options___] := Module[{tmp},
       Graphics3D[{
             (*transparency on sphere causes weird graphical
              issue on my machine when n is high*)
         {Opacity[.07], Sphere[{0, 0, 0}, 1]},
         Translate[cow, {0, 0, 1.2}],
         g /. (h : Line | Polygon)[pts_] :> {
            (*original*)h[{#1, #2, 0} & @@@ pts],
            (*riemann*)h[tmp = toRiemann[pts]],
            (*riemann inverse*)h[tmp = {#1, #2, -#3} & @@@ tmp],
            (*inverted*)h[fromRiemann[tmp]]}},
        options, Lighting -> "Neutral", Boxed -> False,
        Axes -> None, PlotRange -> All]];
    
    next[prev_] := prev /. Polygon[{p1_, p2_, p3_}] :> {
         Polygon[{p1, (p1 + p2)/2, (p1 + p3)/2}],
         Polygon[{p2, (p2 + p3)/2, (p1 + p2)/2}],
         Polygon[{p3, (p1 + p3)/2, (p2 + p3)/2}]};
    
    {n, c} = {5, {0, -5/8; -1/4}};
    axiom = Polygon[c + {Cos[#], Sin[#]} & /@ (Pi/2 - 2 Pi Range[3]/3)];
    
    riemannTableau[
     Graphics[{EdgeForm[Black], Black, Nest[next, N@axiom, n]}],
     ViewVector -> {5 {-1, -1, 1}, {0, 0, 0}}]
    

At the chasm of infinity, our cow glances past its precipice, stares down its abyss. You know that machine in the Hitchhiker's Guide that explodes your mind or whatever by showing you how pathetically insignificant you are compared to the universe? Well this is like a Windows 3.1 version of that. Our poor cow friend's soul is being wrung on the very clothesline of endlessness itself. I think this is the first time I'm happy I'm not a cow.

...is what I would have said if this was any cow but this one.


    1. pat = Graphics[{Black, Disk[{0, 0}, 5], White,
          EdgeForm[{Black, Thickness[.03]}],
          Disk[{0, 0}, # + .07] & /@ Range[4, 1, -1],
          Black, Disk[{0, 0}, .15], Rectangle[{-4, 1.8}, {4, 2.1}],
          Rotate[Rectangle[{-4, 1.8}, {4, 2.1}], -Pi/4, {0, 0}],
          Rectangle[{-.2, -1.3}, {.2, -4}]}];
      
      (jhgn = {Lighting -> "Neutral", #1}) & @@
        SphericalPlot3D[u + v, {u, 0, Pi}, {v, 0, Pi},
         Mesh -> None, TextureCoordinateFunction -> ({#5, #4} &),
         PlotStyle -> Texture[pat]];
      
      
      xf1 = {
         {{0.0017206308062546146`, 0.0012959917814960697`, 0.0025851614902868744`},
          {0.0010674250900446086`, -0.0030803612062593683`, 0.0008337886519470046`},
          {0.0026876120080307113`, 0.0003937069230620502`, -0.0019861927280659755`}},
         {0.3257382788099915`, -0.03759999999999997`, 0.1862691804107692`}};
      
      xf2 = {
         {{0.0017206308062546146`, 0.0012959917814960697`, 0.0025851614902868744`},
          {-0.0010674250900446086`, 0.0030803612062593683`, -0.0008337886519470046`},
          {0.0026876120080307113`, 0.0003937069230620502`, -0.0019861927280659755`}},
         {0.3257382788099915`, 0.03759999999999997`, 0.1862691804107692`}};
      
      cow = {EdgeForm[None], Lighting -> "Neutral",(*Opacity[.999],*) Texture[Graphics[Disk[]]],            
         Append[ExampleData[{"Geometry3D", "Cow"}, "GraphicsComplex"],
          VertexTextureCoordinates -> 1/500 ExampleData[{"Geometry3D", "Cow"}, "PolygonData"]]};
      
      Graphics3D[{cow,
        GeometricTransformation[jhgn, xf1],
        GeometricTransformation[jhgn, xf2]},
       Boxed -> False]
      
    2. pat =(*ColorNegate@*)Graphics[{Black, Disk[{0, 0}, 5], White,
          EdgeForm[{Black, Thickness[.03]}],
          Disk[{0, 0}, # + .07] & /@ Range[4, 1, -1],
          Black, Disk[{0, 0}, .15], Rectangle[{-4, 1.8}, {4, 2.1}],
          Rotate[Rectangle[{-4, 1.8}, {4, 2.1}], -Pi/4, {0, 0}],
          Rectangle[{-.2, -1.3}, {.2, -4}]}];
      
      (jhgn = {Lighting -> "Neutral", #1}) & @@
        SphericalPlot3D[u + v, {u, 0, Pi}, {v, 0, Pi},
         Mesh -> None, TextureCoordinateFunction -> ({#5, #4} &),
         PlotStyle -> Texture[pat], PlotPoints -> 80(*0*)];
      
      sc = 1.5;
      xf1 = {
         sc {{0.010433915075096155`, 0.02050184708176941`, -0.0014526467022039392`},
           {0.014609900826184952`, -0.006252356712311273`, 0.01669614726190017`},
           {0.014456372783472383`, -0.008478490536514085`, -0.01582501134810082`}},
         {0.5447560973768777`, -0.5`, 0.5534681561478464`}};
      
      xf2 = {
         sc {{0.010433915075096155`, 0.02050184708176941`, -0.0014526467022039392`},
           {-0.014609900826184952`, 0.006252356712311273`, -0.01669614726190017`},
           {0.014456372783472383`, -0.008478490536514085`, -0.01582501134810082`}},
         {0.5447560973768777`, 0.5`, 0.5534681561478464`}};
      
      cow = {EdgeForm[None], Lighting -> "Neutral",(*Opacity[.999],*) Texture[ColorNegate@Graphics[Disk[]]],
         Append[ExampleData[{"Geometry3D", "Cow"}, "GraphicsComplex"],
          VertexTextureCoordinates -> 1/500 ExampleData[{"Geometry3D", "Cow"}, "PolygonData"]]};
      
      Graphics3D[{cow,
        GeometricTransformation[jhgn, xf1],
        GeometricTransformation[jhgn, xf2]},
       ViewPoint -> Right, Boxed -> False] // ColorNegate

This cow does not cower. Infinity cannot bully this bull, cannot bloviate this bovine. By all appearances this cow is wearing infinity on its mane. Its horns are probably made of $\small{\aleph _{\aleph_{\tiny{\ddots}}}}$ down 4 or 5 levels, an immutability surpassed only by that of the tusks of the Alephant. Our cow isn't staring into infinity. It's looking down at infinity, observing infinity with detached understanding. If our cow were not so enlightened, and also had the facial muscles, it might betray the subtlest of smiles at infinity's infinity face, for infinity's turbid fractal whirlpools and vast lethargic swamps are but swathes of data like any other to this cow.

Long ago, having mastered the magesterial tetrafecta of science, mathematics, spirituality, and politics, our cow stepped hoof outside Farmer Joe's farm and set out on an adventure of like, just so much awesome. One of its side gigs these days is being the final observer of our domain, preventing our section of the Great Algorithm from backtracking by stellating through the cosmos our most entwined entwinements. I think this is the first time I'm jealous of a cow.

In any case, as you can see the Riemann sphere is pretty useless. But while we're on the subject of 3D let's see how our various approaches do here. Chaos game:

  1. draw[vertices_, numPoints_, options___] :=
      Graphics3D[{Lighter[Green],
        Sphere[FoldList[(#1 + #2)/2 &, First[N@vertices],
          RandomChoice[N@vertices, numPoints]], .001]},
       Lighting -> {{"Point", LightYellow, Scaled[{1, 1, 1}], 5}},
       options, Boxed -> False];
    
    vertices = PolyhedronData[{"Pyramid", 3}, "VertexCoordinates"];
    draw[vertices, 100000,
     ViewPoint -> {0, 0, Infinity},
     ViewVertical -> {1, 0, 0}]
    
  2. draw1[vertices_, numPoints_, options___] :=
      Graphics3D[{Lighter[Green], EdgeForm[None],
        Cuboid[#, # + .01] & /@ FoldList[(#1 + #2)/2 &, First[vertices],
          RandomChoice[N@vertices, numPoints]]},
       Lighting -> {{"Point", LightYellow, Scaled[{1, 1, 1}], 5}},
       options, Boxed -> False];
    
    draw[vertices_, numPoints_, options___] :=
      Graphics3D[{Lighter[Green],
        Sphere[FoldList[(#1 + #2)/2 &, First[N@vertices],
          RandomChoice[N@vertices, numPoints]], .001]},
       Lighting -> {{"Point", LightYellow, Scaled[{1, 1, 1}], 5}},
       options, Boxed -> False];
    
    vertices = PolyhedronData[{"Pyramid", 3}, "VertexCoordinates"];
    
    (*1*)
    Defer[AbsoluteOptions][draw1[vertices, 20000, ImageSize -> Medium]]
    
    (*2*)
    draw[vertices, 2000000,
       (* ViewPoint, ViewVertical from (*1*) *)
       Method -> {"ShrinkWrap" -> True},
       ImageSize -> 2 1280] // Rasterize // ImageResize[#, Scaled[1/4]] &
    

This is using little spheres as the points. You could use pyramids or anything else instead. Even go back to nature and use actual points. It's a bit tricky to get decent images since the chaos game doesn't place points in a regular arrangement, so you need a large number of points. Each of these images uses 2 million spheres and takes about 10 minutes to render on my little laptop.

This top view shows one of the symmetries that appear in the 3D triangle. This side view shows another. And a top view of the 4-corner pyramid. These symmetries are interesting because they appear absolutely no different than 2D renditions (for example). At first this seems mysterious, since the symmetries appear from every which angle. But the reason it happens is because our distance function works on each coordinate independently:

$$p_{i+1}=\frac 1 2 (p_i+v)$$

This formula is applied to the $x$ $y$ $z$ coordinates separately. So we could chop off any one of the coordinates from all the points in the 3D Sierpinski triangle to get a regular 2D Sierpinski triangle. And more to the point, what we're really doing is geometric: finding the point halfway between two given points. This $\small{\frac 1 2 (p_i+v)}$ formula is just a particular algebraic statement of that.

In other words, the geometry of the algorithm doesn't care about our coordinate system, so we're going to get the projected equivalent of a 2D rendering from any angle we pick, not just from the $x$ $y$ $z$ cross sections of our computation (certainly there is a linear algebra term for this).

To make clear what we're talking about, this is the chaos game on a prism, and the same thing from the same viewpoint, except with the 3D projection effect removed. As you can see, the 'hidden dimension' has no offect on what is seen. Au contraire messieur. If it did have an effect, that would be interesting.

Something I noticed though is that while we can remove a coordinate, we can't add a coordinate, in the sense that, for example, there's no way to combine independent $x$ $y$ streams to create a Sierpinski triangle. For our 2D Sierpinski triangle, there's something to the fact that a single point is specified by two coordinates instead of just one.

I think there may be an interesting statistical or information-theoretic interpretation to this. I'm not really familiar with either of these subjects though. Geometric approach:

  1. draw[shapeName_, n_, options___] := Module[{shape, next},
       shape = PolyhedronData[shapeName, "Faces"];
    
       (*scale by 1/2 toward each vertex,in turn*)
       next[prev_] := Scale[prev, 1/2, #] & /@ shape[[1]];
    
       Graphics3D[{EdgeForm[Opacity[.15]], Nest[next, N@shape, n]},
        options, Lighting -> "Neutral", Boxed -> False]];
    
    Grid[Table[
      draw[{"Pyramid", k}, n, ViewPoint -> {0, 0, Infinity}],
      {k, 3, 5}, {n, 0, 3}]]
    

Behold the Lemon Lime Fortress. Throw in a few salt blocks, pour some Corona at the top, join the party at the base. To make our lives one notch easier, our code takes advantage of Mathematica's built-in transformation infrastructure, in this case the symbol Scale. It also pulls the geometry of things from our good friend Mr. PolyhedronData. The nice thing about having such a general setup is that we can readily apply this geometric fractalization on arbitrary shapes:

  1. draw[shapeName_, n_] := Module[{shape, next},
       shape = PolyhedronData[shapeName, "Faces"];
       next[prev_] := Scale[prev, 1/2, #] & /@ shape[[1]];
    
       Graphics3D[Nest[next, N@shape, n],
        Method -> {"ShrinkWrap" -> True},
        Lighting -> "Neutral", Boxed -> False]];
    
    shapes = {"TruncatedIcosahedron", "TriakisIcosahedron", "TetrakisHexahedron",
       "SmallStellatedDodecahedron", "ElongatedPentagonalCupola", "Icosahedron",
       "ElongatedSquareDipyramid", "DuerersSolid"};
    
    Grid[Partition[#, 2]] &[
     Table[Tooltip[Panel[#], shape] &@
       Row[Table[draw[shape, n], {n, 0, 1}], Spacer[30]],
      {shape, shapes}]]
    

Don't ask me what the hell that last shape is. I figure it just managed to stow away into PolyhedronData somehow, like the semiconscious pre-sentient kernel of a future Skynet. The faces of these shapes show very clearly that we get 2D slices for free, like in these perspectives from below (we aren't cheating here). The edges by themselves make pretty patterns:

  1. draw[shapeName_, n_, options___] := Module[{shape, next, axiom},
       shape = PolyhedronData[shapeName, "Faces"];
       next[prev_] := Scale[prev, 1/2, #] & /@ shape[[1]];
    
       axiom = {shape, If[showLittleBalls,
          {FaceForm[{Opacity[.85], White}],
           Glow[Green], Sphere[{0, 0, 0}, .09]}]};
    
       Graphics3D[{Transparent, EdgeForm[{Opacity[opacity], color}],
         Nest[next, N@axiom, n]}, options, Lighting -> "Neutral",
        Boxed -> False]];
    
    shapes = {"TruncatedIcosahedron", "TriakisIcosahedron", "TetrakisHexahedron",
       "SmallStellatedDodecahedron", "ElongatedPentagonalCupola", "Icosahedron",
       "ElongatedSquareDipyramid", "DuerersSolid"};
    
    {color, opacity, showLittleBalls} = {Black, .6, False};
    Grid[Partition[#, 2]] &@
     Table[Tooltip[#, shape] &@
       Row[Table[draw[shape, n, ViewPoint -> Top], {n, 0, 1}], Spacer[30]],
      {shape, shapes}]
    

To make sure that after all this scrolling we're still on the same web page, this is our chaos game algorithm:

1 start at any point. call it p
2 pick a vertex at random
3 find the point halfway between p and that vertex
4 call that point p and draw it
5 goto 2

The only difference between 2D and 3D versions of this algorithm is having 3 coordinates instead of 2. Just as in 2D, we can alter step 3 in various ways. The simplest is to move not halfway towards the chosen vertex, but .25 or .7 of the way, etc:

  1. draw[vertices_, df_, numPoints_, options___] := Graphics3D[{
        Opacity[.5], PointSize[0],
        Point[FoldList[df, First[N@vertices],
          RandomChoice[N@vertices, numPoints]]]},
       (*Method->{"ShrinkWrap"->True},*)
       options, PlotRange -> Automatic, Boxed -> False];
    
    functions = Function[r, r (#1 + #2) &] /@ {1, .96, .6, .5, .2};
    
    Grid[Join[
      {TraditionalForm[Trace[#[a, b]][[2]]] & /@ functions},
      ParallelTable[
       draw[PolyhedronData[{"Pyramid", v}, "VertexCoordinates"],
        df, 50000, ViewPoint -> {Front, Top}],
       {v, 3, 5}, {df, functions}]]]
    

Those odd random walks are because the 4- and 5-pyramids have Mean[vertices] != {0, 0, 0}. One thing I noticed is that random walks resemble the outlines of continents. How curious. I wonder if it boils down to the self-similarity of the Brownian motion of water molecules, or something of the like. I.e. the idea that if our continents were surrounded by materials which did not move Brownianly, our coastlines would have different kinds of shapes. Remember that we can get creative with our distance function:

  1. draw[vertices_, df_, numPoints_, options___] :=
      Graphics3D[{PointSize[0], Opacity[.3],
        Point[FoldList[df, RandomReal[{0, .0001}, 3],
          RandomChoice[N@vertices, numPoints]]]},
       (*Method->{"ShrinkWrap"->True},*)
       options, Boxed -> False];
    
    rotate = RotationTransform;
    functions = {
       (#1 + #2)/RandomChoice[Prime[Range[3]]] &,
       (#1 + #2)/RandomChoice[Prime[Range[3]]!] &,
       (#1 + #2)/RandomChoice[Prime[Range[10]]] &,
       #1 + .5 rotate[10. Degree, {#1, #2}][#2 - #1] &};
    
    Grid[Join[
      {TraditionalForm[Trace[#[a, b]][[2]]] & /@ functions},
      ParallelTable[
       draw[PolyhedronData[{"Pyramid", v}, "VertexCoordinates"], df, 5000],
       {v, 3, 5}, {df, functions}]]]
    

Keep in mind that in Mathematica these are all interactive 3D panes. Since I associate these kinds of sparse fractal distributions with the distribution of matter through the scales of the cosmos, flying through these point structures engages my Björkian semi-spiritual naturalistic side. :) Your mileage may vary. Our old logarithmic distance function can be applied in 3D as well. For two points $a$ and $b$, with $d$ the Euclidean distance and $w$ a specific number between 0 and 1 (though not necessarily), the distance function is:


$$(a+b) \log (d(a,b)+w)$$
  1. game = Compile[{{vertices, _Real, 2}, {w, _Real}, {numpoints, _Integer}},
       Module[{diff},
        FoldList[(diff = #2 - #1;
           Clip[(#1 + #2) Log[Sqrt[diff.diff] + w]]) &,
         {0, 0, 0}, RandomChoice[vertices, numpoints]]]];
    
    draw[vertices_, w_, numPoints_, options___] :=
      Graphics3D[{PointSize[0], Opacity[7 .05],
        Point[game[vertices, w, numPoints]]},
       options, Boxed -> False];
    
    Needs["PolyhedronOperations`"];
    vertices = Stellate[PolyhedronData[{"Pyramid", 5}, "Faces"]][[1]];
    
    draw[vertices, .2, 600000, PlotRange -> All(*,
     Method->{"ShrinkWrap"->True}*)(*,ViewPoint->{Infinity,0,0}*)]
    

These pictures differ by $w$ factor, viewpoint, or the set of vertices on which the game is being played. For most of these I'm using the vertices of regular polyhedra from PolyhedronData. Note that the vertices of the game are not necessarily in proportion to the figure itself.

At this point I should remention that all of the code snippets on this page are self-contained. If you have Mathematica you can copy-paste this and start producing these figures, which, I should also remention, are interactive 3D models. I'm a big fan of black ink on white paper, and these are like being able to change the perspective of a pure ink painting in real time. Teknikara no jutsu.

  1. game = Compile[{{vertices, _Real, 2}, {w, _Real}, {numpoints, _Integer}},
       Module[{diff},
        FoldList[(diff = #2 - #1;
           Clip[(#1 + #2) Log[Sqrt[diff.diff] + w]]) &,
         {0, 0, 0}, RandomChoice[vertices, numpoints]]]];
    
    draw[vertices_, w_, numPoints_, options___] :=
      Graphics3D[{PointSize[0], Opacity[7 .05],
        Point[game[vertices, w, numPoints]]},
       options, Boxed -> False];
    
    proc[img1_, cf_: ColorData[1], mode_: None, blur_: 8] :=
      Module[{img, components, rank, largest, colored},
       img = RemoveAlphaChannel[ColorNegate@ColorConvert[img1, "Grayscale"]];
       components = MorphologicalComponents[img];
    
       Module[{measurements, sorted},
        measurements = ComponentMeasurements[components, "Count"];
        sorted = First /@ Reverse@SortBy[measurements, Last];
        rank[label_] := (rank[label] = Position[sorted, label][[1, 1]])];
    
       colored = Colorize[components,
         ColorFunction -> (cf[rank[#]] &), ColorFunctionScaling -> False];
    
       If[mode == "Angelic",
        colored = ImageMultiply[img, colored]];
    
       ColorNegate[ImageMultiply[ColorNegate[img],
           Blur[#, blur] &@ColorNegate[colored]]] // ImageAdjust];
    
    Needs["PolyhedronOperations`"];
    vertices = OpenTruncate[PolyhedronData[{"Pyramid", 3}, "Faces"]][[1]];
    g = draw[vertices, .5, 600000, Method -> {"ShrinkWrap" -> True}];
    
    proc[g, # /. Join[
        Thread[Range[4] -> {Red, Green, Green, Green}],
        {_ -> Lighter[Green]}] &]
    

Some of thes are like alien Rorschach tests. Like what do you see in this one? I see a mosquito that can suck the lifeblood out of your soul. This one, however, is definitely from an as-yet unreleased Matrix film. And we also have the Minotaur's armor and his shield of Cancer. I'd recognize my buddy's armor in even the most obtuse alien Rorschachs. See also a stereographic projection of one of the rooms of Asterion's maze and an aspect, which needs no explanation.

The originals are 3D but this coloring is a 2D image process. It highlights components of the image based on their sizes. So if your image has 3 large blobs with dozens of tiny blobs all around, you can use, for example, # /. {1 -> Red, 2 -> Green, 3 -> Yellow, _ -> Pink} & to color the big blobs specific colors and all other blobs pink. Though in most of these images I only use one or two colors.

  1. game = Compile[{{vertices, _Real, 2}, {numPoints, _Integer}, {wowzerz, _Real}},
       Module[{diff},
        FoldList[(diff = #2 - #1;
           Clip[(#1 + #2) Log[Sqrt[diff.diff] + wowzerz]]) &,
         {0., 0., 0.}, RandomChoice[vertices, numPoints]]]];
    
    {numFrames, imageSize, numPoints} = {
        {5(*sec*)15(*fps*), 1/2 {640, 480}, 1/3 600000},
        {20(*sec*)15(*fps*), {640, 480}, 600000}}[[1]];
    
    Needs["PolyhedronOperations`"];
    vertices = Stellate[PolyhedronData[{"Pyramid", 5}, "Faces"]][[1]];
    
    frame = Function[w,
       Graphics3D[{Opacity[.1], PointSize[0], Point[game[vertices, numPoints, w]]},
        ImageSize -> imageSize, ViewVertical -> {0, 0, 1}, Boxed -> False, SphericalRegion -> True, PlotRange -> 1,
        ViewVector -> {RotationTransform[2 Pi w, {0, 0, 1}][{1, 0, (w - .25) Pi/2}], {0, 0, 0}}]];
    
    SetDirectory["c:/users/zrp/desktop/frames"];
    
    range = Range[0, 3/4, 3/4/(numFrames - 1)];
    file[w_] := ToString[N@w] <> ".png";
    
    ParallelDo[
      If[! FileExistsQ[file[w]],
       Export[file[w], frame[w]]],
      {w, range}];
    
    Export["mov.avi",
     ColorNegate /@ ImageAdjust /@ Import /@ file /@ range]
    
    Beep[];
    Button[open, SystemOpen["mov.avi"],
     Enabled -> FileExistsQ["mov.avi"]]
    
    1. MovieMaker[frameF_, range : {start_Integer, stop_Integer}, rest___] :=
        MovieMaker[frameF, {start, stop, stop - start}, rest];
      
      (*arithmetic for eg doubling movie length is easier by 'intervals' than by 'frame count'*)
      MovieMaker[frameF_, range : {start_, stop_, numIntervals_}, rest___] :=
        MovieMaker[frameF, List[Range[#1, #2, (#2 - #1)/#3(*(#3-1)*)] & @@ range], rest];
      
      MovieMaker::expqq = "Export is complaining about something. " <>
         "Most likely you're feeding it items with different image sizes.";
      
      MovieMaker::usage =
        "NOTE: copies of this notebook are automatically stored along
        with the generated files. To prevent this, set AutoArchive -> False.
        
        MovieMaker[frameFunction, rangeSpec, options___]
        
        rangeSpec:
        {start, stop, number of intervals}:  {0, 1, 5(*sec*)15(*fps*)}
        {start, stop} integer range:  {1, 20}
        {explicit list}:  {AstronomicalData[\"Earth\",\"OrbitPath\"][[1]]}
        
        The Label option determines the folder name under which the animation
        is created. For example, if changing a variable X makes a different
        animation, then place that variable in the Label spec so that when you
        change that variable, the animation will be generated in a different folder.
        
        Likewise, the first element of the Process spec determines the folder
        and uniqueness of the process function. Processes work in subfolders of
        the main project folder, meaning you can experiment with multiple processes
        in a single project.
        
        MovieMaker[
         {ToLowerCase[#], ToUpperCase[#]} &, {CharacterRange[\" \", \"~\"]},
         Serialization -> Hash, Label -> \"UpperLower\", FileTypes -> {\".mx\", \".png\", \".gif\"},
         Process -> {\"times\", ImageMultiply @@ Map[Rasterize[#, ImageSize -> 400 {1, 1}] &, #] &},
         MovieOptions -> {\"DisplayDurations\" -> 1}, MapFunction -> Map]
        
        Serialization is for converting values to valid file names.
        MapFunction is for when you don't want to use parallelization.
        Directory setting specifies the specific project folder, overriding Label.";
      
      Options[MovieMaker] = {
         Label -> Automatic, Process -> {None, None}, MapFunction -> ParallelMap, AutoArchive -> True,
         FileTypes -> {".png", ".png", ".avi"}, MakeMovie -> True, MovieOptions -> {}, Directory -> Automatic,
         Ordering -> (BlockRandom[RandomSample[#]] &), Serialization -> Composition[List, Chop, N]};
      
    2. (* After I wrote this program, a more powerful approach occurred to me. We could have a
      macro that would be used something like this: *)
      FileBackedProcess[Function[val,
         a = S[1][Rasterize@dirp[val]];
         b = S[2][Rasterize@derp[val]];
         S[3][ImageMultiply[a, b]]]];
      (* where the S[i_][body_] are the momoization points into the file system. If the S finds
      the file corresponding to the [i][body], then the file is imported. Otherwise it executes
      the body and saves the file. The point would be to make the file aspect as
      easy as annotating things with S[i] *)
      
      MovieMaker[frameF_, List[valueList_List], OptionsPattern[]] := Module[{
         tooltip, mainLabel, processLabel, processF, mapF, frameExt, processedExt, movieExt, dir,
         framesDir, processedDir, movieFile, fileMap, numFrames, alive = True, folder0exists,
         foldersExistL, folder1exists, folder2exists, progress1, progress2, movieDone, makeFrames,
         processFrames, makeMovie, serialization, archive, makeMovieA, preview, printPreview, printFileMap},
      
        tooltip[expr_] := Tooltip[#, expr, TooltipDelay -> .25] &;
        {mainLabel, mapF, makeMovieA, serialization} =
         OptionValue[{Label, MapFunction, MakeMovie, Serialization}];
      
        {processLabel, processF} =
         Replace[OptionValue[Process], {
           {pf_} :> {ToString[pf], pf},
           pf : Except[_List] :> {ToString[pf], pf}}];
      
        {frameExt, processedExt, movieExt} = PadRight[
          Flatten[List[OptionValue[FileTypes]]], 3,
          FileTypes /. Options[MovieMaker]];
      
        mainLabel = Replace[mainLabel,
          Automatic -> IntegerString[Hash[{frameF, valueList}, "CRC32"], 36]];
        dir = Replace[OptionValue[Directory], Automatic ->
           FileNameJoin[{NotebookDirectory[], "vids", ToString[mainLabel]}]];
      
        framesDir = FileNameJoin[{dir, "frames"}];
        processedDir = FileNameJoin[{dir, "processed", ToString[processLabel]}];
        movieFile = FileNameJoin[{dir, ToString[{processLabel, mainLabel}] <> movieExt}];
      
        (* main iteration construct *)
        fileMap[f_, vals_: valueList, map_: mapF] := map[Function[val,
           f[
            FileNameJoin[{framesDir,
              ToString[serialization[val]] <> frameExt}],
            FileNameJoin[{processedDir,
              ToString[serialization[val]] <> processedExt}],
            val]],
          vals];
      
        numFrames = Length[valueList];
        progress1 = Total@Boole[fileMap[FileExistsQ[#1] &]];
        progress2 = Total@Boole[fileMap[FileExistsQ[#2] &]];
        foldersExistL = FileExistsQ /@ {dir, framesDir, processedDir};
        movieDone = FileExistsQ[movieFile];
        SetSharedVariable[progress1, progress2];
      
        If[OptionValue[AutoArchive] && FileExistsQ[dir] &&
          ! FileExistsQ[FileNameJoin[{dir, ToString[mainLabel] <> ".nb"}]],
         Export[FileNameJoin[{dir, ToString[mainLabel] <> ".nb"}],
          NotebookGet[EvaluationNotebook[]]]];
      
    3. (**)
      makeFrames[] := (
         Quiet@CreateDirectory[framesDir];
         foldersExistL[[1 ;; 2]] = {True, True};
         If[OptionValue[AutoArchive],
          Export[FileNameJoin[{dir, ToString[mainLabel] <> ".nb"}],
           NotebookGet[EvaluationNotebook[]]]];
      
         fileMap[If[! FileExistsQ[#1],
            Export[#1, frameF[#3]];
            progress1++] &,
          OptionValue[Ordering][valueList]]);
      
      (**)
      processFrames[] := If[
         processF =!= None,
         Quiet@CreateDirectory[processedDir];
         foldersExistL[[3]] = True;
         If[OptionValue[AutoArchive],
          Export[FileNameJoin[{processedDir, ToString[{mainLabel, processLabel}] <> ".nb"}],
           NotebookGet[EvaluationNotebook[]]]];
      
         fileMap[If[! FileExistsQ[#2] && FileExistsQ[#1],
            Export[#2, processF[Import[#1]]];
            progress2++] &,
          OptionValue[Ordering][valueList]]];
      
      (**)
      makeMovie[] := If[makeMovieA,
         If[FileExistsQ[movieFile],
          Print["movie file already exists"],
          With[{ab = If[processF === None, #1, #2]},
           If[And @@ fileMap[FileExistsQ[ab] &],
            Check[
              Export[movieFile, fileMap[Import[ab] &],
               Sequence @@ OptionValue[MovieOptions]];
              movieDone = True, Message[MovieMaker::expqq];
              movieDone = False, {Export::errelem}]]]]];
      
      (**)
      preview[] := preview[RandomChoice[valueList]];
      preview[val_] := Module[{frame, fileName, tempFile},
         tempFile = FileNameJoin[{$TemporaryDirectory, ToString[Hash[val]] <> frameExt}];
         fileName = First@fileMap[#1 &, {val}];
      
         If[FileExistsQ[fileName],
          (**)frame = Import[fileName],
          (**)frame = Import[Export[tempFile, frameF[val]]];
          Print[Labeled[frame, N@val, Right]]; Beep[]];
      
         If[processF =!= None,
          Print[Labeled[processF[frame], N@val, Right]]; Beep[]]];
      
      (**)
      printPreview[] := CellPrint[ExpressionCell[Defer[
           preview[Placeholder["val"]]], "Input"]];
      
      (**)
      printFileMap[] := CellPrint[ExpressionCell[Defer[
           frames2 = fileMap[If[FileExistsQ[#2], Import[#2], Sequence @@ {}] &];],
          "Input"]];
      
      (**)
      archive[] := Module[{fileName},
         fileName = ToString[mainLabel] <> " " <>
           DateString[{"DateShort",
             " (", "Hour12", " ", "Minute", " ", "Second", " ", "AMPM", ")"}];
      
         Quiet@CreateDirectory[dir];
         foldersExistL[[1]] = True;
         Export[FileNameJoin[{dir, fileName <> ".nb"}],
          NotebookGet[EvaluationNotebook[]]];
      
         Beep[]];
      
    4. (*controls*)
        With[{
          btnMakeFrames = Button["frames + process + movie",
            makeFrames[]; Beep[]; processFrames[]; Beep[]; makeMovie[]; Beep[],
            Method -> "Queued", Enabled -> Dynamic[progress1 =!= numFrames]],
          btnProcessFrames = Button["process + movie",
            processFrames[]; Beep[]; makeMovie[]; Beep[],
            Method -> "Queued", Enabled -> Dynamic[
              progress2 =!= numFrames && progress1 =!= 0 && processF =!= None]],
          btnMakeMovie = Button["movie",
            makeMovie[]; Beep[],
            Method -> "Queued", Enabled -> Dynamic[
              (progress2 === numFrames ||
                 (processF === None && progress1 === numFrames)) &&
               ! movieDone && makeMovieA]],
          btnMainFolder = tooltip["open folder"]@
            Button[{mainLabel, processLabel}, SystemOpen[dir],
             Enabled -> Dynamic[foldersExistL[[1]]]],
          btnFramesFolder = tooltip["open folder"]@
            Button[{Dynamic[progress1]/ToString[numFrames],
              ProgressIndicator[Dynamic[progress1/numFrames]]},
             SystemOpen[framesDir],
             Enabled -> Dynamic[foldersExistL[[2]]]],
          btnProcessFolder = tooltip["open folder"]@
            Button[{Dynamic[progress2]/ToString[numFrames],
              ProgressIndicator[Dynamic[progress2/numFrames]]},
             SystemOpen[processedDir],
             Enabled -> Dynamic[processF =!= None && foldersExistL[[3]]]],
          btnMovieFile = tooltip["open movie"]@
            Button[{Dynamic[Boole[movieDone]]/"1",
              ProgressIndicator[Dynamic[Boole[movieDone]/1]]},
             SystemOpen[movieFile], Enabled -> Dynamic[movieDone]]},
      
         (*without going the extra mile, better to have no persistence*)
         Dynamic[If[alive === True,
             Panel[#, FrameMargins -> {{Automatic, Automatic}, {Automatic, 0}}],
             Panel[Tooltip[Overlay[{
                 Style["VWXYZ", Lighter[LightGray, 2/3], FontFamily -> "Wingdings"],
                 Style["dead", Darker[Red, 1/6]]}, All, 2, Alignment -> {Center, Center}],
               "R.I.P. this MovieMaker module"],
              FrameMargins -> 0]]] &@
      
          Manipulate[
           Grid[{
             {btnMainFolder, SpanFromLeft},
             {btnMakeFrames, btnFramesFolder},
             {btnProcessFrames, btnProcessFolder},
             {btnMakeMovie, btnMovieFile}}],
      
           Bookmarks :> {
             "preview" :> AbortProtect[preview[]],
             Overscript[Row[{"print ", Style["preview", Bold], " function"}], ""] :> printPreview[],
             Row[{"print ", Style["fileMap", Bold], " function"}] :> printFileMap[],
             Overscript["write archive", ""] :> archive[],
             "shoot" :> (alive = False)},
           Paneled -> False, FrameMargins -> False]]];
      
  2. game = Compile[{{vertices, _Real, 2}, {numPoints, _Integer}, {wowzerz, _Real}},
       Module[{diff, b},
        (*NestList for less memory usage. i didn't actually verify this*)
        NestList[(
           b = RandomChoice[vertices];
           diff = b - #1;
           Clip[(#1 + b) Log[Sqrt[diff.diff] + wowzerz]]) &,
         {0, 0, 0}, numPoints]]];
    
    proc[img1_, cf_: ColorData[1], mode_: None, blur_: 8] :=
      Module[{img, components, rank, largest, colored},
       img = RemoveAlphaChannel[ColorNegate@ColorConvert[img1, "Grayscale"]];
       components = MorphologicalComponents[img];
    
       Module[{measurements, sorted},
        measurements = ComponentMeasurements[components, "Count"];
        sorted = First /@ Reverse@SortBy[measurements, Last];
        rank[label_] := (rank[label] = Position[sorted, label][[1, 1]])];
    
       colored = Colorize[components,
         ColorFunction -> (cf[rank[#]] &), ColorFunctionScaling -> False];
       If[mode == "Angelic",
        colored = ImageMultiply[img, colored]];
       ColorNegate[ImageMultiply[ColorNegate[img],
           Blur[#, blur] &@ColorNegate[colored]]] // ImageAdjust];
    
    Needs["PolyhedronOperations`"];
    vertices = OpenTruncate[PolyhedronData["Icosahedron", "Faces"]][[1]];
    vertices = Rescale[vertices] - 1/2; (*rescale to 1/2 {-1, 1} range*)
    
    {numFrames, imageSize, numPoints} = {
       {5(*sec*)15(*fps*), {16, 9} (360/9), 600000},
       {5(*sec*)15(*fps*), {16, 9} (1080/9), 10000000}}[[2]];
    
    label = {"NUCLEAR1080P", numPoints, IntegerString[Hash[vertices, "CRC32"], 36]};
    
    process = {
        {"[COLORDATA3]", Composition[
          proc[#, If[# == 1, Blue, ColorData[3][#]] &, "Angelic", 1] &,
          ImageResize[#, Scaled[1/2]] &, Blur[#, 1] &, ImageAdjust]},
        {"[HIGHBLUR]", Composition[
          proc[#, If[# == 1, Blue, ColorData[3][#]] &, "Angelic", 40] &,
          ImageResize[#, Scaled[1/2]] &, ImageAdjust]}}[[2]];
    
    frame[w_] :=
      Graphics3D[{Opacity[.1], PointSize[0],
        Point[game[vertices, numPoints, w]]},
       ImageSize ->(**)2(**)imageSize, ViewVertical -> {0, 0, 1}, Boxed -> False,
       SphericalRegion -> True, Method -> {"ShrinkWrap" -> True},
       ViewVector -> {RotationTransform[2 Pi w, {0, 0, 1}][{1, 0, (w - .25) Pi/2}], {0, 0, 0}}];
    
    MovieMaker[frame, {.4, .75, 4 numFrames},
     Label -> label, Process -> process]
    

This just need some James Horner music. And science majors, witness a dangerous nuclear science experiment gone horribly awesome. These are animations on the $w$ factor. For the source, you can just use the basic code. But if you intend to do more general experimentation, then something like my little MovieMaker utility will be useful. It's a quite general utility. Each of these movies took something on the order of 20 hours for my computer to make. That's why having a minimal-fuss setup is convenient.

As for the renderings and animations themselves, they're basically me chewing a few times on one of the leaves of one of the branches of a tree I happened to run up the side of like a monkey. There's a lot of trees in this jungle to monoperambulate.

What's great about these structures is that they are still fractals. They may look spazzy and some of them may remind you of Vash the Stampede's plant mode arm, but they possess self-similarity throughout. For example, why do the arms of the nest look like that?

It's because the nest as a whole looks like that. And notice that as the big bird flies in from below to explode into the nest, the little birds all around the nest follow along (because adults know best) and explode into their own little nests, and so on, producing the distinctive infinitary echela of simultaneously exploding dinosaur progeny. And notice that the big bird itself is a version of the entire figure. Now, as for the hat, who knows.

The chaos game is an algorithm that we use for the sake of computational convenience. The "real" algorithm doesn't randomly pick among the vertices, it takes every point toward every observer at each step. And it's actually easy to see how the self-similarity of the algorithm comes about. Look here at a house and an observer:

If we run one step of the "real" algorithm, we get this. Something interesting here is that there is no difference between what the observer sees in either case. The little house is exactly blocking his or her view of the bigger house, like an inescapable mathematical version of a really tall person sitting in front of you at a theatre (formally we would say the houses are cosyzygous). If we start with two observers, we get this then this.

So it's clear that the scaled-skewed self-similarity is inevitable. What I find interesting is that for a given set of vertices and distance function, the resulting figure as a whole is also inevitable. You can start the chaos game at any point (or points, because e.g. the $\small{\frac 1 2}$ factor effectively shrinks your whole house into a point) and you will end up with the same figure, just a different approximation of it.

Another way of thinking about it is that the resulting figure is precisely the figure that all observers "agree on":



Because running the full algorithm on the entire figure does nothing. I.e. the figure is the fixed point of the algorithm. This automagic consensusing bonks my head and seems to me to carry a particular philosophical undertone... over which I shalln't digress.

Mathematically, it appears our chaos game shennaneganery as a whole falls under the contraction mapping principle. Tersely complicated explanations of inconfusably simple things not withstanding, I know me some topology but not enough to understand the bigger picture of what's going on.

On the subject of hats, when going from 2D L-systems to 3D L-systems I had to put a hat on the turtle and also give it the ability to do backflips and taco rolls:

Even wearing Mugen's shoes. Wow. Unfortunately, as epic as this is, with our current technology we're limited to e.g. representing the turtle's hat with an abstraction called a "vector", which certainly doesn't connote the same social status or sophistication. Still it's enough for some 3D L-systems, such as this version of the arrowhead construction:

    1. Module[{options = {
           Axiom -> None, Rules -> {}, Iterations -> 1, Definitions -> {},
           DrawStyle -> {}, HatStyle -> {}, Primitive -> Tube, TraceHat -> False,
           HatWorldplaneStyle -> Directive[EdgeForm[None], Opacity[.2]],
           HatPrimitive -> Composition[Arrow, Tube], Angle -> 2. Pi/6,
           RandomStuff -> Sphere[{0, 0, 0}, .05]}},
      
        SetAttributes[Draw, Orderless];
      
        Draw[commands : {Except[_Rule | _RuleDelayed] ..},
          rules : {(_Rule | _RuleDelayed) ..}, rest___] := Draw[Axiom -> commands, Rules -> rules, rest];                
        Draw[rules : {(_Rule | _RuleDelayed) ..}, rest___] := Draw[Rules -> rules, rest];
        Draw[commands : {Except[_Rule | _RuleDelayed] ..}, rest___] := Draw[Axiom -> commands, rest];
        Draw[opts : OptionsPattern[Join[Options[Graphics3D], options]]] :=
         Module[{commands, reshape, states, points, hatTrace, hatWorldplane,
           forwardP, leftP, frontflipP, tacoleftP, flipoutP, pushI, popI, definitionsI},
      
          (*basic parameterized state transfer functions*)
          forwardP[p_][{z_, face_, hat_}] := {z + p face, face, hat};
          leftP[p_][{z_, face_, hat_}] := {z, RotationTransform[p, hat][face], hat};
          tacoleftP[p_][{z_, face_, hat_}] := {z, face, RotationTransform[p, face][hat]};
          frontflipP[p_][{z_, face_, hat_}] := Module[{rot},
            rot = RotationTransform[p, Cross[hat, face]];
            {z, rot[face], rot[hat]}];
          flipoutP[p1_, p2_] := Composition[frontflipP[-p2], tacoleftP[-p1]];
      
          (*general function. fit elements of l1 into structure of l2*)
          reshape[l1_, l2_] := Module[{i = 1, length = Length[l1]},
            Map[l1[[Mod[i++, length, 1]]] &, l2, {-1}]];
      
          (*LIFO stack*)
          {pushI, popI} = Module[{stack = {}},
            {(AppendTo[stack, #]; #) &,
             Module[{val = Last[stack]},
               stack = Most[stack];
               val] &}];
      
          With[{vars = First /@ options},
           Module[vars, vars = OptionValue[vars];
      
      
      
           
    2. 
            If[Axiom === None && Rules =!= {}, Axiom = Rules[[1, 1]]];(*default axiom*)
            Axiom = Flatten[{Axiom}];(*normalize to list/directive*)
            {DrawStyle, HatStyle, HatWorldplaneStyle} = Directive /@ {DrawStyle, HatStyle, HatWorldplaneStyle};
      
            Definitions = Join[Definitions, {
               F -> forward, B -> backward, L -> left, R -> right, FO -> flipout, FO[p_] :> flipout[p],
               FF -> frontflip, BF -> backflip, TL -> tacoleft, TR -> tacoright}];
      
            definitionsI = {
              forward[p_] :> forwardP[p], backward[p_] :> forwardP[-p], left[p_] :> leftP[p],
              right[p_] :> leftP[-p], tacoleft[p_] :> tacoleftP[-p], tacoright[p_] :> tacoleftP[p],
              frontflip[p_] :> frontflipP[p], backflip[p_] :> frontflipP[-p], forward -> forwardP[1],
              backward -> forwardP[-1], left -> leftP[Angle], right -> leftP[-Angle], tacoleft -> tacoleftP[-Angle],
              tacoright -> tacoleftP[Angle], frontflip -> frontflipP[Angle], backflip -> frontflipP[-Angle],
              flipout -> flipoutP[Angle, Angle], flipout[p1_] :> flipoutP[p1, Angle],
              flipout[p1_, p2_] :> flipoutP[p1, p2], push -> pushI,
              pop -> Sequence[popI, Identity](*preadjustment for reshape*)};
      
            (*note no memoization. if you try, keep in mind case of RuleDelayed*)
            commands = Nest[Flatten[Replace[#, Rules, {1}]] &, Axiom, Iterations];
            commands = Flatten[((# /. Definitions) /. definitionsI) & /@ commands];
            states = ComposeList[commands, N@{{0, 0, 0}, {0, 1, 0}, {0, 0, 1}}];
      
            points = reshape[First /@ states, Split[popI === # & /@ Join[{0}, commands]]];(*pop is turtle teleportation*)
            points = Composition[First /@ # &, Split] /@ points;(*delete duplicate points*)
      
            Graphics3D[{
              {RandomStuff /. None -> {}, {DrawStyle, Primitive[points]}},
              If[TraceHat,
               hatTrace = {#1, #1 + 2 #3/5} & @@@ states;
               hatTrace = First /@ Split[hatTrace];(*delete duplicate hats*)
      
               hatWorldplane = Polygon[{#1, #2, #4, #3} & @@ Flatten[#, 1]] & /@ Partition[hatTrace, 2, 1];
               {{HatStyle, HatPrimitive[hatTrace]}, {HatWorldplaneStyle, hatWorldplane}}, {}]},
             Quiet@FilterRules[{opts}, Options[Graphics3D]], Boxed -> False]]]]];
      
      Draw[
       {X, push, BF, L, X, R, R, X, pop, R, X, L, TL, L, X, R, F},
       {F -> {F, BF, push, L, X, R, R, X, pop, R, X, L, L, X, R, F},
        X -> {F, BF, push, L, F, R, R, R, F, pop, R, F, L, L, F, R, F}},
       Iterations -> 3, DrawStyle -> {Opacity[.65], Glow[Darker[Red, 2/3]]},
       Definitions -> {X -> Identity}, Angle -> Pi/8]
      

Whoops. Accidentally X-rayed my heart. Or was that this one? In any case, this is the 3D Sierpinski arrowhead curve. It might not look very 3D, but technically it's 3D because it's made out of a tube instead of a line. All joking aside, try as I might I wasn't able to figure out the construction for the 3D arrowhead curve, sadface.

And though this a crushing defeat, we here at the Sierpinski triangle page are stalwart folk for whom such failure is but a rare trigger of recidivistic saccades to our respective vices, for in the characteristic case we amene our fibrile egos by way of the platitudinous homily that what doesn't kill you makes you stronger. In the process of trying to figure out the 3D arrowhead I ended up making an easy-to-use flexible L-system program.

True story, when I woke up this morning I could have sworn my body was contorting into different LOGO curves, in the hope of trial-and-erroring the arrowhead construction. It was like that dream scene in Fight Club, except instead of a girl it was a LOGO curve. Definitely one of the more Freudiologically-awkward memories I'm going to have to carry around for the rest of my life.

  1. proc[img1_, cf_: ColorData[1], mode_: None, blur_: 8] :=
      Module[{img, components, rank, largest, colored},
       img = RemoveAlphaChannel[ColorNegate@ColorConvert[img1, "Grayscale"]];
       components = MorphologicalComponents[img];
    
       Module[{measurements, sorted},
        measurements = ComponentMeasurements[components, "Count"];
        sorted = First /@ Reverse@SortBy[measurements, Last];
        rank[label_] := (rank[label] = Position[sorted, label][[1, 1]])];
    
       colored = Colorize[components, ColorFunction -> (cf[rank[#]] &), ColorFunctionScaling -> False];
       If[mode == "Angelic", colored = ImageMultiply[img, colored]];
       ColorNegate[ImageMultiply[ColorNegate[img], Blur[#, blur] &@ColorNegate[colored]]] // ImageAdjust];
    
    (**)im = Draw[Iterations -> 17, {F -> {B, left[.020944], B}, B -> {L, F}}, RandomStuff -> None,
        Angle -> Pi/5, ImageSize -> 1280, ViewPoint -> {0, 0, Infinity}] // Rasterize;
    
    GradientFilter[im, 5] // ColorNegate // proc[#, Blue &] & // ColorNegate // ImageResize[#, Scaled[1/2]] &
    
    (**)Draw[{A -> {B, L, B}, B -> {A, R, A}}, Primitive -> (Rotate[Line[#], -Pi/24, {0, 0, 1}] &),
     Iterations -> 13, Angle -> 7 Pi/12, Definitions -> {B -> forward, A -> forward},
     DrawStyle -> Opacity[.5], RandomStuff -> None, ViewPoint -> {0, 0, Infinity}]
    
    (**)Draw[Iterations -> 9, {A -> {B, L, B}, B -> {A, R, A}}, Definitions -> {B -> forward, A -> forward},
       RandomStuff -> {Transparent, Sphere[{0, 0, 0}, .05]}, DrawStyle -> {Opacity[.8], Yellow, Glow[Green]},
       ViewPoint -> {0, 0, Infinity}]
    
    (**)d = Draw[{swirl -> ConstantArray[{BF, F, BF, swirl, FO[Pi/12]}, 5]}, DrawStyle -> Opacity[.9], RandomStuff -> None,
        Primitive -> (Line[First@#,
          VertexColors -> (Darker[#, 1/8] & /@ ColorData["AvocadoColors"] /@ Range[0., 1, 1/(Length[First[#]] - 1)])] &),
        Definitions -> {swirl -> backward}, Iterations -> 6, ImageSize -> 2 1280, Method -> {"ShrinkWrap" -> True},
        Background -> Lighter[LightGray, 7/12]] // Rasterize;
    
    d // ImageResize[#, Scaled[1/4]] & // ImageReflect[#, Top -> Bottom] & // ImagePad[#, 2, Lighter[LightGray, 7/12]] &
    
    (**)Draw[Iterations -> 8, {F :> {F, flipout[.2 RandomReal[], Pi RandomReal[]], F}}]
    
    (**)h = Draw[{R -> {B, R, R, R, F}}, Iterations -> 8, Primitive -> Line, RandomStuff -> None,
        Angle -> 1907/2048, ImageSize -> 2 1280, ViewPoint -> {0, 0, Infinity}] // Rasterize;
    
    proc[h // ImageAdjust, Yellow &, "Anglic", 13] // ImageResize[#, Scaled[1/4]] &
    
    (**)diff = ImageDifference @@ Table[
       Draw[{arc, F, arc}, {F -> {F, F, arc, F, arc, F, arc, F, F}}, Primitive -> (Tube[#, .115] &), Angle -> Pi/6,
         Definitions -> {F -> forward[6], arc -> Flatten[Table[{forward[.1], backflip[.899 .1047], right[1/4 .1047]}, {160}]]},
         Iterations -> 2, DrawStyle -> color, RandomStuff -> None, Lighting -> "Neutral", Method -> {"ShrinkWrap" -> True},
         ViewPoint -> {3, -0.25, -1.5}, ViewVertical -> {0.56, -0.66, -0.7}, ImageSize -> 2 1280] // Rasterize,
       {color, {LightGray, White}}];
    
    diff // ColorNegate // ImageAdjust // ImageResize[#, Scaled[1/4]] &
    

What makes this program great is that even just for 2D L-systems, the 3D perspective makes things more intuitive. The arrowhead problem also demanded debugging features such as keeping track of the turtle's orientation, a definite necessity because of the enormous degrees of freedom that geometric L-systems possess.

To give you an idea of this freedom, all of the items in this table are the same exact L-system at the same exact power. The only difference between them is the base angle specified. (By the way, notice Voltron. This is how you know L-systems are Turing complete.) If you take a couple of these to higher powers you get these images (11th and 13th iterations). It's interesting to wonder what some of these might look like at say the thousandth or billionth iteration. Or even, the millionth.

Sidenote. You may have noticed that I never really explained what L-systems are. In fact what I do and don't explain on this page is pretty much completely arbitrary, largely to annoy people who are already familiar with all of this stuff. "Why aren't you mentioning IFS" I hear them crying. Hilarious. But if you've used Mathematica you know that it's well-suited for replacement schemes such as L-systems in a way that is difficult to convey in the context of other languages. Take a look at a simple function definition in Mathematica:

  1. add[a_, b_] := a + b
  2. add[_, _] := 1

What this is saying is: Whenever something matching the pattern add[a_, b_] is found, replace it by a + b. In other words, function application is a special case of pattern matching. Those _ characters are the analogue of the regex . character, the Kleene proton. So a_ means "match any single thing, and call it a". You can in fact do this, which will make the 'function' return 1 when it is passed any two things, as well as use more involved patterns.

I point this out because it can be difficult to appreciate the fundamental straightforwardness of the Mathematica language, I think even for people who have used it for a while. And especially if you're coming to Mathematica from more mainstream languages where the idea of function application being a special case of something more general would be considered some kind of unreachable koan.

The arrowhead isn't the only L-system that can create the Sierpinski figure. More likely there are an infinite number of distinct L-systems that form the Sierpinski triangle in the limit. When we were fiddling with the Sierpinski triangle as a graph, you may have noticed that the zig zag and criss cross had recursive structure:

    1. With[{v = 5},
        axiom = Polygon[{Cos[#], Sin[#]} & /@ (2 Pi Range[v]/v)]];
      
      next[prev_] := prev /. Polygon[pts_] :>
          (Polygon[ScalingTransform[1/2 {1, 1}, #][pts]] & /@ pts);
      
      draw[n_] := Module[{edges},
         edges = Flatten@Nest[next, N@axiom, n] /. Polygon[pts_] :>
            Sequence @@ UndirectedEdge @@@ Partition[pts, 2, 1, 1];
      
         Graph[edges, VertexCoordinates -> VertexList[Graph[edges]],
          VertexSize -> .25]];
      
      draw[2]
      GraphPlot3D[draw[2]]
      
  1. axiom = Polygon[{{0, 0}, {1, Sqrt[3]}/2, {1, 0}}];
    
    next[prev_] := prev /. Polygon[{p1_, p2_, p3_}] :> {
         Polygon[{p1, (p1 + p2)/2, (p1 + p3)/2}],
         Polygon[{p2, (p2 + p3)/2, (p1 + p2)/2}],
         Polygon[{p3, (p1 + p3)/2, (p2 + p3)/2}]};
    
    draw[n_] := Module[{edges},
       edges = Flatten@Nest[next, N@axiom, n] /. Polygon[pts_] :>
          Sequence @@ UndirectedEdge @@@ Partition[pts, 2, 1, 1];
    
       Graph[edges, VertexCoordinates -> VertexList[Graph[edges]],
        VertexSize -> .25]];
    
    g = draw[2];
    cycle = RandomChoice[{FindHamiltonianCycle, FindEulerianCycle}][g][[1]];
    
    Animate[
     HighlightGraph[g, Graph[cycle[[1 ;; n]]],
      EdgeShapeFunction -> (Line[#1] &),
      VertexShapeFunction -> None,
      GraphHighlightStyle -> "DehighlightHide"],
     {n, 1, Length[cycle], 1}, AnimationRate -> 1]
    

We can find these paths for the 3D Sierpinski graph as well, though not necessarily. In fact all along we could have been grapherizing a lot of our stuff, even things like the different distance functions. My point here however is that we may be able to reverse-engineer an L-system from these structures. And it might not actually be hard at all. It does have the down side however of sounding really boring, so on to nonboringer pastures we skidaddle-prance.

Since cellular automata often have the 'world' array joined at the ends, it makes sense to think of their evolution as being on a cylinder:

  1. draw[array_, options___] := Module[
       {interval, topinterval, width, height, f, coords},
       {height, width} = Dimensions[array];
       interval = 2. Pi/width;
       topinterval = 2. Pi (1 + interval)/width;
       coords = Position[array, 1];
    
       f[{x_, r_}] := Rotate[Translate[
          Cuboid[-#, #] &[.5 topinterval {1, 1, 1}],
          {1, 0, -interval x}], interval r, {0, 0, 1}(*;{1,0,0}*), {0, 0, 0}];
    
       Graphics3D[{{Lighter[LightBlue], Opacity[.5],
          Sphere[{0, 0, -interval height/2}, .5]},
         EdgeForm[None], White, f /@ coords}, options, Boxed -> False]];
    
    draw[CellularAutomaton[22,
      ConstantArray[0, 500]~ReplacePart~{1 -> 1, 251 -> 1}, 125],
     Lighting -> "Neutral"]
  2. draw2[im_Image, options___] := draw2[ImageData[ColorConvert[im, "RGB"]], options];
    draw2[array_, options___] := Module[
       {interval, width, height, f, cubes, coords},
       {height, width} = Dimensions[array][[{1, 2}]];
       interval = 2. Pi/width;
       coords = Position[array, p_ /; p != {0, 0, 0}, {2}];
    
       f[{x_, r_}] := Rotate[Translate[
          Cuboid[-#, #] &[.5 interval {1, 1, 1}],
          {1, 0, -interval x}], interval r, {0, 0, 1}, {0, 0, 0}];
    
       cubes = MapThread[{RGBColor @@ #1, f[#2]} &,
         {array[[##]] & @@@ coords, coords}];
    
       Graphics3D[{{Lighter[LightBlue], Opacity[.5],
          Sphere[{0, 0, -interval height/2}, .5]},
         EdgeForm[None], cubes}, options, Boxed -> False]];
    
    (*this rule from "http://web.cecs.pdx.edu/~mm/evca-review.pdf"*)
    rules = Thread[Tuples[{0, 1}, {7}] ->
        IntegerDigits[FromDigits["0504058705000f77037755837bffb77f", 16], 2, 128]];
    
    arr = FixedPointList[CellularAutomaton[rules], RandomInteger[1, 600]];
    arrEdge = ArrayPlot[arr, PixelConstrained -> 1, Frame -> False] // EdgeDetect // ImageData;
    
    (*ad hoc coloring, originally intended for particle animation*)
    pat1 = {{_, _, _, _, _}, {_, 1, 0, 0, 1}, {_, _, _, _, _}};
    pat2 = {{_, 1, _, _, _}, {_, _, 1, _, _}, {_, _, _, 1, _}};
    pat3 = {{_, _, _, _, _}, {_, 1, 1, 1, _}, {_, _, _, _, _}};
    (f[#1 | Reverse /@ #1, _] = #2) & @@@
      {_ -> {0, 0, 0}, pat1 -> {1, 0, 0}, pat2 -> {0, 1, 0}, pat3 -> {0, 0, 1}};
    
    (*see also ImageFilter, ImageConvolve, a million other things*)
    colored = CellularAutomaton[{f, {}, {1, 2}}, arrEdge];
    
    Image[colored]
    draw2[colored, Lighting -> "Neutral"]
    

This is Rule 22 with two initial black squares. It's a cylindrical mapping of this. The sphere in the center is an homage to the Sega Saturn. Long live Sega Saturn, long live Dreamcast. Neo Geo forever. This is a different projection of the same thing, which might actually be easier to comprehend than the cylindrical projection.

And a plot of a range-7 automaton, described in this paper, that was evolutionarily engineered to discriminate between majority-white and majority-black initial conditions. And a particle plot oNEKO!!! Ka-wa-ii. My hope is the image of this dark hieroglyphic cat infests your dreams with nightmares so mindbendingly horrid your perception of reality and fantasy becomes forever warped. Whoops did I say that out loud. See also my Cellular Automata program.

Of course, there are automata whose evolutions are properly three-dimensional, like these quadrilateral versions of Rule 22:

  1. draw[block_, options___] := Graphics3D[
       {EdgeForm[Gray], Cuboid /@ Position[block, 1]},
       options, ViewVertical -> {-1, 0, 0}, Boxed -> False];
    
    draw[CellularAutomaton[{
    115792089237316195423570985008687907853269984665640564039476030751986839257106
      , 2, {1, 1}}, {{{1}}, 0}, 31],
     Lighting -> "Neutral"]
    
  2. (**)
    grids = Partition[#, 3] & /@ Tuples[{1, 0}, {9}];
    rule = IntegerDigits[#, 2, 512] &@
       115792089237316195423570985008687907853269984665640564039476030751986839257106;
    
    Dynamic[FromDigits[rule, 2]]
    Dynamic[draw[CellularAutomaton[{FromDigits[rule, 2], 2, {1, 1}}, {{{1}}, 0}, 31]]]
    
    With[{plot = Function[c, Magnify[ArrayPlot[#1, FrameStyle -> c], 1/6]]},
     Grid[Partition[#, 32], Spacings -> {.1, .1}] &@
      MapIndexed[
       Toggler[Dynamic[rule[[First@#2]]],
         {0 -> plot[LightGray], 1 -> plot[Red]}] &,
       grids]]
    
    (**)
    z = Import["http://upload.wikimedia.org/wikipedia/commons/thumb/e/e0/Game_of_life_glider_gun.svg/610px-Game_of_life_glider_gun.svg.png"];
    z = ImageData[ImageResize[z, Scaled[1/16], Resampling -> "Nearest"] // Binarize // ColorNegate];
    Image[z] // Magnify
    
    With[{f = Switch[
         {#[[2, 2]], Total[#, 2] - #[[2, 2]]},
         {_, 3} | {1, 2}, 1, _, 0] &},
    
     draw[CellularAutomaton[{f, {}, {1, 1}}, {z, 0}, 100]]]
    

An actual 3D automaton whose evolution would be 4-dimensional:

  1. draw[block_, options___] := Graphics3D[
       {EdgeForm[Darker[Gray]], Cuboid /@ Position[block, 1]},
       options, ViewVertical -> {-1, 0, 0}, Boxed -> False];
    
    f[block_, _] := Switch[
       {block[[2, 2, 2]], Total[block, 3] - block[[2, 2, 2]]},
       {_, 4}(*|{1,2}*), 1, _, 0];
    
    evol = CellularAutomaton[{f, {}, {1, 1, 1}},
       {{{{1, 1}, {1, 1}(*,{1,1},{1,1}*)}}, 0}, 15];
    
    ListAnimate[
     draw[#, Lighting -> "Neutral", ImageSize -> 400 {1, 1}] & /@ evol]
    
  2. draw[block_, options___] := Graphics3D[{
        EdgeForm[None],(*Opacity[.8],*)
        Cuboid /@ Position[block, 1],
        Black, Cuboid /@ Position[block, 2]},
       options, Lighting -> "Neutral", Boxed -> False];
    
    f[block_, _] := Switch[
       {block[[2, 2, 2]], Total[block, 3] - block[[2, 2, 2]]},
       {_, 4}, 1, {0, 3}, 2, _, 0];
    
    evol = CellularAutomaton[{f, {}, {1, 1, 1}},
       {CrossMatrix[1 {1, 1, 1}]~BitXor~1, 0}, 25];
    
    (*can be flashy*)
    (*ListAnimate[draw[#,ViewPoint->Top,ImageSize->400 {1,1}]&/@evol]*)
    
    draw[Last[evol],
       ImageSize -> 2 1280, ViewPoint -> 2 {1, 1, 1},
       Lighting -> {{"Point", Yellow, Scaled[{1, 1, 1}], 5}},
       Method -> {"ShrinkWrap" -> True}] //
      Rasterize // ImageResize[#, Scaled[1/4]] &

And just so we're all clear, time isn't "the fourth dimension." That statement is the conceptual version of eating bagels without cream cheese, namely a manifestation of meaniglessness.

In rectangular 3D each cell is surrounded by $3^3 - 1 = 26$ cells, so the number of even just simple totalistic rules is very large, nevermind starting configurations. This means that finding "interesting" rules and configurations can be a tricky artform. This is another place where I could use, say, a warehouse full of Alienware laptops.

If you have Mathematica 9 (must be nice), its Image3D functionality is perfect for these 3Dified cellular automata. And speaking of grid thingies, let's not forget our unexpectedly-glorious matrix replacement scheme:

  1. 
    
    
    
    
    
    (**)
    Begin["mmx`"];
    
    matrixInput3D1[Dynamic[tensor_], Dynamic[color_], options___] :=
      Dynamic@Module[{grid},
        grid = Position[ArrayPad[tensor, {0, -1}], _?IntegerQ];
    
        EventHandler[#, {"MouseDown", 2} :> {}] &@
           Graphics3D[{#, Transparent, EdgeForm[LightGray], Cuboid /@ grid},
            options,(*Method->{"ShrinkWrap"->True},*)Boxed -> False] &@
    
         Array[With[{loc := tensor[[##]]},
            Mouseover[
             (**){Style[#, Darker[color, .65]] &@
               Text[Dynamic[loc /. 0 -> Style[0, Opacity[.5]]], {##}],
              Opacity[loc /. {0 -> .1, 1 -> .3}], Sphere[{##}, .2]},
             (**){Text[EventHandler[Checkbox[Dynamic[loc], {0, 1}],
                {"MouseDown", 2} :> (loc = 0)], {##}],
              Opacity[.01], Sphere[{##}, .2]}]] &,
          Dimensions[tensor]]];
    
    matrixInput3D2[Dynamic[tensor_], Dynamic[rules_], Dynamic[color_], options___] :=
      Dynamic@DynamicModule[{grid},
        grid = Flatten[Array[List, Dimensions[ArrayPad[tensor, {0, -1}]]], 2];
    
        EventHandler[#, {"MouseDown", 2} :> {}] &@
           Graphics3D[{#, Transparent, EdgeForm[LightGray], Cuboid /@ grid},
            options,(*Method->{"ShrinkWrap"->True},*)Boxed -> False] &@
    
         Array[With[{loc := tensor[[##]]},
            With[{display = Tooltip[Panel[#, FrameMargins -> None],
                 Column[{loc /. rules /. {Reverse -> "R", Transpose -> "T",
                     Composition -> List, Verbatim[Slot][_] :> "m"},
                   "", "Click to cycle", "Right-click to zero"}],
                 TooltipDelay -> .6] &},
    
             Mouseover[
              (**){Style[#, Darker[color, .65]] &@
                Text[Dynamic[loc /. 0 -> Style[0, Opacity[.5]]], {##}],
               Opacity[loc /. {0 -> .1, _ -> .3}], Sphere[{##}, .2]},
              (**){Text[EventHandler[
                 display[
                  Toggler(*PopupMenu*)[Dynamic[loc], First /@ rules,
                   ImageSize -> Automatic]
                  ],
                 {"MouseDown", 2} :> (loc = 0)], {##}],
               Opacity[.01], Sphere[{##}, .2]}]]] &,
          Dimensions[tensor]]];
    
    bg = White;
    dims = # -> If[# > 2, Style[#, Red], #] & /@ Range[5];
    
    rotations = Flatten@Outer[Function[{o, dir},
         Composition[Transpose[#, o] &, dir /@ # &, Transpose[#, o] &]],
        {{1, 2, 3}, {3, 2, 1}, {2, 1, 3}},
        {Composition[Transpose, Reverse],
         Composition[Reverse, Transpose],
         Reverse, Transpose}, 1];
    
    rotations = MapIndexed["S" @@ #2 -> #1 &, rotations];
    defaultRules = Join[{0 -> (0 # &), 1 -> (# &)}, rotations];
    
    iterate[matrix0_, matrixT_, rules_, power_] :=
      Nest[Function[prev,
        ArrayFlatten[Map[#[prev] &,
          Replace[matrixT, rules, {3}], {3}], 3]],
       matrix0, power];
    
    randomMatrix[dimensions_, source_] := With[
       {rv := RandomVariate[ZipfDistribution[Length[source], 1]]},
       Array[source[[rv]] &, dimensions]];
    
    With[{HiPrint := Function[viewpoint,
        With[{pow = power},
         CellPrint[ExpressionCell[
           Defer[
            powzerz = pow;
            With[{objects = Translate[primitive,
                Replace[Position[iterate[
                   matrix0 /. 0 matrix0 -> {{{1}}},
                   matrixT /. 0 matrixT -> {{{1}}},
                   rules, powzerz], If[negativeSpace, 0, 1]],
                 {} -> {1, 1, 1}]]},
             ImageResize[Rasterize[#], Scaled[1/4]] &@
              Defer[Graphics3D][{color, Opacity[opacity],
                Glow[glow], Specularity[specularity],
                EdgeForm[{Opacity[opacity], Darker[color, 4 .15]}], objects},
               Lighting -> "Neutral", Method -> {"ShrinkWrap" -> True},
               ImageSize -> {Automatic, 4 732}, Boxed -> False,
               ViewPoint -> viewpoint, ViewVertical -> vv,
               Background -> background]]],
           "Input"]]]],
    
      printMatrices := Function[
        CellPrint[ExpressionCell[DynamicModule[{
            mtx0 = matrix0, mtxT = matrixT, mtx0o = matrix0, mtxTo = matrixT,
            clr = color, opc = opacity, ns = negativeSpace, pow = power, rls = rules,
            prm = primitive, iter = iterate, bg = background, vp1 = vp, vv1 = vv},
    
           With[{
             btn = Button[DynamicWrapper["print data",
    
                If[mtx0 =!= mtx0o || mtxT =!= mtxTo, mtx0 = mtx0o; mtxT = mtxTo]],
               Print[Grid[{
                  {"kernel matrix", MatrixForm[mtx0o]},
                  {"transformation matrix", MatrixForm[mtxTo]},
                  {"rules", rls}, {"power", pow}}]]],
             mtx0c = matrixInput3D1[Dynamic[mtx0], Dynamic[clr],
               SphericalRegion -> True, ImageSize -> Small,
               Background -> Lighter[bg, .8],
               ViewPoint -> Dynamic[vp1], ViewVertical -> Dynamic[vv1]],
             mtxTc = matrixInput3D2[Dynamic[mtxT], Dynamic[rls], Dynamic[clr],
               SphericalRegion -> True, ImageSize -> Small,
               Background -> Lighter[bg, .8],
               ViewPoint -> Dynamic[vp1], ViewVertical -> Dynamic[vv1]],
             g3d = With[{objects = Translate[prm,
                  Replace[Position[iter[
                     mtx0 /. 0 mtx0 -> {{{1}}},
                     mtxT /. 0 mtxT -> {{{1}}},
                     rls, pow], If[ns, 0, 1]],
                   {} -> {1, 1, 1}]]},
               Graphics3D[{
                 EdgeForm[{Opacity[opc], Darker[clr, 4 .15]}],
                 clr, Opacity[opc], objects},
                ImageSize -> Small, Boxed -> False, SphericalRegion -> True,
                ViewPoint -> Dynamic[vp1], ViewVertical -> Dynamic[vv1],
                Lighting -> "Neutral", Background -> bg]]},
    
            Panel[Grid[{
               {Panel[Placeholder["name"]], SpanFromLeft, btn},
               {mtx0c, mtxTc, g3d}}]]]]]]],
    
      (* controls *)
      dim0C = Control[{{dim0, 1, ""}, dims, ControlType -> PopupMenu}],
      dimTC = Control[{{dimT, 2, ""}, dims, ControlType -> PopupMenu}],
      matrix0C = matrixInput3D1[Dynamic[matrix0], Dynamic[color],
        SphericalRegion -> True, ImageSize -> Dynamic[imgSize1],
        Background -> Dynamic[Lighter[background, .8]],
        ViewPoint -> Dynamic[vp], ViewVertical -> Dynamic[vv]],
      matrixTC = matrixInput3D2[Dynamic[matrixT], Dynamic[rules], Dynamic[color],
        SphericalRegion -> True, ImageSize -> Dynamic[imgSize2],
        Background -> Dynamic[Lighter[background, .8]],
        ViewPoint -> Dynamic[vp], ViewVertical -> Dynamic[vv]],
      rulesC = Pane[Style[#, 10], {400, 200}, Scrollbars -> Automatic] &@
        Control[{{rules, defaultRules, ""},
          InputField, Background -> Dynamic[Lighter[background, .65]],
          FieldSize -> {50, {0., Infinity}}}],
      colorC =
       Control[{{color, RGBColor[.15, .6, 1], "color"}, ColorSlider}],
      backgroundC = Row[{"background   ", Framed[
          ColorSlider[Dynamic[background, (bg = background = #) &],
           AppearanceElements -> "Swatch"],
          FrameStyle -> Gray], " ",
         ColorSlider[Dynamic[background, (bg = background = #) &],
          AppearanceElements -> "Spectrum", ImageSize -> Small]}],
      opacityC = Control@{{opacity, 1, "opacity"}, 0, 1, ImageSize -> Small},
      glowC = Control[{{glow, Black, "glow"}, ColorSlider}],
      specC = Control[{{specularity, Black, "specularity"}, ColorSlider, ImageSize -> Small}],
      primC = Control[{{primitive, Scale[Cuboid[],.99999], "primitive"},
         # -> Graphics3D[{color, #}, Boxed -> False, ImageSize -> 20] & /@
          {{PointSize[0], Point[{0., 0., 0.}]}, Sphere[{0., 0., 0.}, .5],
           {EdgeForm[None], Scale[Cuboid[],.99999]}, Scale[Cuboid[],.99999]}, SetterBar}],
      powerC = Control[{{power, 1, "power"}, 0, 5, 1, Appearance -> "Labeled"}],
      nsC = Control[{{negativeSpace, False,
          Tooltip["negative", "negative space",
           TooltipDelay -> .4]}, {False, True}}]
      },
    
     (*control layout*)
     With[{controls :=
        Row[{
          Column[{
            Row[{dim0C, "   |", dimTC}],
            Row[{"    ", matrix0C, "  ", matrixTC}]}], Spacer[40],
          Column[{
            OpenerView[{"Rules", rulesC}],
            OpenerView[{"Style",
              Column[{
                Row[{
                  Column[{colorC, backgroundC}], Spacer[40],
                  Column[{glowC, specC}]}],
                Row[{opacityC, Spacer[20], nsC, Spacer[20], primC}]}]}],
            powerC}]}],
    
       bookmarks := {
         Overscript["Random kernel matrix", ""] :>
           (matrix0 = randomMatrix[Dimensions[matrix0], {0, 1}]),
         "Random transformation matrix" :>
           (matrixT = randomMatrix[Dimensions[matrixT], First /@ defaultRules]),
         "Random both" :> (
           matrix0 = randomMatrix[Dimensions[matrix0], {0, 1}];
           matrixT = randomMatrix[Dimensions[matrixT], First /@ defaultRules]),
    
         Overscript["Clear kernel matrix", ""] :> (matrix0 = 0 matrix0),
         "Clear transformation matrix" :> (matrixT = 0 matrixT),
         "Clear both" :> ({matrix0, matrixT} = 0 {matrix0, matrixT}),
    
         Overscript["Invert kernel matrix", ""] :> (matrix0 = BitXor[matrix0, 1]),
         "Invert transformation matrix" :> (matrixT = Replace[matrixT, {0 -> 1, _ -> 0}, {3}]),
    
         Overscript["Print matrices", ""] :> printMatrices[],
    
         Overscript["HiPrint", ""] :> HiPrint[vp],
         "HiPrint Far" :> HiPrint[1000 vp]}},
    
      Panel[#, Background -> Dynamic[bg]] &@
       Manipulate[Module[{g3d, side},
    
         If[dim0 {1, 1, 1} =!= Dimensions[matrix0], matrix0 = PadRight[matrix0, dim0 {1, 1, 1}]];
         If[dimT {1, 1, 1} =!= Dimensions[matrixT], matrixT = PadRight[matrixT, dimT {1, 1, 1}]];
         If[bg =!= background, bg = background];
    
         Module[{matrixP},(*remove rules from matrix that no longer exist*)
          matrixP = Map[Function[a, If[a === Replace[a, rules], rules[[1, 1]], a]], matrixT, {3}];
          If[matrixT =!= matrixP, matrixT = matrixP]];
    
         g3d = With[{objects = Translate[primitive,
              Replace[Position[iterate[
                 matrix0 /. 0 matrix0 -> {{{1}}},
                 matrixT /. 0 matrixT -> {{{1}}},
                 rules, power], If[negativeSpace, 0, 1]],
               {} -> {1, 1, 1}]]},
           Graphics3D[{
             Dynamic[EdgeForm[{Opacity[opacity], Darker[color, 4 .15]}]],
             Dynamic[color], Dynamic[Opacity[opacity]], Dynamic[Glow[glow]],
             Dynamic[Specularity[specularity]], objects},
            ImageSize -> {{300, Large}, {300, Large}},
            Lighting -> "Neutral", Background -> Dynamic[background]]];
    
         side = Map[Function[vp1,
            Tooltip[#, ViewPoint -> vp1, TooltipDelay -> .3] &@
    
               EventHandler[#,
                "MouseDown" :> (vp = vp1 /. Infinity -> 4; vv = {0, 0, 1})] &@
             Framed[Deploy[
               Show[g3d, ViewPoint -> vp1, ImageSize -> Small, Boxed -> False]],
              FrameStyle -> Gray, Background -> Dynamic[background]]],
           Permutations[{Infinity, 0, 0}]];
    
         Row[{Column[side,(*Dividers->All,*)FrameStyle -> Gray],
           Show[g3d, Boxed -> False, SphericalRegion -> True,
            (*PlotRangePadding->.001,*)
            ViewPoint -> Dynamic[vp], ViewVertical -> Dynamic[vv]]}]
         ],
    
        {{vv, {0, 0, 1}}, ControlType -> None},
        {{vp, {1.3, -2.4, 2}}, ControlType -> None},
        {{imgSize1, Small},
         ControlType ->
          None},(*prevent matrix controls from autoresizing*)
        {{imgSize2, Small}, ControlType -> None},
        {{background, White}, ControlType -> None},
        {{matrix0,
          If[dim0 < 2, {{{1}}}, randomMatrix[dim0 {1, 1, 1}, {0, 1}]]},
         ControlType -> None},
        {{matrixT,
          If[dimT < 2, {{{1}}},
           randomMatrix[dimT {1, 1, 1}, First /@ defaultRules]]},
         ControlType -> None},
        controls, Bookmarks :> bookmarks,
        LabelStyle -> Darker[Gray], SynchronousUpdating -> Automatic,
        Paneled -> False, SaveDefinitions -> True, Alignment -> Center]]]
    
    (**)
    End[];
    
    
    
    
    
    
    

This scheme clearly shows the projective character of these algorithms. Take for example this nifty 3D plus sign made of 3D plus signs, holy mathphobia inducer. It looks like a 2D fractal plus sign when viewed along each axis, but resembles various 2D constructions when viewed from mixed angles.

What's not obvious from these images is that the matrix controls at the top (the pink spheres) and the output figure share the same viewpoint (twirl one, the other two follow). In Mathematica this is as easy as wrapping a couple of things in Dynamic[ ], after which the system takes care of automatically updating things as necessary. It's pretty much the ideal of what event handling should be, at least for these kinds of applications. The underlying engineering for this on Mathematica's part must be very intricate.

And speaking of intricate, this is probably the most complicated Mathematica program I've so far written, in part because I didn't run it through any last-phase refactoring. If you have the courage to fiddle with this program (and I encourage you to have this courage, as the program has a particular issue I couldn't solve), be prepared to suffer dearly for my laziness.


Give me a moment.


OK, it looks like we're in the inversion section. Where did all this 3D stuff come from? Holy cow. HOLY BRAHMAN DATA COW. Oh I think I know which voice it was. Irregardless, since a bunch of 3D things essentially just programmed themselves into existence while I wasn't looking, this means we can do 3D INVERSIONS!!!! Chaos game.

  1. draw[vertices_, numPts_] :=
      Graphics3D[{PointSize[0], Opacity[.1],
        Point[FoldList[(#1 + #2)/2 &, .5 First[vertices],
          RandomChoice[N@vertices, numPts]]]},
       Boxed -> False];
    
    invert[p_] := p/Norm[p]^2;
    
    vertices = PolyhedronData[{"Pyramid", 3}, "VertexCoordinates"];
    vertices = Normalize /@ (# - Mean[vertices] &) /@ vertices;
    
    Show[
     draw[vertices, 20000],
     draw[vertices, 100000] /. Point[pts_] :> Point[invert /@ pts]]
    

This. Four-headed tri-jawed infinity-mouthed Pac-man langolier. If the world ever decides to give me a nightmare, I hope it picks one of these adorable things to chase me through the dark recesses of my deranged mind. Geometric.

  1. draw[shape_, n_] := Module[{next},
       (*scale by 1/2 toward each vertex, in turn*)
       next[prev_] := Scale[prev, 1/2, #] & /@ shape[[1]];
    
       Graphics3D[{EdgeForm[Opacity[.15]],
         Nest[next, N@shape, n]},
        Lighting -> "Neutral", Boxed -> False]];
    
    invert[p_] := p/Norm[p]^2;
    
    shape = PolyhedronData[{"Pyramid", 3}, "Faces"];
    shape[[1]] = Normalize /@ (# - Mean[shape[[1]]] &) /@ shape[[1]];
    
    Show[
     draw[shape, 3],
     (draw[shape, 4] // Normal) /.
      Polygon[pts_, __] :> Polygon[invert /@ pts]]
    

The ostensive architectonics, quite awesome. c.f. Dyson sphere. The code however is simple. Cobra.

  1. draw[shape_, n_] := Module[{next},
       (*scale by 1/2 toward each vertex, in turn*)
       next[prev_] := Scale[prev, 1/2, #] & /@ shape[[1]];
    
       Graphics3D[{EdgeForm[Opacity[.15]],
         Opacity[.75], Black, Nest[next, N@shape, n]},
        Lighting -> "Neutral", Boxed -> False]];
    
    transform[1][p_] := p^3/Norm[p]^2;
    transform[2][p_] := (Reverse[p].p) p/Norm[p]^2;
    transform[3][p_] := (Reverse[p].Cross[{0, 0, 1}, p]) p/Norm[p]^2;
    
    shape = PolyhedronData[{"Pyramid", 3}, "Faces"];
    
    (draw[shape, 4] // Normal) /.
      Polygon[pts_, __] :> Polygon[transform[1] /@ pts]
    

And fishie! Logarithmic.

  1. game = Compile[{{vertices, _Real, 2}, {w, _Real}, {numpoints, _Integer}},
       Module[{diff},
        FoldList[(diff = #2 - #1;
           Clip[(#1 + #2) Log[Sqrt[diff.diff] + w]]) &,
         {0, 0, 0}, RandomChoice[vertices, numpoints]]]];
    
    invert[p_ /; Norm[p] < .25] := 4 Normalize[p];
    invert[p_] := p/Norm[p]^2;
    
    vertices = PolyhedronData[{"Pyramid", 3}, "VertexCoordinates"];
    (*vertices=Normalize/@(#-Mean[vertices]&)/@vertices;*)
    
    Graphics3D[{PointSize[0], Opacity[.2],
      Point[invert /@ game[vertices, .01, 400000]]},
     Boxed -> False]
    

"Chaos game with logarithmic distance function" is a bit long. We need to give this specific kind of fractal a name. What about "Charlie render"? So I'd be like "here we have an inverted Charlie render at $w$ factor .01" and people would nod comprehendingly while reading that, as if there were an established literature on Charlie renders.

You might object that the contours of this nomenclature don't quite align with the striking yet oft- hauntingly quiescent leylines of its intended referents, but you would be wrong — the matching is nigh onomatopoeial per my linguistic auteurity. Incidentally, you should see what my writing looks like when I really cut loose. Rejoice asplendent my sparing you that paragon 'cross the rubicon, padawan.

Since the originals have a lot of points close to 0, their inverses have a lot of points at very large distances. In this case I've decided to clamp the maximum distance of points to a short range (essentially putting them on a leash, like those ball & chain dogs in Mario Bros. 3). It's another way of dealing with infinities. I like this approach because it preserves the radial texture of the figure, snowglobe-like. Taking this to its conclusion, we normalize all points to the same distance:

  1. game = Compile[{{vertices, _Real, 2}, {w, _Real}, {numpoints, _Integer}},
       Module[{diff},
        FoldList[(diff = #2 - #1;
           Clip[(#1 + #2) Log[Sqrt[diff.diff] + w]]) &,
         {0, 0, 0}, RandomChoice[vertices, numpoints]]]];
    
    vertices = PolyhedronData[{"Pyramid", 3}, "VertexCoordinates"];
    (*vertices=Normalize/@(#-Mean[vertices]&)/@vertices;*)
    
    Module[{pts},
     pts = game[vertices, .01, 400000];
    
     Graphics3D[{
       {Glow[White], Sphere[{0, 0, 0}, .99999]}, PointSize[0], Opacity[.5],
       Point[Normalize /@ pts,
        VertexColors -> (ColorData["AvocadoColors"] /@ Norm /@ pts)]},
      ViewPoint -> {Sqrt[3], -Sqrt[8], 1}, Boxed -> False]]
    

These two are the same, except the first one has an opaque sphere in the interior so that you can't see points beyond the horizon. The extra points in the second one are on the other side of the globe. These points are colored according to their original distance. And the unnormalized figure.


Questions

How many points does the Sierpinski triangle have, besides infinity? Say at a given iteration?

  1. axiom = Polygon[{{0, 0}, {1, Sqrt[3]}/2, {1, 0}}];
    
    next[prev_] := prev /. Polygon[{p1_, p2_, p3_}] :> {
         Polygon[{p1, (p1 + p2)/2, (p1 + p3)/2}],
         Polygon[{p2, (p2 + p3)/2, (p1 + p2)/2}],
         Polygon[{p3, (p1 + p3)/2, (p2 + p3)/2}]};
    
    draw[n_] := Module[{g},
       g = Graphics[{
          White, EdgeForm[Black],
          Nest[next, axiom, n]}];
    
       Show[
        (g // next) /. p : Polygon[pts_] :>
          {p, Black, Disk[#, (1/2)^(n + 5)] & /@ pts},
        g /. Polygon[pts_] :>
          {EdgeForm[None], Disk[#, (1/2)^(n + 5.075)] & /@ pts}]];
    
  2. FindSequenceFunction[
      Length@DeleteDuplicates@
          Cases[draw[#], Disk[p_, ___] :> p, Infinity] & /@ Range[6]
      ][n]

All additions are in powers of 3. So at a given iteration we have $\Sigma \space 3^k$ total points. There's all sorts of ways to find the closed form of this sum, not least of which is to use the internet. I'm a fan of the algebraic approach:

$$ \begin{array}{ccccccccccc} S & = & 3^1 & + & 3^2 & + & 3^3 & + & \cdots & + & 3^n \\ -3S & = & & - & 3^2 & - & 3^3 & - & \cdots & - & 3^{n+1} \\ -2S & = & 3^1 & - & 3^{n+1} \\ S & = & \frac 1 2 (3^{n+1}-3) \\ S & = & \frac 3 2 (3^n-1) \\ \end{array} $$

The nice thing about this kind of manual deduction is that it gives us an excuse to plaster more math on our page, giving perusers who don't know any better the impression that we're really smart. This sum accounts for the additions. We also need to account for the first 3 points. For a given iteration, we have a total of $\frac 3 2 (3^n-1) + 3=\frac 3 2 (3^n + 1)$ points. The arithmetic works out better if we count the polygons instead of the points.

If my web search kune do hasn't failed me, this would make most of our algorithms "geometric space and therefore time" (GSATT) algorithms. Actually I just made that up, I don't know what they're called. It's not really relevant for us since the geometricness also means we get a large number of points with few iterations.

What does the "integration" of the Sierpinski triangle look like? There's various ways to interpret this in 2D, but I'm curious about how the number of points of the triangle increases along a straight line, as if the triangle were a single-variable function:

    1. next[prev_] := prev /. Interval[{a_, b_}] :> {
           Interval[{a, a + (b - a)/3}],
           Interval[{a + 2 (b - a)/3, b}]};
      
      cantor[n_] := IntervalUnion @@ Flatten@
          Nest[next, N@Interval[{0, 1}], n];
      
      rectangles[n_, h_: .02, scale_: 1] :=
        Nest[next, Interval[{0, 1}], n] /. Interval[{a_, b_}] :>
          Rectangle[{a, -h (n + 10 h) scale}, {b, -h (n + 1) scale}];
      
      (*this "integration" depends on the "curve" being "uniformly sampled"*)
      int[pts_] := MapIndexed[{##} /.
           {{x_, y_}, {i_}} :> {x, i} &, SortBy[pts, First]];
      
      set = cantor[16];(*this is 2^16 intervals*)
      {null, {pts}} = Reap[Do[
         If[IntervalMemberQ[set, a], Sow[{a, 0}]],
         {a, 0., 1., 1/1000000}]];
      
      Graphics[rectangles /@ Range[6]]
      
      Show[Graphics[rectangles[6] /.
         Rectangle[{x1_, y1_}, {x2_, y2_}] :> Rectangle[{x1, 0}, {x2, .02}]],
       ListLinePlot[{#1, #2/Length[pts]} & @@@ int[pts],
        PlotStyle -> Black]]
      
  1. axiom = Polygon[{{0, 0}, {1, 0}, {1, 1}}];
    
    next[prev_] := prev /. Polygon[{p1_, p2_, p3_}] :> {
         Polygon[{p1, (p1 + p2)/2, (p1 + p3)/2}],
         Polygon[{p2, (p2 + p3)/2, (p1 + p2)/2}],
         Polygon[{p3, (p1 + p3)/2, (p2 + p3)/2}]};
    
    points[n_] := DeleteDuplicates[Flatten[
        Nest[next, N@axiom, n] /. Polygon -> Sequence, n]];
    
    (*this "integration" depends on the "curve" being "uniformly sampled"*)
    int[pts_] := MapIndexed[{##} /.
         {{x_, y_}, {i_}} :> {x, i} &, SortBy[pts, First]];
    
    pts = points[10];
    
    Show[
     Graphics[{Opacity[.1], PointSize[0], Black, Point[pts]}],
     ListLinePlot[{#1, #2/Length[pts]} & @@@ int[pts], PlotStyle -> Black]]
    

Hmm. I was hoping it would look something like the so-called Devil's Stairscase, which is the same thing for the Cantor set. You can just feel the Staircase's ragged darkness filling you with joy. But this, this looks like the underside of a fluffy cloud. I think I will call it Lumpy Space Satan's Hairline. Not as dark and morally grimy a name as I was hoping to coin, but not bad either.

My original reason for inverting the Sierpinksi triangle was to see how it might magnify the inner texture. I.e. turning the triangle inside out to make the inside more visible. "You could have explained that in the actual inversion section" you say. Indeed, but let's not hark on couldas and shouldas. The point is there is an intuition behind these things, and we can ask other questions in the same spirit. For example, what if we extend the 2D Sierpinski triangle into 3D, with each point a different $z$ coordinate (depth) depending on its distance from the center of the triangle?

  1. draw[v_, n_] := Module[{ring, figure},
       ring[c_, r_, depth_] := Module[{ps},
         ps = c + r {Cos[#], Sin[#]} & /@ (2 Pi Range[v]/v);
    
         If[depth == 0, Polygon[ps],
          ring[(c + #)/2, r/2, depth - 1] & /@ ps]];
    
       figure = ring[0., 1., n] /. Polygon[pts_] :>
          Polygon[{#1, #2, Norm[{#1, #2}]} & @@@ pts];
    
       (*figure=ring[0.,1.,n]/.Polygon[pts_]:>
       Polygon[Normalize[#]~Append~Norm[#]&/@pts];*)
    
       Graphics3D[{Transparent,
         EdgeForm[{Opacity[.5], Black}], figure}]];
    
    draw[3, 5]
    

We get what we expect, a boomerang-looking thing. And look at this lovely demonic-looking Moire pattern, surely the universe's recompense for that fluffy cloud nonsense above. We can also normalize the points so that all we see is the radial detail. That produces a coronet-looking thing, which can be unrolled:

  1. draw[v_, n_, s_: 2, cutoff_: 0, width_: 1] := Module[{ring, figure},
       ring[c_, r_, depth_] := Module[{ps},
         ps = c + r {Cos[#], Sin[#]} & /@ (2 Pi Range[v]/v);
    
         If[depth == 0, Polygon[ps],
          ring[(c + #)/2, r/2, depth - 1] & /@ ps]];
    
       figure = ring[0., 1., n] /. Polygon[pts_] :>
          Polygon[{ArcTan @@ (# /. {0., 0.} -> {1., 0.}), s Norm[#]} & /@ pts];
       figure = Flatten[figure];
    
       figure = Cases[figure, Polygon[pts_] /;
          Mean[Norm /@ Differences[pts]] < .5
           (*&&MemberQ[First/@pts,a_/;-width Pi/v<a<width Pi/v]*)
           && MemberQ[Last /@ pts, y_ /; y > cutoff s]];
    
       Graphics[{Opacity[.5],
         EdgeForm[{Opacity[.13], JoinForm["Round"]}],
         figure}, ImageSize -> Large]];
    
    draw[6, 4, 2, .4] /. Polygon[pts_] :>
      {Opacity[.5], EdgeForm[{Opacity[.01], LightGray}],
       Hue[.05 Norm[Mean[pts]]], Polygon[pts]}
    

What does a radial histogram of the Sierpinski triangle look like?

  1. axiom = Polygon[{Cos[#], Sin[#]} & /@ (Pi/2 + 2 Pi Range[3]/3)];
    
    next[prev_] := prev /. Polygon[{p1_, p2_, p3_}] :> {
         Polygon[{p1, (p1 + p2)/2, (p1 + p3)/2}],
         Polygon[{p2, (p2 + p3)/2, (p1 + p2)/2}],
         Polygon[{p3, (p1 + p3)/2, (p2 + p3)/2}]};
    
    points[n_] := DeleteDuplicates[Flatten[
        Nest[next, N@axiom, n] /. Polygon -> Sequence, n]];
    
    pts = points[8];
    stats = Transpose@MapAt[Partition[#, 2, 1] &, #, 1] &@
       HistogramList[ArcTan @@@ pts, "Knuth"];
    
    max = Max[Last /@ stats];
    polys = Polygon[#2/max {{0, 0},
           {Cos[#1[[1]] + .005], Sin[#1[[1]] + .005]},
           {Cos[#1[[2]]], Sin[#1[[2]]]}}] & @@@ stats;
    (*poly=Polygon[#2/max{Cos[#1[[1]]],Sin[#1[[1]]]}&@@@stats];*)
    
    Graphics[{
      {PointSize[0], Opacity[.1], Point[pts]},
      {ColorData[1][1], polys}}]
    

What happens if we run the Game of Life on the Sierpinski triangle?

  1. sier[n_] := Mod[Array[Binomial, {n, n}, 0], 2];
    
    s = ArrayPad[sier[2^9], 2^5];
    (*s=sier[2^9]+Transpose[sier[2^9]]/. 2->1;*)
    i = 0;
    
    PrintTemporary[Dynamic[i]];
    PrintTemporary[Dynamic[Image[s]]];
    
    With[{lifeSpec = {224, {2, {{2, 2, 2}, {2, 1, 2}, {2, 2, 2}}}, {1, 1}}},
      While[i++ < Infinity,
       s = CellularAutomaton[lifeSpec, s]]];
    
    i
    Image[s]
    

Basically nothing. The triangle does this and that, shoots a couple gliders, settles. Larger versions do more or less the same thing but take longer to settle. Not very interesting, but it raises the idea of using fractals as starting configurations. However, on the internet I found that lines produce Sierpinski triangles:

  1. Export["c:/users/zrp/desktop/line.bmp",
     Image[{ConstantArray[0, 2^13]}]]
  2. With[{lifeSpec = {224, {2, {{2, 2, 2}, {2, 1, 2}, {2, 2, 2}}}, {1, 1}}},
      frames = CellularAutomaton[lifeSpec,
        {{ConstantArray[1, 2^8(*-14*)]}, 0}, 130]];
    
    Export["c:/users/zrp/desktop/zrp.gif",
     ColorNegate /@ Image /@ frames, "DisplayDurations" -> .17]
    
  3. With[{lifeSpec = {224, {2, {{2, 2, 2}, {2, 1, 2}, {2, 2, 2}}}, {1, 1}}},
      frames = CellularAutomaton[lifeSpec,
        Boole@Array[#1 == 2^7 &, 2^8 {1, 1}], 2^7]];
    
    Export["c:/users/zrp/desktop/zrp.gif",
     ColorNegate /@ Image /@ frames, "DisplayDurations" -> .17]

I didn't even have to add the horns. This is one end of a line after some iterations. The pattern continues propagating forever and ever as long as there is line left and becomes more distinguished at larger scales. It appears to be driven entirely by the line itself. Consider the evolution of a line that is infinitely long, something you can actually witness in the Game of Life by connecting the edges of the board.

As the finite line splits, it leaves debris due to the circumstances of the ends. The pattern you end up with is a trace of the line's subdivisions. It's because the line splits cleanly and does so in a Sierpinski recursion that you end up with clear Sierpinski triangles at larger scales.

If you want to play with large Game of Life constructions, the easiest way is to export them as images and open them in a dedicated Game of Life program, as those can run the game at very high speeds.

What does a random walk on the Sierpinski graph look like?

  1. axiom = Polygon[{{0, 0}, {1, Sqrt[3]}/2, {1, 0}}];
    
    next[prev_] := prev /. Polygon[{p1_, p2_, p3_}] :> {
         Polygon[{p1, (p1 + p2)/2, (p1 + p3)/2}],
         Polygon[{p2, (p2 + p3)/2, (p1 + p2)/2}],
         Polygon[{p3, (p1 + p3)/2, (p2 + p3)/2}]};
    
    draw[n_] := Module[{edges},
       edges = Flatten@Nest[next, N@axiom, n] /. Polygon[pts_] :>
          Sequence @@ UndirectedEdge @@@ Partition[pts, 2, 1, 1];
    
       Graph[edges, VertexCoordinates -> VertexList[Graph[edges]],
        VertexSize -> .25]];
    
    graphWalkPath[g_Graph, steps_: 15] := Module[{neighbors},
       neighbors[vertex_] := (neighbors[vertex] =
          Complement[VertexList[NeighborhoodGraph[g, vertex]], {vertex}]);
    
       NestList[RandomChoice[neighbors[#]] &, First[VertexList[g]], steps]];
    
    SetAttributes[UndirectedEdge, Orderless];
    graphWalk[args__] := Graph[DeleteDuplicates[
        UndirectedEdge @@@ Partition[graphWalkPath[args], 2, 1]]];
    
    g = draw[3];
    Grid[Partition[#, 10]] &@Table[
      Graphics[{Opacity[.8], JoinForm["Round"], Line[graphWalkPath[g, 50]]},
       ImageSize -> 50 {1, 1}], {100}]
    
    Grid[Partition[#, 10]] &@Table[
      HighlightGraph[g, graphWalk[g, 50],
       ImageSize -> 50 {1, 1}], {100}]
    

About what you would expect. I'll leave the stats to those whose laziness is bounded from above, instead of below. What does a "circle" look like on the Sierpinski graph?

  1. axiom = Polygon[{{0, 0}, {1, Sqrt[3]}/2, {1, 0}}];
    
    next[prev_] := prev /. Polygon[{p1_, p2_, p3_}] :> {
         Polygon[{p1, (p1 + p2)/2, (p1 + p3)/2}],
         Polygon[{p2, (p2 + p3)/2, (p1 + p2)/2}],
         Polygon[{p3, (p1 + p3)/2, (p2 + p3)/2}]};
    
    draw[n_] := Module[{edges},
       edges = Flatten@Nest[next, N@axiom, n] /. Polygon[pts_] :>
          Sequence @@ UndirectedEdge @@@ Partition[pts, 2, 1, 1];
    
       Graph[edges, VertexCoordinates -> VertexList[Graph[edges]],
        VertexSize -> .25]];
    
    style = Sequence[EdgeStyle -> Orange];
    circles[g_Graph, r_: 1] := (circles[g, r] =
        Module[{vs = VertexList[g]},
         DeleteDuplicates[
          NeighborhoodGraph[g, #, r, style] & /@ vs,
          IsomorphicGraphQ]]);
    
    Pane[#, 600] &@Column[
      Row[Prepend[
          circles[draw[5], #],
          Style[#, Lighter[Black, 1/6]]], " "] & /@ Range[1, 3],
      Alignment -> Center, Spacings -> 1]
    

This brings to light a more important question: What the hell is this? It's like the ugly duckling of radius 3 Sierpinski subgraphs. Just look at it. LOL. But OK, I mightn't myself be the most handsomest chap on the block, and graphs are people too after all.

There's a good chance that subgraph is hideous because it contains one of the 3 end vertices of the graph as a whole, though I'm too lazy to check this. Those vertices are in part pathological because they have degree 2, whereas all the other vertices have degree 4. But really I think the Sierpinski graph itself is contrived. At least, the finite version seems contrived to me.

Perhaps because the Sierpinski pattern might actually be a grid, in the sense that the empty space is an integral part of its characterization a la our infinite quadrilateral descent construction. If we base a graph on the pattern produced by the mod 2 binomial, we get this graph:

  1. locs[n_] := RotationTransform[-Pi/2] /@
       Position[Mod[Array[Binomial, {n, n}, 0], 2], 1];
    
    draw[n_, grid_: True] :=
      Module[{edges, subsets = Subsets[locs[n], {2}]},
       edges = UndirectedEdge @@@
         Pick[subsets,
          ManhattanDistance(*;ChessboardDistance*)@@@ subsets, 1];
    
       Graph[edges, VertexCoordinates -> If[grid,
          {.5, .5} + # & /@ VertexList[Graph[edges]]]]];
    
    Show[
     Graphics[Rectangle /@ locs[2^4]],
     draw[2^4]]
    

Which looks like this in a tiered layout. Maybe the "real" Sierpinski graph is a binary tree of this sort, and it's only connected on all sides in the infinite case. And maybe right now I'm making mathematicians bash their heads against walls, which would be awesome.

The binomial mod 2 construction was one of the approaches that went AWOL during our 3Dification blitz. Does it have a 3 dimensional version? Yes, the multinomial mod 2. The code is almost as pretty as it is for the 2D version:

  1. array[n_] := Mod[Array[Multinomial, n {1, 1, 1}, 0], 2];
    
    draw[n_] := Graphics3D[Cuboid /@ Position[array[n], 1],
       Lighting -> "Neutral", Boxed -> False];
    
  2. array[n_] := Mod[Array[Multinomial, n {1, 1, 1}, 0], 2];
    
    draw[n_] := Module[{edges, subsets},
       subsets = Subsets[Position[array[n], 1], {2}];
    
       edges = UndirectedEdge @@@
         Pick[subsets,
          ManhattanDistance(*;ChessboardDistance*)@@@ subsets, 1];
    
       Graph[edges]];
    
    draw[2^3]
    GraphPlot3D[draw[2^3], PlotStyle -> ColorData[1][1]]
    

Poor Boxed, always being set to False. What happens to our chaos game algorithm if we implement some notion of momentum for the active point?

  1. game[v_, numPoints_] := Module[{vertices, update, vl},
       vertices = {Cos[#], Sin[#]} & /@ (2 Pi Range[v]/v);
    
       update[{pos0_, vl0_}, nextVertex_] := (
         (*vl=vl0+(pos0+nextVertex)/650-pos0;*)
         (*vl = vl0+Normalize[nextVertex - pos0];*)
         vl = vl0 + .0001 (nextVertex - pos0);
         {pos0 + vl, vl});
    
       First /@ FoldList[update, N@{{0, 0}, {0, 0}},
         RandomChoice[N@vertices, numPoints]]];
    
    Graphics[{PointSize[0], Opacity[.1],
      Point[game[3, 100000]]}]
    
  2. game[v_, numPoints_] := Module[{vertices, update, vl},
       vertices = {Cos[#], Sin[#]} & /@ (2 Pi Range[v]/v);
    
       update[{pos0_, vl0_}, nextVertex_] := (
         vl = vl0 + Clip[(nextVertex - pos0), 20 {-1, 1}];
         {pos0 + vl, vl});
    
       First /@ FoldList[update, N@{{0, 0}, {0, 0}},
         RandomChoice[N@vertices, numPoints]]];
    
    zrp = ParallelTable[Module[{pts, max, maxx, maxy},
        pts = game[4, n];
        max = Sqrt[2] Max[Abs /@ pts];
    
        SeedRandom[400, Method -> "ExtendedCA"];
        Rasterize@Graphics[{Opacity[.9], Line[pts]},
          PlotRange -> max, ImageSize -> {150, 150}]],
       {n, 1, 10000, 83}];
    
    Export["c:/users/zrp/desktop/zrp.gif",
     ColorQuantize[#, 4] & /@ Most[zrp]]
    

I didn't find any interesting formulas, but still I managed to get a variety of figures by fiddling with numbers. Probably I would have to use math to find something more interesting.

The figures have precise symmetries (180 degree rotation), apparently because the particle eventually overshoots far enough that the randomness becomes a small jitter component (because the vertices of the game become very distant), so it accumulates a near-linear velocity/path on its way back. I'm not sure about this explanation though.

Even the ones that look like random walks are symmetric. They aren't standard random walks, rather the particle is overshooting back and forth. This raises the idea of symmetrizing random walks:

  1. draw[v_, r_, numSteps_: 100] := Module[{directions, walk},
       directions = {Cos[#], Sin[#]} & /@ (2. Pi Range[v]/v);
    
       walk = Accumulate[RandomChoice[directions, numSteps + 1]];
       Graphics[Rotate[Line[walk], #] & /@ (2. Pi Range[r]/r)]];
    

Awesome possum. What do you see here? When fiddling with momentum I found a simple variation on the logarithmic distance function:

  1. Module[{gameC},
      gameC = Compile[{{spec, _Real, 1}, {numPoints, _Integer, 0}},
        Module[{vertices, diff, v, w, s, b, r, p, z, c},
         {v, w, s, b, r, p, z, c} = spec;
         vertices = {Cos[#], Sin[#]} & /@ (2 Pi Range[v]/v);
    
         FoldList[(
            diff = #2 - #1;
            p #1 + Clip[(#1 + z #2) Log[b, (diff.diff)^(1/r) + w], c {-1, 1}]) &,
          {0, 0}, RandomChoice[s vertices, numPoints]]]];
    
      game[spec_, numPoints_] :=
       gameC[PadRight[spec, 8, {3, .5, 1, E, 2, 0, 1, 2}], numPoints]
      ];
    
    draw[spec_, numPoints_: 100000, rot_: 0, style___] := Graphics[{
        PointSize[0], Opacity[.1], style,
        Rotate[Point[game[spec, numPoints]], rot]}];
    
    draw[{5, .8}]
    
    (*{vertex count, w factor, scale, base, root, prefix, polarity, clip range}*)
    (*Show[draw[{3, .12, .3, 3.9, .9, 1, -.5 E, 1}, 2000000, -Pi/6],
       ImageSize -> 2 1280] // Rasterize // ImageResize[#, Scaled[1/4]] &*)
    

Or something of a generalization. I blindly parameterized several parts of the formula. Some of the parameters are sensitive, but in any case it's easy to find spiffy images. The hard part is deciding which of them to put here. 3D version:

  1. Module[{gameC},
      gameC =
       Compile[{{vertices, _Real, 2}, {spec, _Real, 1}, {numPoints, _Integer, 0}},
        Module[{diff, w, r, s, b, p, z, c},
         {w, s, b, r, p, z, c} = spec;
    
         FoldList[(
            diff = #2 - #1;
            p #1 + Clip[(#1 + z #2) Log[b, (diff.diff)^(1/r) + w], c {-1, 1}]) &,
          Mean[vertices] RandomReal[], RandomChoice[s vertices, numPoints]]]];
    
      game[vertices_, spec_, numPoints_: 100000] :=
       gameC[vertices, PadRight[spec, 7, {.5, 1, E, 2, 0, 1, 2}], numPoints]
      ];
    
    draw[vertices_, spec_, numPoints_: 100000, options___] :=
      Graphics3D[{PointSize[0], Opacity[.1],
        Point[game[vertices, spec, numPoints]]},
       options, Boxed -> False];
    
    (*{w factor, scale, base, root, prefix, polarity, clip range}*)
    
    vertices = PolyhedronData[{"Dipyramid", 5}, "VertexCoordinates"];
    draw[vertices, {.8}, 200000]
    
    (*note.*)
    v2D = {Cos[#], Sin[#]} & /@ (2 Pi Range[5]/5);
    {Graphics[Point[game[v2D, {.8}]]],
     Graphics3D[Point[game[{##, 0} & @@@ v2D, {.8}]]]}
    

I love how some of these look like sketches. You'd expect to find this as an illustration in a wizard's journal, but it's actually from a Graphics3D pane in my Mathematica notebook. This opportune box, beside being the final confine of a truculent force, is the result of the clipping I use to keep the point from escaping. I set the clipping as a parameter because it can be used to effect.

This clipping restriction isn't always necessary, and it might not be necessary for all points within a given figure, which raises an interesting prospect: What if we try to identify the points that fly off into infinity and those that don't?

    1. check = Compile[{{c, _Complex, 0}},
         Module[{i = 0, z = 0 I},
          While[
           Abs[z] < Sqrt[2] && i++ < 240,
           z = z^2 + c];
          -i]];
      
      ImageAdjust@Image@
        ParallelTable[check[x + y I],
         {y, -1.1, 1.1, .0035}, {x, -1.55, .6, .0035}]
      
  1. vertices = (**)5(**) {Cos[#], Sin[#]} & /@ (2 Pi Range[3]/3);
    
    check = Compile[{{x, _Real, 0}, {y, _Real, 0}},
       Module[{i, b, diff, z = {0., 0.}, vertices = vertices},
        Total@Table[
          i = 0; z = {x, y};
          While[z.z < 40 && i++ < 120,
           b = RandomChoice[vertices];
           diff = b - z;
           z = (z + b) Log[Sqrt[diff.diff] + .01]];
          -i,
          {20(*0*)(*number of trials*)}]]];
    
    img = ImageAdjust@Image@
        ParallelTable[check[x, y],
         {y, -6.5, 6.5, .01}, {x, -6.5, 6.5, .01}];
    
    img // Colorize // ImageResize[#, 550] &
    

Here the white points go off into infinity quickly. The black points don't, or at least they take a lot longer to escape. There are certainly patterns here, but they're much less pronounced and computationally harder to reveal to than they are for the Mandelbrot set, which is the same idea for Julia iterations. But if you spin a few knobs you can find interesting figures irregardless. The different colors/shades are different escape speeds. It may not be immediately apparent, but these are fractals also.

A lot of fractals have scaled/skewed characteristics, including the Mandelbrot set. I wonder if there's a non-trivial chaos game that can create the Mandelbrot set. Since we're skittering around complex numbers, is there an interesting complex-valued version of the logarithmic chaos game?

  1. game = Compile[{{v, _Integer, 0}, {numPoints, _Integer, 0}},
       Module[{vertices},
        vertices =(*1.5*)E^(I 2 Pi Range[v]/v);
    
        FoldList[(*(Log[#1]+#2)/2&*)
         (#1 + (#1 + #2) Log[Sqrt[#2 - #1] + .7])/2.1 &, .1,
         RandomChoice[N@vertices, numPoints]]]];
    
    Graphics[{Opacity[.1], PointSize[0],
      Point[{Im[#], Re[#]} & /@ game[2, 400000]]}]
    

I don't know. Strictly speaking there wouldn't be a difference, but if you put on a blindfold and chuck logarithms at Mathematica helter-skelter, pretty pictures eventually come out. So I guess the answer is some form of yes. Formally these might still be considered Julia sets.

I've mentioned before that there are a lot of crazy distance functions out there for us to use in our chaos game, and there isn't anything special about the logarithm function. What do plots using other functions look like?

    1. vertices = (**)5(**) {Cos[#], Sin[#]} & /@ (2 Pi Range[3]/3);
      
      check = Compile[{{x, _Real, 0}, {y, _Real, 0}},
         Module[{i, b, diff, z = {0., 0.}, vertices = vertices},
          Total@Table[
            i = 0; z = {x, y};
            While[z.z < 40 && i++ < 120,
             b = RandomChoice[vertices];
             diff = b - z;
             z = (z + b) Sin[Sqrt[diff.diff] + .01]];
            -i,
            {20(*0*)(*number of trials*)}]]];
      
      img = ImageAdjust@Image@
          ParallelTable[check[x, y],
           {y, -6.5, 6.5, .01}, {x, -6.5, 6.5, .01}];
      
      img // Colorize // ImageAdjust // ColorConvert[#, "Grayscale"] & // ImageAdjust //
          ImageResize[#, Scaled[1/2]] & // ImageRotate[#, -Pi/2] & // ColorNegate //
       ImageApply[#^(1/1.3) &, #] &
  1. game[f_, rest__] := rest //
       Compile[{{v, _Integer}, {w, _Real}, {numPoints, _Integer}, {rot, _Real}},
        Module[{diff, tmp, vertices},
         vertices = {Cos[#], Sin[#]} & /@ (rot + 2 Pi Range[v]/v);
    
         FoldList[(
            diff = #2 - #1;
            tmp = f[Sqrt[diff.diff] + w];
            Clip[(#1 + #2) tmp, 24 {-1, 1}]) &,
          {0, 0}, RandomChoice[vertices, numPoints]]]];
    
    draw[{args__}, rot_: 0, options___] :=
      Graphics[{PointSize[0], Opacity[.25],
        Point[game[args, rot]]}, options];
    
    draw[{Sin, 5, -1.23, 100000}, -Pi/10]
    draw[{Cos, 3, Pi/2.675, 100000}, -Pi/6]
    draw[{RamanujanTauL, 5, 3.1, 10000}, -Pi/10]
    

They look just as awesome, of course. Here we have plots using the sine, cosine, and, you guessed it, Ramanujan tau Dirichlet L-function. And this is using the same basic form as the logarithm version, without us even having to put on Loki's mask and get real buckwild. Speaking of masks.

Usually I don't pick favorites, but I like the cosine image (of course sin/cos are just offsets of eachother) because it has an infinite number of folded sheet things that seem to have precise contours. I'll leave the 3Dification as an exercise, but not the how-fast-points-go-to-infinity plot.

The reason it's easy to get all these pictures without trying very hard is that the self-similarity is almost guaranteed by the chaos game algorithm. As we saw earlier, "move toward a point" amounts to the same thing as "make a resized copy of everything toward the perspective of that point."

This is a simplification, but the point is that you essentially get the skeleton of self-similarity for free, or perhaps something a bit more broad. And more abstractly, I think some remarks could be made about the real number system itself.

What does the Sierpinski triangle sound like? One easy interpretation is to consider the L-system construction for the triangle and convert different angles to different frequencies as the turtle makes the triangle:



  1. mp3  midi


  2. mp3  midi
  3. axiom = A;
    rules = {A -> {B, R, A, R, B}, B -> {A, L, B, L, A}};
    conversions = {A -> forward, B -> forward, L -> left, R -> right};
    
    (*state transformations*)
    forward[{z_, theta_}] := {z + E^(I theta), theta};
    left[{z_, theta_}] := {z, theta + 2. Pi/6};
    right[{z_, theta_}] := {z, theta - 2. Pi/6};
    
    sier[n_] := Module[{program, zs},
       program = Flatten[Nest[# /. rules &, axiom, n]] /. conversions;
       zs = First /@ ComposeList[program, {0, 0}];
       First /@ Split[{Re[#], Im[#]} & /@ zs]];
    
    (*convert angle into the given frequency range*)
    freq[min_, max_][angle_] := angle (max - min)/Pi + min;
    
    wave[coords_, dur_: 10, freq_: freq[6, 30]] := Module[{angles, freqs},
       angles = Abs[ArcTan @@@ Differences[coords]];
       freqs = Round /@ freq /@ angles;
       Sound[SoundNote /@ freqs, dur]];
    
    wave[sier[3], 5, freq[8, 15]]
    
    (*overtones zomg*)
    Sound[{"NewAge",
      wave[#, {0, 20}] & /@ Table[
        RotationTransform[2 Pi i/4] /@ sier[i], {i, 4}]},
     SoundVolume -> .8]
    

It sounds totally lame. Not surprising since the L-system construction is simple. There is real power here though. This tonifier operates on coordinate lists of any kind, not just those produced by this particular L-system. And if you do things like layer different iterations on top of each other, you can get nifty chord thingies.

A variation of this would be to determine the waveform directly from the L-system. In a past life I made such a program in C#/WPF. It was around 1200 lines of code. In Mathematica it would be around 30 lines of code, and maybe around 150 lines total with a solid UI around it. It would also be about a million times more powerful/general/flexible. There's a lot of reasons for this, none of which have to do with math.

Luckily for me I don't have to explain. The goddess of finishing projects has finally crawled out of her cave and seen fit to smite her lightning bolt through my ears and across my temporal lobe, for this tune has sated and sedated the voices and quelled their cantankerous echoes. And so ends part 1.



This page made while drinking Starbucks and listening to CoLD SToRAGE. If you know programming, consider contributing to Mathics. If you're having trouble with the code snippets, try clearing all variables or restarting the kernel. If you're losing the formatting of copy-pasted code, and that annoys you, right click -> insert code cell. For general Mathematica inquiries, visit Mathematica.SE or the Wolfram Community.

With the exception of Mr. Scruples who is under "CC BY NC SA", some of his companion source code stolen from StackExchange under "CC BY SA", and the geometry of our Brahman data cow whose license status is unknown, all images/animations/video/audio/source code on this page is in the public domain. It would be the best thing ever if you made money from my work. :D Also see this post by Vitaliy Kaurov for some fun info.