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: 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): axiom = {A}; rules = {A -> {B, R, A, R, B}, B -> {A, L, B, L, A}}; conversions = {A -> forward, B -> forward, L -> left, R -> right}; 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: 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: 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: Pascal's triangle mod 4 Pascal's triangle $x\equiv 2\space(\text{mod }4)$ 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: 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: 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: 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: $$ \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) $$ 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): 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)$$ $$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)$$ $$\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: $$ \begin{array}{cc} \binom{0}{0}=1 & \binom{0}{1}=0 \\ \binom{1}{0}=1 & \binom{1}{1}=1 \\ \end{array} $$ 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)!}$$ $$\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: $

eg(

eg n\land k)=n\lor

eg 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. 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: 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

eg 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 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]]]}]; 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: 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}]; 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: triangle :> Polygon triangle :> Line triangle [ pts_ ] :> Line [ RandomChoice [pts, RandomInteger [{ 2 , 3 }]]] triangle [ pts_ ] :> Disk [ Mean [pts], 1 / 2 ^ (n + 1 )] 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) $$

$$ \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. kitty = WolframAlpha [ "cat picture" , "PodImages" ][[ 2 ]]; 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 $$ $$ 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!}$$ $$\frac {6!} {4!2!}$$ From the factorial definition of the binomial coefficient: $$\binom{n}{k} = \frac {n!} {k!(n-k)!}$$ $$\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: 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: 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: polyGraph [ poly_ , options___ ] := Graph [ UndirectedEdge @@@ PolyhedronData [poly, "EdgeIndices" ], options]; sierpinskiMatrixGraph [ n_ , options___ ] := Module [{a}, "inner" "external" 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: $$ \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} $$ 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: $$ \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) $$ 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: 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: 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]}]; 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 "SpringEmbedding" ]]; 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: hanoiGraph [ state [{}, {}, {}, Range [ 4 ]]] hanoiGraph [ state [{}, {}, {}, {}, Range [ 4 ]]] hanoiGraph [ state [{}, {}, Range [ 3 ], Range [ 3 ]]] hanoiGraph [ state [{}, {}, { 1 }, Range [ 3 ]]] hanoiGraph [ state [{}, {}, { 2 }, Range [ 3 ]]] hanoiGraph [ state [{}, {}, { 3 }, Range [ 3 ]]] validQ [ s_state ] := And @@ Equal @@@ s; hanoiGraph [ state [{}, {}, ConstantArray [ 1 , 5 ]]] validQ [ s_state ] := And @@ LessEqual @@@ s; hanoiGraph [ state [{}, ConstantArray [ 2 , 3 ], ConstantArray [ 1 , 3 ]]] 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: 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_ , options___ ] := Module [{vertices, edges, n}, n = Count [s, _Integer , Infinity ]; {vertices, {edges}} = Reap [ Nest [neighbors, {s}, 2 ^ n]]; SetAttributes [UndirectedEdge, Orderless ]; Graph [ DeleteDuplicates [edges]]]; 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 ] }]]; {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}: 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 ]; 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 }]}; "DynamicSetting" 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 ], 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: 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: 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 ] 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 ]]]; 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: 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 }]}; 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). 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: 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 ]] 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: 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: 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: 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: 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]}; 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 ]] 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]}; 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 / 1 0 . }]; 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: 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: game = Compile [{{v, _Integer }, {wowzerz, _Real }, {numPoints, _Integer }}, Module [{diff, vertices}, vertices = { Cos [ # ], Sin [ # ]} & /@ ( 2 Pi Range [v] / v); FoldList [( diff = #2 - #1 ; 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 ]]} ] 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: 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] 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 [ 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: 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:





chars = Characters [ GenomeData [{ "ChromosomeX" , { 28000000 , 36000000 }}]]; tallies = Tally [chars]; BarChart [ Last /@ tallies, ChartLabels -> First /@ tallies] 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 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: 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)$$ $$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$$ $$p_{i-1}=2p_i-v$$ So, reversing all the points in a particular subquadrant: 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: 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: 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: 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:





ListPlot [ FoldList [( #1 + #2 ) / 2 & , 1 , Mod [ Range [ 1000 ], 2 ]], PlotRange -> 1 , Ticks -> { Automatic , Range [ 0 , 1 , 1 / 3 ]}] 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 ] ~ 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: chars = Characters [ GenomeData [{ "ChromosomeX" , { 28000000 , 36000000 }}]]; tallies = Sort [ Tally [ Partition [chars, 3 , 1 ]]]; BarChart [ Last /@ tallies, ChartLabels -> Column /@ First /@ tallies] 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: $$ \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} $$ 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 ] & 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: 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: 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: 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: 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: 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 ] matrixInput [ Dynamic [ m_ ], Dynamic [ rot_ ]] := Dynamic [ Rotate [ Deploy [ MatrixForm [ # , TableSpacing -> { 0 , 0 }]], rot] & @ Array [ 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}]}]]; 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], {{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: 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 ] matrixInput1 [ Dynamic [ m_ ], Dynamic [ rot_ ]] := Dynamic [ Rotate [ Deploy [ MatrixForm [ # , TableSpacing -> { 0 , 0 }]], rot] & @ Array [ 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

Right-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]; 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)}; 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}]]; 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]]], {{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$$ $$\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}$$ 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]] $$z_\text{inv} = \frac 1 {|z|}\frac z {|z|} = \frac z {|z|^2}$$ 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: 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: 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|}$$ $$\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|}$$ $$\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: 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_ ] := p / Norm [p] ^ 2 ; Column [ Panel [ Row [ # ]] & /@ Table [ draw [v, n] /. Polygon [ pts_ ] :> Polygon [ invert /@ pts], {v, 3 , 6 }, {n, 0 , 4 }]] 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: invert [ p_ ] := p / Norm [p] ^ 2 ; 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]]]; drawInv [ v_ , n_ ] := Module [{g = draw [v, n]}, Show [g /. poly_Polygon :> invertPoly [poly], PlotRange -> 1.1 plotRangeInv [g]]]; 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}}}$$ $$\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: 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: 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: 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: 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]]; 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; transformationFunctions = { # / d [ # ] ^ 2 & , # ^ 3 / d [ # ] ^ 2 & , ( Reverse [ # ] . # ) # / d [ # ] & , ( # + # ^ 3 / d [ # ] ^ 2 ) / 2 & , ( # + # / d [ # ] ^ 2 ) / 2 & , Round [ # , . 05 ] &