Generating a Sierpinski carpet

Version 11.1 introduces MengerMesh:

MengerMesh[3]

---

---

This seems the most natural to me:

carpet[n_] := Nest[ArrayFlatten[{{#, #, #}, {#, 0, #}, {#, #, #}}] &, 1, n]

ArrayPlot[carpet @ 5, PixelConstrained -> 1]

Shorter (in InputForm), but perhaps harder to read and slightly slower, though speed hardly matters given the geometric memory usage:

carpet[n_] := Nest[ArrayFlatten @ ArrayPad[{{0}}, 1, {{#}}] &, 1, n]

---

Style by level

With a minor change we can increment the values with each fractal level allowing identification such as styling, or other processing.

Wild colors are but a few commands away:

carpet2[n_] := Nest[ArrayFlatten[{{#, #, #}, {#, 0, #}, {#, #, #}}] &[1 + #] &, 1, n]

Table[
  ArrayPlot[carpet2 @ 4, PixelConstrained -> 1, ColorFunction -> color],
  {color, ColorData["Gradients"]}
]

---

Extension to three dimensions

A Menger sponge courtesy of chyanog, with refinements:

carpet3D[n_] :=
   With[{m = # (1 - CrossMatrix[{1,1,1}])}, Nest[ArrayFlatten[m, 3] &, 1, n]]

Image3D[ carpet3D[4] ]

---

Element coordinates

If you wish to get coordinates for display with graphics primitives or analysis this can be done efficiently using SparseArray Properties:

coords = SparseArray[#]["NonzeroPositions"] &;

Example usages:

Graphics @ Point @ coords @ carpet @ 4

Graphics3D[Cuboid /@ coords @ carpet3D @ 3]

Assuming you want a vector-based image, it's more efficient to cut holes:

translations = {#, #} & /@ Complement[Tuples[{-1, 0, 1}, 2], {{0, 0}}];
shrink[{{x0_, y0_}, {x1_, y1_}}] := {{2 x0 + x1, 2 y0 + y1}, {x0 + 2 x1, y0 + 2 y1}}/3
children[sq : {{x0_, y0_}, {x1_, y1_}}] := 
  With[{side = x1 - x0, newsq = shrink[sq]},
    (newsq + #) & /@ (side translations)]
gen = NestList[Join @@ children /@ # &, {{{1, 1}, {2, 2}}/3}, 3];
Graphics[Rectangle @@@ Join @@ gen]

For the same generation you need to draw $1/8$ of the rectangles.

ngen = 4;
gen = NestList[Join @@ children /@ # &, {{{1, 1}, {2, 2}}/3}, ngen];
colors = Table[Blend[{RGBColor[0.5, 0.5, 1], White}, i/ngen], {i, 0, ngen}];
Graphics[{Black, Rectangle[{0, 0}, {1, 1}], White, Rectangle @@@ Join @@ gen}]
Graphics[{Red, Rectangle[{0, 0}, {1, 1}], MapThread[Prepend, {Apply[Rectangle, gen, {2}], colors}]}]

Here is the 3D version:

rule = 0 -> CrossMatrix[{1, 1, 1}];
Graphics3D[Cuboid /@ Position[Nest[ArrayFlatten[# /. rule, 3] &, 0, 3], 0]]

Here are two methods using rules, shamelessly modified from a MathGroup posting (http://forums.wolfram.com/mathgroup/archive/2007/May/msg01356.html).

rules = # -> ArrayPad[{{0}}, 1, #] & /@ {0, 1}
f1[m_] := ArrayFlatten[m /. rules]
drawSerp[n_] := ArrayPlot[Nest[f1, 1, n], Frame -> False]
drawSerp[3]

An alternative cute ASCII implementation, borrowed from http://rosettacode.org/wiki/Sierpinski_carpet#Mathematica

n = 3;
Grid[Nest[ArrayFlatten[# /. rules] &, {{1}}, 
   n] //. {0 -> " ", 1 -> "#"}]